Clear Filters
Clear Filters

filtfilt provides excessive transient

73 views (last 30 days)
alessandro lapini
alessandro lapini on 2 Jan 2018
Commented: Vivian on 6 Sep 2024 at 1:21
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?
  3 Comments
Paul
Paul on 31 Aug 2024
Edited: Paul on 31 Aug 2024
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
Vivian on 6 Sep 2024 at 1:21
@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.

Sign in to comment.

Answers (1)

Anh Tran
Anh Tran on 4 Jan 2018
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!