15
$\begingroup$

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');
|improve this question
$\endgroup$
  • 4
    $\begingroup$ Perhaps, for the benefit of those readers who do not know Matlab, you could explain what the more general expression means? $\endgroup$ – whuber Mar 14 '13 at 5:12
  • $\begingroup$ Could you please add an example to your question that makes clear what the a parameter of that function does? $\endgroup$ – Mr.Wizard♦ Mar 14 '13 at 5:33
26
$\begingroup$

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.

|improve this answer
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Mr.Wizard♦ Mar 14 '13 at 6:18
  • $\begingroup$ 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. $\endgroup$ – rm -rf♦ Mar 14 '13 at 6:21
  • $\begingroup$ 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) $\endgroup$ – Mr.Wizard♦ Mar 14 '13 at 6:38
  • $\begingroup$ 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). $\endgroup$ – rm -rf♦ Mar 14 '13 at 6:43
16
$\begingroup$

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

and return the same answers.

|improve this answer
$\endgroup$
  • $\begingroup$ Nice! I should look into more of the new functions in v9 $\endgroup$ – rm -rf♦ Mar 15 '13 at 18:24
  • 1
    $\begingroup$ 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. $\endgroup$ – bill s Mar 16 '13 at 7:05
  • $\begingroup$ @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? $\endgroup$ – chyanog Mar 24 '13 at 6:00
  • 1
    $\begingroup$ 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. $\endgroup$ – bill s Mar 24 '13 at 14:33
8
$\begingroup$

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.

|improve this answer
$\endgroup$
  • 1
    $\begingroup$ I suspect the intent expressed in the Matlab reference might be a ratio of convolutions, but it's far from clear. $\endgroup$ – whuber Mar 14 '13 at 5:43

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.