Why doesn't 'filtfilt' match Gustafsson's plot?
4 views (last 30 days)
Show older comments
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 Comments
Answers (0)
See Also
Categories
Find more on Filter Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!