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"} ]
NSamples = 1024
. \$\endgroup\$