IIR digital filter: samples not matching the frequency response

I am working on a small project in Mathematica involving a Chebyshev 4th order highpass filter, trying to obtain its digital IIR filter with the bilinear transformation.

I implemented the differential equation of the filter and calculated the coefficients of the impulse response.

To verify if the implementation is correct, I applied the inverse Fourier transform and then plotted the gain and the points over it.

I've used this code for other filters and no problems appeared.

Here is the example I did for a bandpass filter the other day:

As we can see the samples match the gain obtained with the frequency response perfectly.

I have no idea what is going on with this specific case as the first 3 samples are off:

At first I thought "well maybe because this is a highpass filter there is some problem with the convergence of the the inverse discrete fourier Transform or something like that," but then I tried a Butterworth highpass filter, similar code adapted, everything was fine:

There is some detail that I am missing but I can't find it. Can someone help me please?

Below is the full workable code in Mathematica:

    ClearAll["Global*"]
FSamp = 22050;
(* PROJETO 1 CHEBYSHEV PASSA-ALTO ORDEM 4 FP= 2KHZ AP= 2 DB *)
TCheb = A0/(B0 + B1 S + B2 S^2 + B3 S^3 + S^4);
DCheb = Collect[(S^2 + 0.50644045 S + 0.22156843) (S^2 +
0.20977450 S + 0.92867521), S];
K = 0.16344503;
A0 = K;
aux = CoefficientList[DCheb, S];
B0 = aux[[1]];
B1 = aux[[2]];
B2 = aux[[3]];
B3 = aux[[4]];
TCheb
TChebDenorm = TCheb /. S -> wp/s;
wp = 2 \[Pi] 2000;
wp = 2 \[Pi]*FSamp/\[Pi]*
Tan[(\[Pi]*wp/(2 \[Pi]))/FSamp];(*pre-distorção da frequência*)
DenTFilter = Simplify[Denominator[TChebDenorm]*s^4]
NumTFilter = Simplify[Numerator[TChebDenorm]*s^4]
aux = CoefficientList[DenTFilter, s];
DenTFilter = Simplify[DenTFilter/aux[[5]]]
NumTFilter = s^4
TFilter1 = NumTFilter/DenTFilter
aux = FactorList[Denominator[TFilter1]];
DenTFilter11 = aux[[2]][[1]];
DenTFilter12 = aux[[3]][[1]];
TFilter11 = s^2/DenTFilter11
TFilter12 = s^2/DenTFilter12
TDigital11 = Simplify[TFilter11 /. s -> (2*FSamp*(z - 1)/(z + 1))];
Den = Simplify[Denominator[TDigital11]/z^2];
Num = Expand[Numerator[TDigital11]/z^2];
TDigital11 = Num/Den
TDigital12 = Simplify[TFilter12 /. s -> (2*FSamp*(z - 1)/(z + 1))];
Den = Simplify[Denominator[TDigital12]/z^2];
Num = Expand[Numerator[TDigital12]/z^2];
TDigital12 = Num/Den
TDigital1 = TDigital11*TDigital12;
graphicGAINDigital1 =
Plot[20 Log[10, Abs[TDigital1 /. z -> E^(I 2 \[Pi] f/FSamp)]], {f,
1, FSamp/2}, PlotRange -> {{1, FSamp/2}, {-200, 10}},
PlotPoints -> 1000,
PlotStyle -> {{AbsoluteThickness[0.5], RGBColor[0, 0, 1]}},
GridLines -> Automatic, Frame -> True,
BaseStyle -> {FontFamily -> "Times", FontSize -> 12} ,
PlotLegends -> {"|T(\!$$\*SuperscriptBox[\(e$$, \
$$j\[Omega]$$]\))|"}];
NSamples = 128;
Clear[c01, c11, c21, d11, d21, c02, c12, c22, d12, d22];
Table[input[n] = 0, {n, 0, NSamples - 1}];
Table[outputW1[n] = 0, {n, 0, NSamples - 1}];
Table[outputY1[n] = 0, {n, 0, NSamples - 1}];
Table[outputW2[n] = 0, {n, 0, NSamples - 1}];
Table[outputY2[n] = 0, {n, 0, NSamples - 1}];
Table[GainDFT[n] = 0, {n, 0, NSamples - 1}];
TDigital11;
aux = CoefficientList[Numerator[TDigital11]*z^2, z];
c01 = aux[[3]];
c11 = aux[[2]];
c21 = aux[[1]];
aux = CoefficientList[Denominator[TDigital11]*z^2, z];
d11 = -aux[[2]];
d21 = -aux[[1]];
TDigital12;
aux = CoefficientList[Numerator[TDigital12]*z^2, z];
c02 = aux[[3]];
c12 = aux[[2]];
c22 = aux[[1]];
aux = CoefficientList[Denominator[TDigital12]*z^2, z];
d12 = -aux[[2]];
d22 = -aux[[1]] ;
input[0] = 1;
Table[input[n] = 0, {n, 1, NSamples - 1}];
outputW1[0] = input[0];
outputY1[0] = c01*outputW1[0];
outputW2[0] = outputY1[0];
outputY2[0] = c02*outputW2[0];
outputW1[1] = input[1] + d11*outputW1[0];
outputY1[1] = c01*outputW1[1] + c11*outputW1[0];
outputW2[1] = outputY1[1] + d12*outputW2[0];
outputY2[1] = c02*outputW2[1] + c12*outputW2[0];
Table[outputW1[n] =
input[n] + d11*outputW1[n - 1] + d21*outputW1[n - 2], {n, 2,
NSamples - 1}];
Table[outputY1[n] =
c01*outputW1[n] + c11*outputW1[n - 1] + c21*outputW1[n - 2], {n,
2, NSamples - 1}];
Table[outputW2[n] =
outputY1[n] + d12*outputW2[n - 1] + d22*outputW2[n - 2], {n, 2,
NSamples - 1}];
Table[outputY2[n] =
c02*outputW2[n] + c12*outputW2[n - 1] + c22*outputW2[n - 2], {n,
2, NSamples - 1}];
Table[GainDFT[n] =
20*Log[10,
Abs[Sum[outputY2[m]*E^(-I*m*n* ((2 \[Pi])/NSamples) ), {m, 0,
NSamples - 1}]]], {n, 0, NSamples - 1}]
GainDFTPlot =
Table[{n/(NSamples/2)*FSamp/2, GainDFT[n]}, {n, 0, NSamples/2}];
GainDFTGRAPH =
ListPlot[GainDFTPlot,
PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.010]},
Filling -> Axis, GridLines -> Automatic,
FrameLabel -> {"n", "h(n)"}, Frame -> True, PlotRange -> All,
ImageSize -> Medium];
Show[graphicGAINDigital1, GainDFTGRAPH,
FrameLabel -> {"f(Hz)", "dB"} ]

• Since the filter has infinite samples in the impulse response and ifft requires finite number of samples, how many samples of the impulse response did you use for ifft? Can you show a result with double the number of samples? My intuition is that, for high pass filter with infinite attenuation at DC frequency, the positive samples in the impulse response must exactly cancel the negative samples to represent infinite attenuation. Please report back.
– AJN
Aug 16, 2021 at 0:58
• Probably try NSamples = 1024`.
– AJN
Aug 16, 2021 at 1:00
• I've tried 256 samples and it already looked better. I will also try 1024. Thank you for your explanation, can you please elaborate more on why the infinite attenuation at low frequency does that cancellation of terms? Aug 16, 2021 at 2:43
• Since it seems to improve the result, i will write it up as an answer.
– AJN
Aug 16, 2021 at 12:51

TLDR

It is the effect of truncation of an IIR filter to an FIR filter.

Longer explanation

For a high pass filter with zero gain at DC, the average value of the impulse response is exactly zero. This is by definition. For an IIR (discrete-time) filter, since the impulse response is infinitely long, the average value may be reach zero only after infinite samples.

If the impulse response is truncated arbitrarily, the average value need not be zero for the chosen number of samples used for truncation. More over, truncation also effectively converted the filter to FIR whose response need not perfectly match the original IIR frequency response. When this truncated impulse response is FFT'd to recreate the frequency response you are actually viewing the frequency response of the FIR filter obtained via truncation of the infinite impulse response.

Below is an example I constructed.

The high pass filter frequency and impulse response are shown in the first two plots. The running average of the impulse response is shown in the third plot (It is actually the step response).

The points were the running average crosses zero result in the truncated response having lower gain at zero frequency. In the example shown they are 12, 32 and 192 samples.

The point were the running average has a peak shows that, the response if truncated there, will have a significant DC component. In the example it is 22 samples.

The bottom 4 plots show the corresponding FFT'd frequency responses superimposed over the original response.

note

Even though the zero frequency / DC value is easiest to intuitively understand, other frequencies are also clearly affected. It is clearly seen in the "22 samples" example.

• In short, if the DC gain is important, use the running average to find the ideal number of samples to truncate at. Unfortunately it will may not be a power of 2. Pad with zeros to the next power of two. Zero padding won't shift the response average value away from zero. Alternately, you can find the average of the truncation and add a single additional point manually to the truncated response which will make the new average zero.
– AJN
Aug 16, 2021 at 14:02
• Thank you for such a detailed explanation, I totally got it. One of the experiments I did with one of the filters I was analyzing (another highpass filter with zero DC gain) always had the first sample with the incorrect phase even with 1024 samples. But now I get why: all samples sum to zero but not exactly zero, and in that case it sums to a negative value very close to zero -> phase comes at 180 degrees instead of zero degrees. Aug 16, 2021 at 14:31