Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

New to MATLAB?

Use Filter Constants to Hard Code Filter

Asked by Donald Hume

Donald Hume (view profile)

on 20 Jun 2011
Latest activity Commented on by TJ Sartain

TJ Sartain (view profile)

on 11 Mar 2015 at 21:41

Hey,

      I am trying to implement a real-time filter so am using MATLAB's butter() function to generate the needed [b,a] vectors
    [b,a] = butter(4,4/35,'low');
      Just to be clear I have used these generated vectors with the filter(a,b,data) function to successfully filter my data and it looks quite desirable. But as I am in the end trying to create a real time filter, I am trying to code this into a for-loop (for testing purposes). My code is as follows:
    for n=5:1:length(x)
        y(n) = b(1)*x(n)+b(2)*x(n-1)+b(3)*x(n-2)+b(4)*x(n-3)+b(5)*x(n-4)-a(2)*y(n-1)-a(3)*y(n-2)+a(4)*y(n-3)+a(5)*y(n-4);
    end

This is the mathematical representation as far as I can gather from the doc: http://www.mathworks.com/help/techdoc/ref/filter.html

Can anyone tell me how I am incorrectly modeling the filter() command? I have also switched the a, b, column vectors (in case that was an issue). The above method just goes to infinity, and with a<->b the data just seems to be amplified.

Thanks for the help in advance.

0 Comments

Donald Hume

Donald Hume (view profile)

2 Answers

Answer by Jan Simon

Jan Simon (view profile)

on 20 Jun 2011
Edited by Jan Simon

Jan Simon (view profile)

on 26 Oct 2014
Accepted answer

The difference equation looks ok, but you do not show how e.g. "y(n-4)" is initialized.

Matlab's FILTER uses the "direct form II transposed" implementation, which is more efficient. Together with inital and final conditions:

function [Y, z] = myFilter(b, a, X, z)
% Author: Jan Simon, Heidelberg, (C) 2011
n    = length(a);
z(n) = 0;      % Creates zeros if input z is omitted
b = b / a(1);  % [Edited, Jan, 26-Oct-2014, normalize parameters]
a = a / a(1);
Y    = zeros(size(X));
for m = 1:length(Y)
   Y(m) = b(1) * X(m) + z(1);
   for i = 2:n
      z(i - 1) = b(i) * X(m) + z(i) - a(i) * Y(m);
   end
end
z = z(1:n - 1);

 

[EDITED]: A C-Mex implementation which handles arrays also: FEX: FilterM.

12 Comments

Harry

Harry (view profile)

on 26 Oct 2014

The difference of 0.00001 does not mean that the algorithm I've posted differs from Matlab's substantially. It means it differs by 0.00001 in this example. I would also prefer the direct form II transposed implementation for numerical stability.

Can we put the normalisation inside the myFilter() function? I tried inserting these two lines at the top:

a = a/a(1);
b = b/a(1);

but still get very large errors.

Jan Simon

Jan Simon (view profile)

on 26 Oct 2014

Seriously?

After you have normalized a by a(1), a(1) is 1.0 . Afterwards normalizing b by a(1) does not change the value of b. Therefore you have to normalize b at first.

b = b / a(1);
a = a / a(1);

And immediately I get 0.0 difference between Matlab's filter() and the above posted M-version.

When the output of the built-in function and the M-version does not differ by rounding errors for random input data, it is extremly likely that the same algorithm is applied. If the algorithm you've posted replies a difference by 1e-5 it does not seem to be the same algorithm, obviously. And even Matlab's documentation mentions the direct form II approach, which is more stable and faster, but it is not explicitly explained, that this algorithm is implemented. But I do not have any doubts about this.

TJ Sartain

TJ Sartain (view profile)

on 11 Mar 2015 at 21:41

I'm getting slight discrepancies on a 6000+ set of numbers in the range of -.0003 to .0003. Very accurate but not exactly the same result as Matlab. I'm implementing the algorithm as part of a custom situation using Obj-C. There may be other things causing the problem but can anyone verify this pseudocode of Jan's code here:

 Set n to be the number of elements in a.
 Make sure z is of length n. Add on as many 0 elements as needed.
 Normalize b by dividing every elemnt by the first element in a.
 Next, normalize a by also dividing by the first element in a.
 Create array Y to be as long as X and filled with 0.
 Loop with m indexing the entire length of Y
    Set Y(m) to be the first element of b times X(m) plus the first element of z.
    Loop with i indexing the length of a but starting at the 2nd position.
        Set z(i-1) = b(i) * X(m) + z(i) - a(i) * Y(m)
 Truncate z down to size n-1
Jan Simon

Jan Simon (view profile)

Answer by khatereh

khatereh (view profile)

on 6 Jan 2012

Hi, I want to use your function instead of matlab filter function. I calculated the filter coefficient in b matrix and it is FIR filter so all the a values are 1. What should be the value for z in my case? I am confused how should I use z.

Thanks so much. Regards, KHatereh

1 Comment

Jan Simon

Jan Simon (view profile)

on 26 Oct 2014

The meaning of z is explained in the documentation: doc filter.

The initial conditions for the internal state of the filter can be set such, that the transient effects are damped. Look into the code of the filtfilt function for a method to do this automatically.

Set z to zero, if you do not have any information about the signal.

For me the meaning of z got clear, when I examined this example: Imagine a long signal X, which is divided in 2 parts X1 and X2. Now the complete signal X is filtered with certain parameters and the initial settings z=0 (this means zeros(1,n-1) with n is the length of the filter parameters):

z = zeros(1, length(b) - 1);
Y = filter(b, a, X, z);

Now we do this for the first part:

z = zeros(1, length(b) - 1);
[Y1, z1] = filter(b, a, X1, z);

Now the output z1 is the internal state of the filter, a kind of history over the last elements. If we use the output z1 of the 1st part as input of the 2nd, we get exactly the same outpt as for the full signal:

Y2 = filter(b, a, X2, z1);
isequal(Y, [Y1, Y2])   % TRUE

But if we omit z1 as input for filtering X2, there is a small difference mostly at the start of Y2 due to the transient effects.

In this case, we do have some knowledge about the history of the internal filter state for X2, but for X1 this state is not defined and zeros are a fair guess, but not necessarily smart.

khatereh

khatereh (view profile)

Contact us