Why is the second-order-section cascade approach giving "worse" results?

Question

It seems that as a general rule, for digital IIR filters of order more than 2, one is encouraged to use a cascade of second order sections, for reasons connected to the poor stability of high-order filters. See, for example, my question What is the largest "safe" order for the digital Butterworth filter of a given signal? and the answers given there, or the "Limitations" section of https://uk.mathworks.com/help/signal/ref/butter.html.

Now suppose I have a digital signal $x(t)=\sin(2\pi t)$, $t \in \{0,0.08,0.16,\ldots,99.92,100\}$.

Suppose I want to filter this signal in the passband $[f_\mathrm{min},f_\mathrm{max}]=[0.5,3]$ using a zero-phase Butterworth bandpass filter, with filter order $60$; I believe the transfer function of the filter should be $$ H(2\pi if) \ = \ \frac{f^{30}}{f^{30}+\left( 0.4f^2 - 0.6 \right)^{30}} $$ where 0.4=1/(3-0.5) and 0.6=1/((1/0.5)-(1/3)).

Presumably, the correct result for the filtered signal will be almost identical to the original signal.

In MATLAB, I can perform the filter by a cascade of second order sections, using the code

n=15; Wn=[0.5 3]/6.25; [z,p,k]=butter(n,Wn); [sos,g] = zp2sos(z,p,k); d=filtfilt(sos,g,x);

or I can perform the filter directly from the transfer function coefficients, using the code

n=15; Wn=[0.5 3]/6.25; [b,a]=butter(n,Wn); d=filtfilt(b,a,x);

The results for the two approaches are as follows:

The upper graph gives the result of the SOS approach, and the lower graph gives the result of the direct approach.

It seems that for quite a significant proportion of the duration of the signal (around the start and end), the SOS method is giving considerably inaccurate results.

Why is the direct approach giving (what I presume to be) a much more accurate result overall than the SOS approach? Should some caveats be added to the general advice that one should always use SOS?

Answer 1

I don't think that has anything to do with SOS vs transfer function representation, but it appears to be an artifact how filtfilt() handles initial conditions. The documentation is annoyingly vague about this and simply says it

filtfilt minimizes start-up and ending transients by matching initial conditions.

If you implement filtfilt() manually you get the expected result: both results are basically correct, but SOS has less numerical noise. See code snippet below.

In general, using filtfilt() is tricky: ideally you have an input signal that has lots of zeros at the beginning and the end, where "lots" means more than the length of the impulse response. Otherwise the results are highly dependent on the initial states of the filters and how you manage them.

Numerically, SOS is always better. Take a look at your "a" coefficients: the ratio of the largest to smallest coefficient is in the order of 1e10. Adding numbers with that different in magnitude creates numerically issues, even if you use double precision. For the SOS representation, that ratios is only about 14, that's about 9 orders of magnitude better!

%% signal
t = (0:.08:100)';
x = sin(2*pi*t);

%% original
n=15;
Wn=[0.5 3]/6.25;
[z,p,k]=butter(n,Wn);
[sos,g] = zp2sos(z,p,k);
d1=filtfilt(sos,g,x);
[b,a]=butter(n,Wn);
d2=filtfilt(b,a,x);

%% manual implementation of filtfilt()
y = filter(b,a,x);
y1 = filter(b,a,flipud(y));
y = g*sosfilt(sos,x);
y2 = g*sosfilt(sos,flipud(y));
plot([y2 y2-y1]); legend('result','difference');