IIR digital filter: samples not matching the frequency response

Question

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"} ]

Answer 1

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.