1
$\begingroup$

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:

enter image description here

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:

enter image description here

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:

enter image description here

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"} ]
$\endgroup$

1 Answer 1

2
$\begingroup$

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]

enter image description here

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]

enter image description here

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

enter image description here

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]

enter image description here

The error is small and due to numerical calculations.

I hope that helps.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.