filtfilt provides excessive transient

9 次查看（过去 30 天）
alessandro lapini
评论： Vivian 2024-9-6
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 个评论
显示 1更早的评论
Paul
Paul 2024-8-31
编辑：Paul 2024-8-31
Ran in:
@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
Vivian 2024-9-6
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

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

类别

在 Help Center 和 File Exchange 中查找有关 Digital Filter Analysis 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by