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)
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 !
