filtfilt provides excessive transient
73 views (last 30 days)
Show older comments
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
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))
[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.
Answers (1)
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.
0 Comments
See Also
Categories
Find more on Digital Filter Analysis in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!