Why does the “filtfilt” function cause a large transient response?

I have observed that the "filtfilt" function exhibits undesired behavior when I use it with IIR bandpass filters that have steep transition bands. Specifically, the output signal shows an excessive transient response.

Here is a code to generate an example

%% filter design
bp = designfilt('bandpassiir', 'StopbandFrequency1', 50, 'PassbandFrequency1', 250,...
  'PassbandFrequency2', 3600, 'StopbandFrequency2', 3700, 'StopbandAttenuation1', 30,...
  'PassbandRipple', 0.1, 'StopbandAttenuation2', 30, 'SampleRate', 8000, 'DesignMethod', 'cheby2');

%% apply filtfilt to a random (white) long input signal; output signal shows an undesirable transient
rng(100);
x = randn(2^20,1);
ysos = filtfilt(bp,x);

% plot result with sos form
figure(1)
plot(ysos)
axis padded

The above code generates the following results, which shows an excessive transient response at a level of 1 × 10^7 at the beginning and the end of filtered signal.

Answer 1

[MathWorks Support Team]

The issue lies in the ordering of the SOS sections in this filter design - specifically, the first two sections of the SOS matrix have significantly larger gains compared to the rest, causing “filtfilt” transients to become unstable.

To resolve this issue, you could consider these two workarounds

Workaround 1: Reorder SOS Sections

The first workaround is to reorder SOS sections to avoid large gain sections. This can be done by calling “reorder” function. Specifically, if you have the DSP System Toolbox, set “SystemObject” option as “true” when designing the filter with “designfilt”. This ensures the digital filter returns a “dsp.SOSFilter” object. Then, call the “reorder” function to avoid leading SOS sections with large gains.

Please refer to the sample code below

%% WORKAROUND 1: reorder SOS sections
clear
% create a random long input signal
rng(100);
x = randn(2^20,1);

% filter before reordering
bp = designfilt('bandpassiir', 'StopbandFrequency1', 50, 'PassbandFrequency1', 250,...
  'PassbandFrequency2', 3600, 'StopbandFrequency2', 3700, 'StopbandAttenuation1', 30,...
  'PassbandRipple', 0.1, 'StopbandAttenuation2', 30, 'SampleRate', 8000, 'DesignMethod', 'cheby2',...
  'SystemObject',true);
filterAnalyzer(bp,CTFAnalysisMode="individual")

% filter after reordering
reorder(bp,'bandpass')
filterAnalyzer(bp,CTFAnalysisMode="individual")

NUM = bp.Numerator;
DEN = bp.Denominator;
g = bp.ScaleValues;

ysos = filtfilt({NUM,DEN,g},x);

% plot result with SOS form
figure(1)
plot(ysos)
axis padded

Workaround 2: Using Zero-Pole-Gain Form then SOS Form

The second workaround is to convert the filter to zero-pole-gain form and then to cascaded transfer function form. This will make sure the sections are ordered to make the filter more numerically stable.

Please refer to the sample code below

%% UPDATED WORKAGROUND 2: using zero-pole-gain form and then SOS form
clear
% create a random long input signal
rng(100);
x = randn(2^20,1);

% degisn filter
bp = designfilt('bandpassiir', 'StopbandFrequency1', 50, 'PassbandFrequency1', 250,...
  'PassbandFrequency2', 3600, 'StopbandFrequency2', 3700, 'StopbandAttenuation1', 30,...
  'PassbandRipple', 0.1, 'StopbandAttenuation2', 30, 'SampleRate', 8000, 'DesignMethod', 'cheby2');

[z,p,k] = zpk(bp);
[sos,g] = zp2sos(z,p,k); 

ysos = filtfilt(sos,g,x);

[b,a] = bp.tf;
ytf = filtfilt(b,a,x);

% plot result with ZPK approach
figure(1)
plot(ysos)
axis padded

% plot result with TF form
figure(2)
plot(ytf)
axis padded

In both workarounds, the filtered signal no longer have large transients.

Reordering involves examining the zeros and poles of the filter and placing them in an optimal order to avoid sections with very large or very small gains, effectively equalizing all filter gains. Therefore, the pole/zero ordering and scale value recomputing achieve this purpose. Reordering can be further optimized if the filter response is known, which is why the response name was passed to the “reorder” function in workaround1. Once section gains are equalized, “filtfilt” can better estimate the transients.

The “reorder” function converts the SOS sections to zeros, poles, and gains, then reorders them using “zp2sos(z, p, k, 'down')”. The “zp2sos” function reorders poles following the algorithm described in this reference page.

Afterwards, the reorder function further optimizes the sections based on the filter type. You can find a detailed description in the reference below:

D. Schlichtharle. Digital Filters: Basics and Design. Springer-Verlag, Berlin, 2000.

[Paul]

Further discussion at the originating question