Use Filter Constants to Hard Code Filter

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 1

[Jan]

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.

[Donald Hume]

I just initialize y(1) to y(4) to 0 before the code. sorry for not including! Your code looks very similar to mine, I'll take a closer look tomorrow at work.

[Donald Hume]

Wow this is a great answer. Could you shed some insight into why my code wasn't working properly? I set initial conditions for y(1)-y(4) to zero (not shown). Thanks so much for your help. I will be sure to leave your name in the final implementation.

[Donald Hume]

Also, I'm not sure what the z array is.

[Jan]

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]

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?

[Harry]

I also couldn't get this to produce the same results as the built-in function. In fact, after a lot of searching online, I couldn't find anything.

However, in GNU Octave's help file for the filter() function, it gives a simple equation that is easy to implement (see comments at top of my code, below).

Here is some code, optimised for readability (not speed or generality) that seems to give me the same results as Matlab's filter() function:

function y = harry_filter(b, a, x)
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
% function y = harry_filter(b, a, x)
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
% From GNU Octave's filter() function help:
%
% An equivalent form of the equation is:
%
%            N                   M
%  y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k)  for 1<=n<=length(x)
%           k=1                 k=0
%
% where N=length(a)-1, M=length(b)-1, c = a/a(1) and d = b/a(1).
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

% Useful variables
N = length(a)-1;
M = length(b)-1;
c = a(:)/a(1);
d = b(:)/a(1);

% Pre-allocate y
y = zeros(1, length(x));

for n = 1:length(x)
    % First sum term
    y1 = 0;
    for k = 1:N
        if n-k >= 1
            y1 = y1 + c(k+1) * y(n-k);
        end
    end
    % Second sum term
    y2 = 0;
    for k = 0:M
        if n-k >= 1
            y2 = y2 + d(k+1) * x(n-k);
        end
    end    
    y(n) = -y1 + y2;
end

[Jan]

When I read the documentation of filter, I find under "More about": Rational transfer function the pseudo code for the direct form II transposed implementation. At least until Matlab R2011b the code above created exactly the same output as Matlab's filter(). If you observe differences, please post some example data and explain, if this concern rounding errors.

What do you observe for this code:

B = [0.000416546, 0.00124964, 0.00124964, 0.000416546];
A = [1, -2.68616, 2.41966, -0.730165];      % Lowpass Butterworth
t = linspace(0, 2*pi, 400);
x = sin(t) + rand(1, numel(t)) * 0.1;
y1 = filter(B,A,x);
y2 = myFilter(B,A,x);

max(abs(y1-y2))

Octave's FILTER cannot handle initial conditions as Matlab's or my M-coded version.

[Harry]

Octave's filter() can handle initial conditions just fine. However, like Matlab, Octave provides a fast pre-compiled implementation (with no Matlab source code), so I just gave some stripped-down source code because that's what I just spent several hours looking for. Maybe it'll be useful to someone else.

The specific example you gave produced zero error. However, if we instead use a few pseudorandom numbers:

randn('state',0);
B = randn(1,4);
A = randn(1,4);
x = randn(1,50);
y1 = filter(B,A,x);
y2 = myFilter(B,A,x);
max(abs(y1-y2))

we get an error of 2.7526e+10. (The code I posted gives 0.00001 on the same data). It is possibly also worth mentioning that your code can't handle length(a) > length(b).

[Jan]

Thanks for your reply. Yes, of course this must produce a big difference. :-) I've posted the absolute core of the filter function only and omitted all input testing and other work, which ist required for a convenient and productive code. Here you observe, that my code expect A to be normalized, although Matlab's filter cares about this internally, see doc filter:

 If a(1) is not equal to 1, then filter normalizes the filter coefficients by a(1). Therefore, a(1) must be nonzero.

Please try this:

randn('state',0);
B = randn(1,4);
A = [1.0, randn(1,3)];   % Normalized! A(1)==1.0
x = randn(1,50);
y1 = filter(B,A,x);
y2 = myFilter(B,A,x);
max(abs(y1-y2))          % I get 0 difference

The code you've posted normalizes A also. The difference of 0.00001 means, that the algorithm you've posted differs from Matlab's substantially. The direct form II transposed approach has a larger numerical stability and therefore I'd prefer it, and for me it still seems like Matlab uses it also.

[Harry]

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]

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]

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]

You cannot expect that the code creates exactly the same results when you run it in a different language. Notice that Matlab sets the precision of the floating point engine to 53 bits, while this might be set to 64 bits with your compiler. The order of the computations matters also, but a compiler can reorder the terms for an internal optimization. Therefore it is expected, that the shown algorithm creates slightly different results in the magnitude of EPS(input values), and of course these rounding errors can accumulate for long signals. Remember, that in this case Matlab's filter is not the gold standard. Only a computation with e.g. quadrupel precision would clarify, which solution is "better".

[Han Zerui]

Could you please explain why you set z(n) = 0 at first and only return z(1:n-1) at last?

[Jan]

@Han Zerui: z has one element less than a, because it contains the "former" states. Therefore z(n) is always 0. Therefore the inner loop in my code:

   for i = 2:n
      z(i - 1) = b(i) * X(m) + z(i) - a(i) * Y(m);
   end   

could be written as:

   for i = 2:n - 1
      z(i - 1) = b(i) * X(m) + z(i) - a(i) * Y(m);
   end
   z(n - 1) = b(n) * X(m) - a(n) * Y(m);

[Han Zerui]

@Jan Thank you for the answer. That makes so much more sense.

[Francely Guzmán Otagri]

Hello, I have two questions.

First, I am assuming that this equation

is referenced by this line of code:

Y(m) = b(1) * X(m) + z(1);

I am right? If yes, I don't understand why you only evaluate the first component of 'b'. (b(1))

And the second question is regarding 'z'. Why do you take as initial condition for 'Y(m)' the first component of the resulting vector 'z' that appears at the end of the code?

Thanks.I hope you can solve my question as soon as possible.

[Martin Goubej]

@Francely Guzmán Otagri The provided code (as well as the Matlab filter command) implements the Transposed-Direct-Form-II structure of the IIR filter, check e.g. this page for explanation https://ccrma.stanford.edu/~jos/filters/Transposed_Direct_Forms.html#17507. The first line implements the upper summation term from the Fig. 9.2 and the rest deals with the update of the delay blocks.

As for the 'z' question - this comes form the filter structure, the uppermost summation term depends directly only on the output of the first internal delay. This is why it is assigned to the actual output value Y(m).

[Jan]

With some more features:

function [Y, z] = myFilter(b, a, X, z)
% Author: Jan, Heidelberg, (C) 2022
na = numel(a);  % Equilize size of a and b
n  = numel(b);
if na > n
    b(na) = 0;
    n     = na;
elseif na < n
    a(n)  = 0;
end

z(n) = 0;       % Creates zeros if input z is omitted
b = b / a(1);   % Normalize a and b
a = a / a(1);

Y = zeros(size(X));
if na > 1       % IIR filter
   for m = 1:length(Y)
      Y(m) = b(1) * X(m) + z(1);
      for i = 2:n-1
         z(i - 1) = b(i) * X(m) + z(i) - a(i) * Y(m);
      end
      z(n - 1) = b(n) * X(m) - a(n) * Y(m);  % Omit z(n), which is 0
   end

else            % FIR filter: a(2:n) = 0
   for m = 1:length(Y)
      Y(m) = b(1) * X(m) + z(1);
      for i = 2:n-1
         z(i - 1) = b(i) * X(m) + z(i);
      end
      z(n - 1) = b(n) * X(m);  % Omit z(n), which is 0
   end    
end
z = z(1:n - 1);
end

Answer 2

[khatereh]

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

[Jan]

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.

Answer 3

[Yves]

Can someone please comment on whether this z/z1 or zi/zf - initial/final condition (delay) of the digital filter is equivalent to the state variables in state-space model (ABCD matrix) of the filter?

Answer 4

[Michal]

Probably fully functional naive version of filter:

function [Y, z] = myFilter(b, a, X, z)
    % a and b should have same order
    na = length(a);
    nb = length(b);
    if na > nb
        b = [b,zeros(1,na-nb)];
        n = na;
    elseif na < nb
        a = [a,zeros(1,nb-na)];
        n = nb;
    else
        n = na;
    end
    % naive filter implementation
    z(n-1) = 0;      % Creates zeros if input z is omitted

    b = b / a(1);    % normalize parameters
    a = a / a(1);
    
    Y    = zeros(size(X));
    for m = 1:length(Y)
        Xm = X(m);
        Y(m) = b(1) * Xm + z(1);
        Ym = Y(m);
        for i = 2:n-1
            z(i - 1) = b(i) * Xm + z(i) - a(i) * Ym;
        end
        z(n - 1) = b(n) * Xm - a(n) * Ym;
    end
end

[Jan]

A more compact method to equilize the length of a and b:

n = max(numel(a)m numel(b));
if numel(a) < n
    a(n) = 0;
elseif numel(b) < n
    b(n) = 0;
end

Bug: z(n-1) = 0; This does not only create z if it is not provided, but overwrites the n-1.th element by 0. Better: z(n)=0, because the n.th element 0 in this method.

[Michal]

Your equilization length (a and b) method is good! But, the proposed modification:

z(n) = 0

broke the compatibility with original Matlab built-in function filter in a case of use z as an initial condition.

Error using filter

Initial conditions must be a vector of length max(length(a),length(b))-1,

or an array with the leading dimension of size max(length(a),length(b))-1

and with remaining dimensions matching those of x.

[Jan]

@Michal: No, z(n)=0 does not break the compatibility with Matlab's filter, if it is applied inside the function and cropped finally: z = z(1:n - 1) as in my implementation above. The error message appears, if you expand z before calling filter, but this was not suggested.

Your method to expand z by the line:

z(n-1) = 0;

is a bug, because it sets the last element of z to 0. If z was [1,2,3] initially (then a and b have 4 elements), your code sets z to [1,2,0]. Simply try it:

x = rand(1, 100);
B = [0.000416546, 0.00124964, 0.00124964, 0.000416546];
A = [1, -2.68616, 2.41966, -0.730165];      % Lowpass Butterworth
z = [1,2,3];

y1 = filter(B, A, x, z);
y2 = michaelsFilter(B, A, x, z);
plot(y1, 'r')
hold('on')
plot(y2, 'b')  % Obviously different

function [Y, z] = michaelsFilter(b, a, X, z)
    % a and b should have same order
    na = length(a);
    nb = length(b);
    if na > nb
        b = [b,zeros(1,na-nb)];
        n = na;
    elseif na < nb
        a = [a,zeros(1,nb-na)];
        n = nb;
    else
        n = na;
    end
    % naive filter implementation
    z(n - 1) = 0;    % BUG! Must be z(n)=0, then crop z(n) after the loop. 
    % Or:
    % if length(z) < n-1
    %    z(n - 1) = 0;
    % end

    b = b / a(1);    % normalize parameters
    a = a / a(1);
    
    Y = zeros(size(X));
    for m = 1:length(Y)
        Xm = X(m);
        Y(m) = b(1) * Xm + z(1);
        Ym = Y(m);
        for i = 2:n-1
            z(i - 1) = b(i) * Xm + z(i) - a(i) * Ym;
        end
        z(n - 1) = b(n) * Xm - a(n) * Ym;
    end
end

Do you see it now?