MATLAB: Use Filter Constants to Hard Code Filter

butterworthfilter designMATLABreal timeSignal Processing Toolbox

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.

Best 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.
Related Question