filtfilt provides excessive transient

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 zerophase 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', 3700, '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 accptable 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?

[Vivian]

I realize that this post is several years old, but since it’s still getting views…. I was looking into the OP’s question and noticed that I can, in fact, duplicate his result. However, if instead of:

y = filtfilt(bp,x);

I do:

[b,a] = sos2tf(bp.Coefficients);
y = filtfilt(b,a,x);

I do not get the large transient.

The OP’s filter does have poles very close to the unit circle. However, one would think that using filtfilt with bp (sos) would be less susceptible to numerical issues. So blaming this all on the filter doesn’t seem quite right.

I have observed the same behavior with other narrow transition band IIR filters (i.e., large transient with sos; none with b,a). Could there be a bug in the sos implementation of Matlab’s filtfilt?

See also this related post: https://dsp.stackexchange.com/questions/46945/why-is-the-second-order-section-cascade-approach-giving-worse-results

[Paul]

@alessandro lapini and @Vivian

These results are very interesting because the sos form is typically recommended over the tf form (which I guess is why designfilt returns an sos filter), and filtfilt is often recommended to perform filtering (even though it doesn't really perform zero-phase filtering).

Illustrate the results described above

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

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

Applying filtfilt on the sos form yields enormous transients at the ends of the output

figure
plot(ysos)
axis padded

whereas applying filtfilt to the tf form does not have those large transients

figure
plot(ytf)
axis padded

Plot the difference, I'm a bit surprised the difference isn't smaller between the transients

figure
semilogy(abs(ysos-ytf))
axis padded

Comparing the frequency responses shows that filtfilt with the tf form essentially yields the expected result, whereas the output with the sos form is wildly off. I suspect, but don't know for sure, that this result is driven by those wild transients.

[hpb,wn] = freqz(bp,4096);
hx       = freqz(x,1,wn);
hsos     = freqz(ysos,1,wn);
htf      = freqz(ytf,1,wn);

figure
plot(wn,db([abs(hpb).^2, abs(hsos./hx), abs(htf./hx)])),grid
legend('hpb','sos','tf')

With the tf form, filtfilt makes one forward-backward pass over the input, with some extra sauce "to minimize startup and ending transients" (thought the way the doc page filtfilt describes this is incorrect).

With the sos form, filtfilt makes one forward-backward pass for each section (eight in this case) and applies the extra sauce to "minimize ... transients" on each each pass for each section. Unclear why filtfilt operates this way, as opposed to do just doing a single forward-backward pass, or if it's implemented correctly. Whether or not the approach used in filtfilt for the sos form is even mathematically justified is an open question in my mind, particularly considering these results.

I guess that @Anh Tran is correct that the result with the sos form is driven by how filtfilt tries to "minimize ... transients," though I'm not convinced these results have anything do with the how close the filter is to being unstable (which is based on the poles, not the zeros) or "numerical conditioning" for which the sos form is supposed to be superior.

Also, filtfilt uses a different algorithm for inputs that have less than 10000 elements, though I presume that it's functionally equivalent to the algorithm used for longer inputs.

The difference between the sos and tf result becomes more obvious for shorter inputs for this filter

xsmall = randn(11000,1);  % more than 1e4 points to ensure filtfilt uses the same algorithm as above
ysos = filtfilt(bp,xsmall);
ytf  = filtfilt(b,a,xsmall);
figure
semilogy(abs(ysos-ytf))

Am curious about whether or not bandpass exhibits the same behavior

[y,d] = bandpass(x,[250 3600],8000,ImpulseResponse = "iir");
hd = freqz(d,wn);

The filters aren't quite the same

figure
subplot(211);plot(wn,abs(hpb),wn,abs(hd)),grid
subplot(212);plot(wn,180/pi*angle(hpb),wn,180/pi*angle(hd)),grid

and the filter obtained from bandpass doesn't result in the transients, even though bandpass calls filtfilt for an iir filter (according to its doc page).

figure
plot(y)

@Vivian, I think you should contact Tech Support about this issue that you've found regarding the difference in output from filtfilt between the sos and tf forms. If you do, please post back here with a summary of their response.

[Vivian]

@Paul, thanks for providing the worked out examples. I reached out to Tech Support today and will let you know if/when they get back to me.

Answer 1

[Anh Tran]

The transients observed are due to a combination of using a marginally stable filter coupled with the initial condition matching performed by "filtfilt" to minimize start-up and ending transients, as described in the documentation. This is the reason behind the difference in behavior between using "filter" and using "filtfilt". It is an inherent limitation of the "filtfilt" function.

Regarding the stability of the filter, you may view the filter's stability using "fvtool", (fvtool(bp) in your case). In the Filter Visualization Tool window, select "Analysis > Pole/Zero Plot", and you will see a discrete system pole/zero map If you notice, your zeroes are very close to the edge of the unit circle, indicating the marginally stable behavior of this filter. When this type of filter is combined with the initial condition matching, the numerical conditioning yields the transient results that you observed.