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

Question

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?

If you want to "increase the order to higher much higher than 2", remove the word "simple" from your question. — andre314, Jul 30 '19 at 06:34
@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. — Astor Florida, Jul 30 '19 at 06:38
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. — Hugh, Jul 30 '19 at 13:28

Hugh · Answer 1 · 2019-07-30T14:00:53.620

6

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

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.

edited Jul 30 '19 at 14:00

answered Jul 30 '19 at 13:55

Hugh

16,387
3
31
83

Is changing the order from 2 to 200 is as trivial as it looks? – Astor Florida Jul 30 '19 at 14:20
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. – Hugh Jul 30 '19 at 14:23
Thanks for the comment. It sounds like I need to go back and read the introductory text again! – Astor Florida Jul 30 '19 at 14:26
1

@Hugh If you cascade n Butterworth filters, the whole is not a Butterwoth filter (this is allready true for =2). – andre314 Jul 30 '19 at 19:21
@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? – Hugh Jul 30 '19 at 19:25
@Hugh a Butterworth filter is a filter wich module of the transfert function is 1/(1 + w^n) : There are no terms below n. If you cascade 2 filters you get 1/(1 + 2 w^n + w^(2n)). The 2 w^n shouldn't be there for a Butterworth filter. Here I am in the analog world (ie not discrete) but that doesn't change the problem. – andre314 Jul 30 '19 at 19:32
@andre314 I have just had a quick look at the wikipedia entry for Butterworth filter. The original filter was certainly of the form you state. So if you do cascade the filters it is not a Butterworth filter but just a IIR filter with good properties. Does the cascade maintain the maximal flatness property? – Hugh Jul 30 '19 at 19:45
No, the cascade doesn't maintain the maximal property. For example 1/(1+w^4) is more flat that 1/(1+ 2 w^2 + w^4). Just plot them and see (do not forget to normalise the 3dB rolloff frequency). This is also true for other filters (Tsebyscheff...), though there are some subtleties in the vocabulary. More if you are interested... – andre314 Jul 30 '19 at 19:57
@andre314 Thanks I have learnt something. I suppose if you did want to implement a 200 order filter you would have to factor the polynomial into smaller orders and then cascade each smaller order. This might then give you the numerical stability. – Hugh Jul 30 '19 at 19:59
Yes it would be good to move to a dedicated chatroom, but How ? – andre314 Jul 30 '19 at 20:01
Let us continue this discussion in chat. – Hugh Jul 30 '19 at 20:01

DrMrstheMonarch · Answer 2 · 2019-07-30T20:48:22.863

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

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

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

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

2 Answers2