Why my own direct form II filter function performs different with the Matlab function 'filter'

I'm trying to write my own sequential filter function and here is my code

function [yn,w_buff] = Sequential_filter(b,a,xn,w_buff)

% Input: 
% b     : 1*(N+1) vector that contains the numerator coefficients for a N-order
%         filter.
% a     : 1*(N+1) vector that contains the denominator coefficients for the
%         N-order filter, where a(1) is assumed to be 1.
% xn    : The input data for the nth sample, which is supposed to be a scalar
% w_buff: (N+1)*1 vector that contains the buff data. 

% Output:
% yn    : The filtered data for the nth sample.
% w_buff: (N+1)*1 vector that returns the updated buff data.

% The transfer function of system is 
% H[z] = \frac{ \sum_{k=0}^{N} b_{k} z^{-k} }{ \sum_{k=0}^{N} a_{k}z^{-k} }
% and it can be divided into two subsystems as
% H_{1}[z] = \sum_{k=0}^{N} b_{k} z^{-k}
% H_{2}[z] = \frac{ 1 }{ \sum_{k=0}^{N} a_{k} z^{-k} }
% Then the signal flow graph is
% x[n] --> H_{2} --> w[n] --> H_{1} --> y[n]

w_buff      = circshift(w_buff,1,1);

wn          = xn - a(1,2:end)*w_buff(2:end,1); 

w_buff(1,1) = wn;

yn          = b*w_buff;

end

Then I test my sequential filter function and compared it with the Matlab function filter as follows

[b,a]       = sos2tf(SOS,G);
y_out1      = filter(b,a,x_in);

w_buff      = zeros(N+1,1);
y_out2      = zeros(length(x_in),1);

for n = 1:length(x_in)

    [y_out2(n,1),w_buff] = Sequential_filter(b,a,x_in(n,1),w_buff)

end

Err         = y_out1 - y_out2;

But the performance is significant different with the Matlab function filter, especially for some filters which holds the very close together poles. Why this happened? Is there any logic bug in my program? Or Matlab has optimized filter function to reduce the cumulative error?

I sincerely look forward to your instruction and thank you very very much.

Answer 1

[Shivam Gothi]

Hello @Wentong

The logic used by you is perfectly correct and works well.

Cause of issue:

I tried to use your function on a third order Butterworth filter (refer : https://www.mathworks.com/help/signal/ref/sos2tf.html#mw_001efa27-2134-4877-99e0-b15c426487c0)

whose 3db frequency is 0.25 Hz. The [b,a] coefficients of the filter were evaluated as shown in the below code:

I obtained the following outputs.

It is seen that your filter is becoming unstable.

SOLUTION:

Please note that in the transfer function of a Butterworth filter the coefficient ( a(1) ) is not equal to 1. However, the user-defined function "Sequential_filter" that you provided requires ( a(1) = 1 ) to comply with the "direct form 2" syntax. Therefore, it is necessary to normalize the coefficients ([b, a]) so that ( a(1) = 1 ), as demonstrated below:

This resolves the issue you are facing.

I am attaching the code below:

sos = [2 4 2 6 0 2;3 3 0 6 0 0];
[b,a] = sos2tf(sos,1)

b = 1×4

6 18 18 6

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

a = 1×4

36 0 12 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x_in=sin(2*pi*0.25*(0:0.001:4));
N=3;
y_out1      = filter(b,a,x_in);
w_buff      = zeros(N+1,1);
y_out2      = zeros(1,length(x_in));
for n = 1:length(x_in)
[y_out2(n),w_buff] = Sequential_filter(b,a,x_in(n),w_buff);
end
%%%%%  PLOT THE ERROR BETWEEN TWO FUNCTIONS %%%%%%%%%%%%%%%%%%
plot(y_out2-y_out1);

function [yn,w_buff] = Sequential_filter(b,a,xn,w_buff)
    % Input: 
    % b     : 1*(N+1) vector that contains the numerator coefficients for a N-order
    %         filter.
    % a     : 1*(N+1) vector that contains the denominator coefficients for the
    %         N-order filter, where a(1) is assumed to be 1.
    % xn    : The input data for the nth sample, which is supposed to be a scalar
    % w_buff: (N+1)*1 vector that contains the buff data. 
    % Output:
    % yn    : The filtered data for the nth sample.
    % w_buff: (N+1)*1 vector that returns the updated buff data.
    % The transfer function of system is 
    % H[z] = \frac{ \sum_{k=0}^{N} b_{k} z^{-k} }{ \sum_{k=0}^{N} a_{k}z^{-k} }
    % and it can be divided into two subsystems as
    % H_{1}[z] = \sum_{k=0}^{N} b_{k} z^{-k}
    % H_{2}[z] = \frac{ 1 }{ \sum_{k=0}^{N} a_{k} z^{-k} }
    % Then the signal flow graph is
    % x[n] --> H_{2} --> w[n] --> H_{1} --> y[n]
    b=b/a(1);
    a=a/a(1);
    w_buff      = circshift(w_buff,1,1);
    wn          = xn - a(1,2:end)*w_buff(2:end,1);
    w_buff(1,1) = wn;
    yn          = b*w_buff;
end

I hope this helps !

[Wentong]

Dear @Shivam Gothi

I greatly appreciate your nice and enthusiastic response, it is very helpfull. However, I can't close this question yet because this might not be the issue.

I'm trying to reproduce a highpass filter mentioned in reference https://www.sciencedirect.com/science/article/pii/S1474667016432192

and the output of my sequential filter function is still very strange though all coefficients have been normalized by a(1).

I attached the input data and the Matlab code to this comment, and maybe this will make it more explicit for you and other people to comprehend the issue I'm having.

Once again, thank you for your wonderful assistance ! I'm very grateful.

[Shivam Gothi]

Hello @Wentong

According to my understanding, the code for the function "Sequence_filter" is correct.

I used the same file as attached by you in the above connemt, but instead of passing "x_in" input from "x_in.mat" file, I generated "x_in" by below code:

x_in = sin(2*pi*0.25*[0:0.01:12])

The output I got is shown below:

We can see that the code is working as per expectation and error is also negligible.

CASE 2: PASSING "x_in" FROM "x_in.mat" FILE.

I have plotted your "x_in" signal (attached in below image). We can see that it is unbounded and has a very large amplitude. Therefore, as per my understanding, the difference between "y_out2" and "y_out1" in this case might be due to "quantization error". I am suggesting one of the way to reduce this error.

Workaround and Possible solution:

It is possible to achieve the desired results by passing "x_in.mat" as an input to your "sequence_filter".

You can form an equivalent system consisting of two cascade blocks of second order system. This approach is demonstrated in below figure.

Now, use your "Sequence_filter" function on each of the blocks above. I have tried this approach and attaching the results obtained.

We can see that using this approach error is significantly low as compared to what it was earlier for same input "x_in.mat". I think this should validate the working of "sequential_filter" code.

I am also attaching the code for implementing the filter using above approach. You can manipulate it according to your need and also try it on different inputs.

[Wentong]

Forgive me for taking your answer a little late.

I appreciate your thoughtful response and hope the best for you!