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.

[Paul]

@Vivian

Did you get a response from Tech Support? If so, would you mind sharing it (or at least a summary)?

[alelap]

After 7 years I posted this question, I have never found a test case where I had to use filtfilt instead the 'home-made' version I wrote in the post. I think Mathworks should rethink filtfilt or offer an alternative implementation for zero-phase filtering.

[Paul]

Hi alelap,

I think your homemade version could be improved by zero-padding the input to filter on both passes. I've shown what I mean at this thread, which is a more general dicussion on the "zero-phase" response of filtfilt.

Answer 1

[Wentao]

Hello, this issue has been answered by MathWorks Technical Support Department. Please refer to this link for detailed explaination.

https://www.mathworks.com/matlabcentral/answers/2172783-why-does-the-filtfilt-function-cause-a-large-transient-response

Thank you for bringing this issue to our attention.

[Paul]

@Wentao

I'm going to continue the discussion here because, in my experience, the Mathworks Support Team does not respond to comments on their Q&A.

The linked answer basically states that the issue is with the order of the second order sections, and that order causes the filtfilt method of the digitalFilter object, which calls filtfilt for this IIR filter, to be not as "numerically stable" as one might expect.

I'm afraid I have to disagree with that claim.

The results don't really bear the hallmark of such instability.

The transient in question is a lightly damped sinusoid, which looks very clean.

Furthermore, and more importantly, increasing the nfact parameter in filtfilt completely removes the transient from the output. So it seems that the transient we are seeing is entirely due to the filtfilt implementation (the nfact value is proably a heuristic that doesn't work for this problem?) and not inherently due to the filter structure.

Perhaps you can comment on why filtfilt is implemented as it is, where it yields different outputs for the same input depending on if the filter is implemented in tf or sos form, and now we see that there is also a sensitivity to the ordering of the sos sections.

Wouldn't it be better if filtfilt were implemented by making a single forward and backward pass and simply padding the input signal by an appropriate amount going both ways? Such an approach would, mathematically, be agnostic to the filter implementation and would have the advantage of actually being closer to the zero-phase response that filtfilt purports to achieve (and does not as is apparent for the sos case due to the change in result by reordering the sections).

Having said all of that, I do agree that the zp2sos workaround is a good thing to do regardless, and I was quite surprised to learn that designfilt doesn't already take that approach.

Answer 2

[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.