Why doesn't 'filtfilt' match Gustafsson's plot?
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Filter Analysis finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!