# IIR Digital Filter: samples not matching the frequency response

So I am a small project in Mathematica involving a Chebyshev 4th order highpass filter, and to obtain its digital IIR filter with the bilinear transformation. Then I implement the differences equation of the filter and calculate the coefficients of the impulse response. To verify if the implementation is correct, I apply the inverse Fourier transform and then plot 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. However 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 pull it off. 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"} ]


Your code is long and rather difficult for me to sort through. I thought I would give you the Mathematica way of generating filters. Hope this helps.

First I choose a Chebyshev filter from those provided. I then choose a plot range and inspect.

filt = Chebyshev1FilterModel[{"Highpass", 4, 2 \[Pi] 2000.}];
pr = {{0, 11000}, {10^-4, 10}}; (* plot range *)
LogPlot[Abs[filt[I 2 \[Pi] f]], {f, 0, 11000}, PlotRange -> pr,
Frame -> True]


Next we apply the bilinear transform (the default option, others are available) and get a discrete filter.

FSamp = 22050;
disFilt = ToDiscreteTimeModel[filt, 1/FSamp];


Now to check the filter. I generate a time history which is all zeros except the first point which is 1.0. Then I apply the filter to get the impulse response function. The plot looks sensible.

nn = 2^14;
th = ConstantArray[0, nn];
th[[1]] = 1;
ir = RecurrenceFilter[disFilt, th];
ListLinePlot[ir[[1 ;; 200]], PlotRange -> All, Frame -> True,
Mesh -> All]


Now I take the Fourier transform of the impulse response function and plot that, the discrete filter and the original filter function. Here ff are the frequency values

ft = Fourier[ir, FourierParameters -> {1, -1}];
ff = Table[(n - 1) FSamp/nn, {n, 2^14}];
pr = {{0, 11000}, {10^-4, 10}};
Show[
ListLogPlot[Transpose[{ff, Abs[ft]}],
PlotStyle -> {PointSize[0.005]}, PlotRange -> pr,
PlotLegends -> LineLegend[{"Fourier"}]],
LogPlot[Abs[disFilt[Exp[2 \[Pi]  I f /FSamp]]], {f, 0, 11000},
PlotStyle -> Green, PlotRange -> pr,
PlotLegends -> LineLegend[{"Bilinear"}]],
LogPlot[Abs[filt[I 2 \[Pi] f]], {f, 0, 11000}, PlotStyle -> Red,
PlotRange -> pr, PlotLegends -> LineLegend[{"Filter Function"}]]
]


The check with Fourier overlays the discreet filter generated with the bilinear transform. To look at this in more detail I calculate the error between the two.

err = ft - Flatten[(dis[Exp[2 \[Pi]  I # /FSamp]] & /@ ff)];
ListLinePlot[Transpose[{ff, Abs[err]}], PlotRange -> All,
Frame -> True]
`

The error is small and due to numerical calculations.

I hope that helps.