There is a misunderstanding of what filter
really does in the MATLAB community, largely because of its widespread use as a cheap moving average/smoother (because the actual moving average function is in a paid toolbox).
The function filter(b, a, x)
convolves the input list x
with a digital filter whose transfer function (TF) is described by the lists b
and a
. If a
is 1, then the filter is an FIR filter and can be easily implemented using ListConvolve
. If a
is a list, then the filter is an IIR filter whose response is a little more involved.
In either case, the output is given by the following difference equation (I'm using the notation in the IIR wiki page I linked to, for reference):
$$y[n] = \frac{1}{a_{0}} \left(\sum_{i=0}^P b_{i}x[n-i] - \sum_{j=1}^Q a_{j} y[n-j]\right)$$
This can be implemented in Mathematica as:
Clear@filter
filter[b_List, a_List, x_List] :=
Module[{y, L = Length@x, P = Length@b - 1, Q = Length@a - 1, X},
MapIndexed[(X[#2[[1]]] = #) &, x];
X[_] = 0;
y[0 | 0. | _?Negative] = 0;
y[n_] := y[n] = (Total[b Table[X[n - i], {i, 0, P}]] -
Total[Rest@a Table[y[n - j], {j, Q}]])/First@a;
Table[y[n], {n, 1, L}]
]
Normally, this could be solved with RecurrenceTable
(and indeed, it works for certain cases), but it doesn't sit well with arbitrary b
and a
. You can verify the results against MATLAB's filter
:
MATLAB:
filter([1,2],1,1:6)
% 1 4 7 10 13 16
filter([1,3,1],[3,2],1:6)
% 0.3333 1.4444 2.3704 3.4198 4.3868 5.4088
Mathematica:
filter[{1, 2}, {1}, Range@6]
(* {1, 4, 7, 10, 13, 16} *)
filter[{1, 3, 1}, {3, 2}, Range@6] // N
(* {0.333333, 1.44444, 2.37037, 3.41975, 4.38683, 5.40878} *)
Note that I don't do any error checking against the length of b
and a
, etc. That can be easily added, if so desired.
a
parameter of that function does? $\endgroup$