Matlab `filtfilt` provides excessive transient

Question

I observed that filtfilt suffers from an undesired behavior when I provide IIR bandpass filters having steep transition bands. Specifically, the output signal exhibits excessive transient response; nevertheless, this behaviour does not emerge if I use a 'home-made' version of zero phase filtering based on filter. I guess that it is not a numerical issue involving the filter coefficients or structure, nor the filter function.

% design bandpass filter having transition bandwidth of 200 Hz (Fs = 8000)
bp = designfilt('bandpassiir', 'StopbandFrequency1', 50, 'PassbandFrequency1', 250,...
    'PassbandFrequency2', 3600, 'StopbandFrequency2', 3800, 'StopbandAttenuation1', 30,...
    'PassbandRipple', 0.1, 'StopbandAttenuation2', 30, 'SampleRate', 8000, 'DesignMethod', 'cheby2');
% check stability
assert(isstable(bp),'Unstable filter');
% apply filtfilt to a random (white) long input signal; output signal shows an undesirable transient
x = randn(2^20,1);
y = filtfilt(bp,x);
% apply 'home-made' filtfilt to the same input; output signal shows a more acceptable transient
y2 = flipud(filter(bp,flipud(filter(bp,x))));
% compare effects
figure; semilogy(abs(y-y2));

As a rule of thumb, the effect grows as the transition bands get narrower, while it tends to vanish as they get broader.

Where is the problem? Have I missed some recommendations or hints in the function's help?

Answer 1

It must have to do with the initial conditions used by the function filtfilt.m. The idea is to match initial conditions in a way such that startup and end transients are minimized. This, however, doesn't always seem to work, and it appears that for your filter specifications it actually does more harm than good. As far as I know there is no way to tell filtfilt.m not to mess with initial conditions. What remains in this case is to simply implement a straightforward version, just as you did.

Also take a look at this related answer.

EDIT:

I tried to reproduce your results with Octave, but I don't get any of the effects that you describe. I use a $10^{th}$ order Chebyshev 2 filter (which is the lowest order satisfying your specs), and I used a modified version of filtfilt allowing me to both switch on and off the use of initial conditions and the extrapolation of the input signal. There are two explanations for this: either your version of Matlab's filtfilt has a bug, or your filter is substantially different from mine and it introduces numerical problems (due to poles very close to the unit circle), which happen to manifest themselves when the input signal is extrapolated, and/or when special initial conditions are used (as is the case with filtfilt). If you provide your filter coefficients we can probably solve this mystery.

EDIT 2:

Tendero kindly reproduced the filter coefficients using the Matlab commands provided in the question:

b = [0.682585947577945 -0.012906667355225 -3.354515884719406 0.025370108447479 6.651816094305592 -0.000000000000000 -6.651816094305592 -0.025370108447479 3.354515884719406 0.012906667355225 -0.682585947577945];
a = [1.000000000000000 -0.194376467633997 -4.144930396202897 0.592472615834755 7.013126532356848 -0.691483043645839 -6.033921428993976 0.364443043103532 2.633924737223088 -0.072944768631229 -0.465923573748158];

I used these coefficients in Octave with the function filtfilt.m with a random input signal x = randn(2^20,1) just like in the question. As mentioned above, I used a modified version of filtfilt.m to be able to switch on and off the different measures for reducing transients:

1. no extrapolation of the input signal, zero initial conditions
2. no extrapolation of the input signal, optimal initial conditions (as implemented in the standard version of filtfilt.m)
3. extrapolation of the input signal, zero initial conditions
4. extrapolation of the input signal, optimal initial conditions

The last combination is the processing implemented in the standard version of filtfilt.m.

The plot below shows the first and last $50$ samples of the $4$ different output signals. For all other samples the signals are virtually identical.

As expected, one can see slight differences during startup and at the end, but these differences are very small, and the problem of the OP cannot be reproduced using the original filter coefficients and the Octave function filtfilt.m. There must be either a bug in the Matlab version of filtfilt.m used by the OP, or another problem that is neither related to the filter coefficients nor to the function filtfilt.m.

Answer 2

Looking at the script filtfilt.m, there are two things that differ from the straightforward implementation you did in one single line.

First of all, there is what Matt stated in his answer. The function filtfilt uses initial conditions when it calls the function filter (you can see Mathworks' documentation).

The other difference is that when filtfilt calls filter, it does not use x as input but another array created within the function. I don't understand why or what that signal is, but I believe that's the reason why we are getting that weird peak. In Octave's version, what filtfilt does is the following (putting filter delays initialization to zero, i.e. the modified filtfilt):

lx = size(x,1);
lb = length(b);
la = length(a);
n = max(lb, la);
lrefl = 3 * (n - 1);
if la < n, a(n) = 0; endif
if lb < n, b(n) = 0; endif

v = [ 2*x(1) - x((lrefl+1):-1:2);
      x(:);
      2*x(end) - x((end-1):-1:end-lrefl)];

v = filter(b,a,v);                   # forward filter
v = flipud(filter(b,a,flipud(v)));   # reverse filter
y = v((lrefl+1):(lx+lrefl));

Maybe someone can help us understand what that v represents and why it is useful, but that's the only difference I can spot between filtfilt and your 'homemade' implementation, and so what I expect to be the reason for the undesired transient.