Why doesn't 'filtfilt' match Gustafsson's plot?

4 次查看(过去 30 天)
The MATLAB filtfilt() function references F. Gustafsson, 1996:
Gustafsson, Fredrik. "Determining the initial states in forward-backward filtering." IEEE Transactions on signal processing 44.4 (1996): 988-992.
Figure 1 of that paper shows some example plots of forward/backware filtering a white noise sequence. I can reproduce the optimized plot (lower-left) and the zero initial conditions plot (upper-right), but not the example he cites using MATLAB's filtfilt() function (upper-left). Now, it's possible that the implementation of filtfilt() has changed sometime in the last 20+ years -- is this related to his footnote #2 about v4.1 and earlier having a bug, or did an older implementation allow passing initial conditions? If not the case, then why does the following code not reproduce the upper-left plot?
figure(1);
clf(1);
% Initialize a white-noise sequence of length 200 with seed 0
N = 200;
rng(0,'v4');
u = randn(N,1);
% Design a band-pass Butterworth filter of tenth order with
% (normalized) cut-off frequencies of 0.2 and 0.25
n = 10;
[b,a] = tf(designfilt('bandpassiir', ...
'DesignMethod', 'butter', ...
'FilterOrder', n, ...
'HalfPowerFrequency1', 0.20, ...
'HalfPowerFrequency2', 0.25 ));
% Filter the signal using Matlab's FILTFILT() function
Y_fb = filtfilt(b,a,u);
Y_bf = flip(filtfilt(b,a,flip(u)));
% This plot should match upper-left subplot of Gustafsson, Fig. 1
subplot(1,2,1);
box('on');
xlim([1 N]);
ylim(0.6*[-1 1]);
if exist('upper_left.jpg','file')
imagesc(xlim,ylim,flip(imread('upper_left.jpg'),1));
set(gca,'ydir','normal');
end
hold('on');
plot(Y_fb, 'Color', 'b', 'LineStyle', '--', 'LineWidth', 1);
plot(Y_bf, 'Color', 'r', 'LineStyle', '--', 'LineWidth', 2);
% Filter the signal using zero initial conditions
Y_fb_0 = flip(filter(b,a,flip(filter(b,a,u,zeros(1,n)))));
Y_bf_0 = filter(b,a,flip(filter(b,a,flip(u),zeros(1,n))));
% This plot should match upper-right subplot of Gustafsson, Fig. 1
subplot(1,2,2);
box('on');
xlim([1 N]);
ylim(0.6*[-1 1]);
if exist('upper_right.jpg','file')
imagesc(xlim,ylim,flip(imread('upper_right.jpg'),1));
set(gca,'ydir','normal');
end
hold('on');
plot(Y_fb_0, 'Color', 'b', 'LineStyle', '--', 'LineWidth', 1);
plot(Y_bf_0, 'Color', 'r', 'LineStyle', '--', 'LineWidth', 2);

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Filter Analysis 的更多信息

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by