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

3 次查看（过去 30 天）

David De Lorenzo 2020-1-6

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/499158-why-doesn-t-filtfilt-match-gustafsson-s-plot

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.

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.54.870&rep=rep1&type=pdf

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);