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

Learn moreOpportunities for recent engineering grads.

Apply Today
Asked by Donald Hume on 20 Jun 2011

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.

Answer by Jan Simon on 20 Jun 2011

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; 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.

Show 2 older comments

Jan Simon on 23 Jun 2011

The array z contains the current status of the filter. It allows to run FILTER piecewise and getting the same result:

x = rand(200, 1);

y1 = filter(b, a, x)

[y2a, z] = filter(b, a, x(1:100)); y2b = filter(b, a, x(101:200), z);

Now y1 and [y2a, y2b] are equal, because the final conditions of the first half are used as initial conditions of the second half.

See the help of FILTER and look in the code of FILTFILT.

I did not invent the direct form II algorithm. It is shown in "doc filter" also.

Brent on 26 Feb 2014

I'm curious, is this function definition still relevant? Upon running the code, it appears that Matlab's filter() function does not output the same results. It appears that the newer implementation takes into account a window which has the same size as the vector 'b' that is passed as an argument.

Could someone update this code above to reflect the behavior of the current Matlab filter function?

Answer by khatereh 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

## 0 Comments