# What is Mathematica's equivalent to MATLAB's filter function?

The MATLAB code

filter(0.5,[1, -0.5], [1:10])


is equivalent to

Rest@FoldList[(#1 + #2)/2. &, 0, Range[10]]


I don't know how to implement something more general like，filter([1,2,3],[4,5,6], [1:10]) in Mathematica.

I'm trying to rewrite a snippet of MATLAB code to Mathematica. I'm just interested in the filter function and there is no other purpose. What is its equivalent or how can I implement it?

v = [0.0 + 2 j; -sqrt (3) - 1 j; sqrt (3) - 1 j];
n = randi (3, 1, 10000000);
p = filter (0.5, [1 - 0.5], v (n));
plot (p, '.b');

-
Perhaps, for the benefit of those readers who do not know Matlab, you could explain what the more general expression means? – whuber Mar 14 '13 at 5:12
Could you please add an example to your question that makes clear what the a parameter of that function does? – Mr.Wizard Mar 14 '13 at 5:33

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.

-
I'm still trying to understand this, and your code. Does a have to be shorter than b by a certain factor? I get RecursionLimit errors with some of the input I tried. – Mr.Wizard Mar 14 '13 at 6:18
Ah, good catch. Add y[0 | 0. | _?Negative] := 0; (edited). And to answer your question, no. I'm not sure it can be optimized any more without risking losing clarity. The reason for the X definition is so that I can easily set X[_] = 0, whereas if you stuck with manipulating the list x, you'll end up doing other kinds of code gymnastics when it gets to x[[-1]], etc. – R. M. Mar 14 '13 at 6:21
Yes, it's not as hairy as it looks all balled up in the post. Nevertheless I think I'll try doing this with FoldList or something a bit later, but I do imagine clarity will be lost, even if I succeed. At the moment I have something else to work on. (and you already have my +1) – Mr.Wizard Mar 14 '13 at 6:38
Clarity will indeed be lost even if you succeed with FoldList (and I don't have high hopes of that :P). The problem with using FoldList is that you can have arbitrary length filters (with arbitrary a which will decide the number of previous steps to keep). – R. M. Mar 14 '13 at 6:43

In Mathematica 9, there is a more direct way using the function RecurrenceFilter.

For example, the two examples from MATLAB can be done straightforwardly as:

 RecurrenceFilter[{{1}, {1, 2}}, Range[6]]


and

RecurrenceFilter[{{3, 2}, {1, 3, 1}}, Range[6]] // N


-
Nice! I should look into more of the new functions in v9 – R. M. Mar 15 '13 at 18:24
There are a number of new functions (and improved functionality) in 9, especially in the signal processing area. Functions like RecurrenceFilter will make it easier for people to transition to Mathematica from Matlab. – bill s Mar 16 '13 at 7:05
@bills, I have another question, for two dimension list, RecurrenceFilter[{{1, -0.5}, {0.5}}, {{1, 2}, {3, 4}, {5, 6}}] isn't equivalent to matlab code filter([0.5], [1,-0.5], [[1,2]; [3,4]; [5,6]]), RecurrenceFilter[{{1, -0.5}, {0.5}}, {x, y, z}] /. {x -> {1, 2}, y -> {3, 4}, z -> {5, 6}} is equivalent to it, but I think it's not efficient , is there a more direct way? – chyaong Mar 24 '13 at 6:00
chyanog -- I think you can just use RecurrenceFilter[{{1, -0.5}, {0.5}}, Transpose[{{1, 2}, {3, 4}, {5, 6}}]] . Neither of these (neither Matlab nor Mathematica) are doing 2D recursion, they are both just applying the filter to the rows (Mathematica) or the columns (Matlab) separately. – bill s Mar 24 '13 at 14:33

I think you want ListConvolve or ListCorrelate.
You can implement the example on the linked page like this:

data = Range[1, 4, 0.2];
windowSize = 5;
ones = ConstantArray;
filter = ListCorrelate[#, #3, -1, 0] &;
filter[ones[1, windowSize]/windowSize, 1, data]

{0.2, 0.44, 0.72, 1.04, 1.4, 1.6, 1.8, 2., 2.2, 2.4, 2.6, 2.8, 3., 3.2, 3.4, 3.6}


Please note this is not complete: I don't know how the "denominator coefficient vector a" is supposed to be used and the example on that page doesn't make it clear to me yet. Also, I guessed as what ones does and seem to have guessed right, but I didn't look it up.

-
I suspect the intent expressed in the Matlab reference might be a ratio of convolutions, but it's far from clear. – whuber Mar 14 '13 at 5:43