2
$\begingroup$

I want to design a simple IIR filter with mathematica. For example, consider the filter deigned with this online IIR calculator:

y[n] = x[n- 2] + 2*x[n- 1] + x[n- 0] - 0.8752145483*y[n- 2] + 1.8668922797*y[n- 1]

Butterworth, Lowpass, filter order: 2, sample rate: 1000, corner frequency 1: 15

I can certainly envision just using entering this equation into Mathematica with a For[] loop, but that doesn't sound satisfactory to me. Especially for a higher order IIR filter.

Needless to say the mathematica page about this is way above my pay grade. I just want to implement a simple filer. It is easy in Matlab, so it should be easier in Mathematica, right?

My question is:

How can I implement a simple low pass IIR filter with Mathematica?

The correct answer will allow all the order to increase to a higher number than 2. It will can use the in-built Mathematica functions, or some other custom means. The correct answer will not have a procedural loop. I assume the equation above can somehow be entered using a Mathematica style list?

$\endgroup$
  • $\begingroup$ If you want to "increase the order to higher much higher than 2", remove the word "simple" from your question. $\endgroup$ – andre314 Jul 30 '19 at 6:34
  • 1
    $\begingroup$ @andre314 It would not be so hard to do and order > 2 with a procedural loop. The online IIR calculator linked above even gives you the code (in 'c') for a loop up to 10th order. Check it out.I could just type it in. $\endgroup$ – axsvl77 Jul 30 '19 at 6:38
  • $\begingroup$ This is straightforward in mathematica. Above you say corner frequency 1:15. Do you mean 1.15 Hz? I will give you some code shortly. $\endgroup$ – Hugh Jul 30 '19 at 13:28
5
$\begingroup$

Here is some simple code to do what you want

  sr = 1000.;(* Sample rate Hz*)
  cf = 1.15; (* Corner frequency Hz *)
  order = 2; (* Filter order *)
  filt = ButterworthFilterModel[{"Lowpass", order ,2 π cf}]; (* make filter *)
  rfilt = ToDiscreteTimeModel[filt,1/sr]; (* Convert to recurrence filter *)

The ButterworthFilterModel makes the filter in the s-plain. The ToDiscreteTimeModel makes the recurrence filter. You don't have to go to a site that gives you filter coefficients.

Now let's make some data and then filter it. The data is a sine wave with lots of noise;

nn = 5000; (* number of points in data set *)
data = Table[Sin[2 π 0.5 t], {t, 0, (nn - 1)/sr , 1/sr}] + 
   RandomReal[{-1, 1}, nn];
fdata = RecurrenceFilter[rfilt, data];

The RecurrenceFilter does what you were suggesting one could do in a For loop.

Here is the data before and after filtering.

ListLinePlot[data]
ListLinePlot[fdata]

Mathematica graphics

Mathematica graphics

There are some extra minor points, we should pre-warp, but that may not be necessary for you. I may add these later.

Hope that helps.

$\endgroup$
  • $\begingroup$ Is changing the order from 2 to 200 is as trivial as it looks? $\endgroup$ – axsvl77 Jul 30 '19 at 14:20
  • $\begingroup$ Yes but the filter will probably be unstable due to numerical accuracy. Do you really want 200? This would be 200 poles in the filter! What exactly are you trying to do? The order is the order of the filter. If you really want an order of 200 you could make a filter of order of 2 and then filter it 100 times. $\endgroup$ – Hugh Jul 30 '19 at 14:23
  • $\begingroup$ Thanks for the comment. It sounds like I need to go back and read the introductory text again! $\endgroup$ – axsvl77 Jul 30 '19 at 14:26
  • $\begingroup$ @Hugh If you cascade n Butterworth filters, the whole is not a Butterwoth filter (this is allready true for =2). $\endgroup$ – andre314 Jul 30 '19 at 19:21
  • $\begingroup$ @andre314 Where are you getting that from? It is normal to cascade filters rather than going to higher orders. The standard unit is the second order filter. What do you take as a Butterworth filter? $\endgroup$ – Hugh Jul 30 '19 at 19:25
1
$\begingroup$

Testing/designing filters can be pretty 'trivial' if you're not too worried about understanding what you're doing, and can quickly tested with different filters.

We can find filters via:

ToExpression /@ Names["*FilterModel"]
(*{BesselFilterModel, BiquadraticFilterModel, ButterworthFilterModel, 
Chebyshev1FilterModel, Chebyshev2FilterModel, EllipticFilterModel}*)

You wanted a BesselFilterModel[]

We now generate a signal, pass our filter params we want to test, and chop them up a bit:

sp = 0.001
u = With[{\[Omega] = 3/2}, Sin[\[Omega] t] + 1/2 Sin[5 \[Omega] t] + 1/4 Sin[8 \[Omega] t]];
signal = Table[u, {t, 0, 12, sp}];


fillers = Table[BesselFilterModel[{i, j}], {i, 1, 2}, {j, 1, 2}] // Flatten;
p1 = Chop[TransferFunctionExpand /@ N[fillers]];
discrete = Chop[ToDiscreteTimeModel[#, sp]] & /@ p1;
things = Table[Flatten[OutputResponse[discrete[[i]], u, {t, 0, 12}]], {i, 1, Length[discrete]}];
Show[ListLinePlot[signal], ListLinePlot[things]]

thing

And to add an extra bit of fun, since you're looking for a iir, I assume you may want to eventually code it on an embedded platform...i do this quite a bit lately, using MicrocontrollerKit we can easily generate some code ( I use Wiring, but C code can also be generated ).

Needs["MicrocontrollerKit`"]

{\[ScriptCapitalM]1, \[ScriptCapitalM]2, \[ScriptCapitalM]3, \[ScriptCapitalM]4} = Table[MicrocontrollerEmbedCode[sys, <|"Target" -> "ArduinoUno", "Inputs" -> {"A0" -> "Analog"}, "Outputs" -> {"Serial"}|>, <|"ConnectionPort" -> None|>, <|"Language" -> "Wiring"|>], {sys, {discrete[[1]], discrete[[2]], discrete[[3]], discrete[[4]]}}]

Text[Grid[{Table[M["SourceCode"], {M, {\[ScriptCapitalM]1, \[ScriptCapitalM]2, \[ScriptCapitalM]3, \[ScriptCapitalM]4}}]}, Frame -> True]]

codes

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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