Can FILTFILT(B, A, X) yeild the same results as FILTFILT(SOS, G, X) assuming inputs are generated via [z,p,k] = butter(3,Wn) + [sos,g] = zp2sos(z,p,k) and [A,B] = butter(3,Wn)

6 views (last 30 days)
Chris Heyman
Chris Heyman on 18 Apr 2022
Answered: Paul on 18 Apr 2022
Hi all,
I am working on a signal processing and bandpass filtering problem which uses a butterworth filter. I was recently advised to try using the [zero, pole, gain] output instead of the [A,B] output (this output returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator)) for increased reliability between matlab versions.
I do not understand the nuances of these different methods. I am getting different result (tradfilt, h_signal, thephase) with each method and want to understand why.
Method 1 Code
Wn = [Fpass(1,1)/(Fs/2) Fpass(1,2)/(Fs/2)]; % normalized by the nyquist frequency
[B,A]= butter(3,Wn)
tradfilt = filtfilt(B,A,signal);
h_signal=hilbert(tradfilt);
thephase=atan2(tradfilt, imag(h_signal));
Method 2 code
Wn = [Fpass(1,1)/(Fs/2) Fpass(1,2)/(Fs/2)]; % normalized by the nyquist frequency
[z,p,k] = butter(3,Wn);
[sos,g] = zp2sos(z,p,k);
tradfilt= filtfilt(sos,g,signal);
h_signal=hilbert(tradfilt);
thephase=atan2(tradfilt, imag(h_signal));
Thanks!

Accepted Answer

Paul
Paul on 18 Apr 2022
Hello Chris,
Here's the code I ran (unknown to me, butter() assumes a bandpass by default with a two-element frequency vector, which is not documented as far as I can tell).
Fs = 1;
Fpass = [0.1 .3];
rng(100);
signal = randn(1,100);
Wn = [Fpass(1,1)/(Fs/2) Fpass(1,2)/(Fs/2)]; % normalized by the nyquist frequency
[B,A]= butter(3,Wn);
tradfilt1 = filtfilt(B,A,signal);
%h_signal=hilbert(tradfilt);
%thephase=atan2(tradfilt, imag(h_signal));
%Wn = [Fpass(1,1)/(Fs/2) Fpass(1,2)/(Fs/2)]; % normalized by the nyquist frequency
[z,p,k] = butter(3,Wn);
[sos,g] = zp2sos(z,p,k);
tradfilt2= filtfilt(sos,g,signal);
figure;
plot(0:99,tradfilt1,'-o',0:99,tradfilt2,'-x')
figure
plot(0:99,tradfilt1 - tradfilt2)
Sure enough, the filtfilt() outputs are different, particularly at the edges. I think the reason for this is that with b/a inputs, filtfilt() makes a single forward pass over the input and then a single reverse pass. For each pass, filtfilt() prepends the input signal with initial conditions based on the input signal to minimize transients around the start/end of the result. But for the sos input, with three sections, filtifilt() does the forward/reverse pass three times (once for each section) and comes up with new initial conditions for each pass, which depend on the results of the previous pass, or the input signal in the case of the first forward/reverse pass. So it stands to reason that the response at the edges won't match between the b/a and sos inputs. I don't have an opinion on whether or not that's how filtfilt() should work, but I'm quire sure that's how filtfilt() does work, at least for this example.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!