1 ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
2 ## Copyright (C) 2007 Francesco Potortì <pot@gnu.org>
3 ## Copyright (C) 2008 Luca Citi <lciti@essex.ac.uk>
5 ## This program is free software; you can redistribute it and/or modify it under
6 ## the terms of the GNU General Public License as published by the Free Software
7 ## Foundation; either version 3 of the License, or (at your option) any later
10 ## This program is distributed in the hope that it will be useful, but WITHOUT
11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15 ## You should have received a copy of the GNU General Public License along with
16 ## this program; if not, see <http://www.gnu.org/licenses/>.
18 ## usage: y = filtfilt(b, a, x)
20 ## Forward and reverse filter the signal. This corrects for phase
21 ## distortion introduced by a one-pass filter, though it does square the
22 ## magnitude response in the process. That's the theory at least. In
23 ## practice the phase correction is not perfect, and magnitude response
24 ## is distorted, particularly in the stop band.
27 ## [b, a]=butter(3, 0.1); % 10 Hz low-pass filter
28 ## t = 0:0.01:1.0; % 1 second sample
29 ## x=sin(2*pi*t*2.3)+0.25*randn(size(t)); % 2.3 Hz sinusoid+noise
30 ## y = filtfilt(b,a,x); z = filter(b,a,x); % apply filter
31 ## plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')
34 ## TODO: (pkienzle) My version seems to have similar quality to matlab,
35 ## but both are pretty bad. They do remove gross lag errors, though.
38 function y = filtfilt(b, a, x)
42 rotate = (size(x,1)==1);
43 if rotate, # a row vector
44 x = x(:); # make it a column vector
54 if la < n, a(n) = 0; end
55 if lb < n, b(n) = 0; end
57 ## Compute a the initial state taking inspiration from
58 ## Likhterov & Kopeika, 2003. "Hardware-efficient technique for
59 ## minimizing startup transients in Direct Form II digital filters"
60 kdc = sum(b) / sum(a);
61 if (abs(kdc) < inf) # neither NaN nor +/- Inf
62 si = fliplr(cumsum(fliplr(b - kdc * a)));
64 si = zeros(size(a)); # fall back to zero initialization
68 for (c = 1:size(x,2)) # filter all columns, one by one
69 v = [2*x(1,c)-x((lrefl+1):-1:2,c); x(:,c);
70 2*x(end,c)-x((end-1):-1:end-lrefl,c)]; # a column vector
72 ## Do forward and reverse filtering
73 v = filter(b,a,v,si*v(1)); # forward filter
74 v = flipud(filter(b,a,flipud(v),si*v(end))); # reverse filter
75 y(:,c) = v((lrefl+1):(lx+lrefl));
78 if (rotate) # x was a row vector
79 y = rot90(y); # rotate it back
86 %!error filtfilt (1, 2, 3, 4);
91 %! [b,a] = butter(10, [.2, .25]);
92 %! yfb = filtfilt(b, a, r);
93 %! assert (size(r), size(yfb));
94 %! assert (mean(abs(yfb)) < 1e3);
95 %! assert (mean(abs(yfb)) < mean(abs(r)));
96 %! ybf = fliplr(filtfilt(b, a, fliplr(r)));
97 %! assert (mean(abs(ybf)) < 1e3);
98 %! assert (mean(abs(ybf)) < mean(abs(r)));
102 %! r = randn(1,1000);
103 %! s = 10 * sin(pi * 4e-2 * (1:length(r)));
104 %! [b,a] = cheby1(2, .5, [4e-4 8e-2]);
105 %! y = filtfilt(b, a, r+s);
106 %! assert (size(r), size(y));
107 %! assert (mean(abs(y)) < 1e3);
108 %! assert (corr(s(250:750), y(250:750)) > .95)
109 %! [b,a] = butter(2, [4e-4 8e-2]);
110 %! yb = filtfilt(b, a, r+s);
111 %! assert (mean(abs(yb)) < 1e3);
112 %! assert (corr(y, yb) > .99)
116 %! r = randn(1,1000);
117 %! s = 10 * sin(pi * 4e-2 * (1:length(r)));
118 %! [b,a] = butter(2, [4e-4 8e-2]);
119 %! y = filtfilt(b, a, [r.' s.']);
120 %! yr = filtfilt(b, a, r);
121 %! ys = filtfilt(b, a, s);
122 %! assert (y, [yr.' ys.']);
123 %! y2 = filtfilt(b.', a.', [r.' s.']);