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

4 ビュー (過去 30 日間)
David De Lorenzo
David De Lorenzo 2020 年 1 月 6 日
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 Center および File ExchangeFilter Analysis についてさらに検索

製品


リリース

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by