|
0
|
1 ## Copyright (C) 1999 Paul Kienzle |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify |
|
|
4 ## it under the terms of the GNU General Public License as published by |
|
|
5 ## the Free Software Foundation; either version 2 of the License, or |
|
|
6 ## (at your option) any later version. |
|
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, |
|
|
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
11 ## GNU General Public License for more details. |
|
|
12 ## |
|
|
13 ## You should have received a copy of the GNU General Public License |
|
|
14 ## along with this program; if not, write to the Free Software |
|
|
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
|
16 |
|
|
17 ## usage: y = filtfilt(b, a, x) |
|
|
18 ## |
|
|
19 ## Forward and reverse filter the signal. This corrects for phase |
|
|
20 ## distortion introduced by a one-pass filter, though it does square the |
|
|
21 ## magnitude response in the process. That's the theory at least. In |
|
|
22 ## practice the phase correction is not perfect, and magnitude response |
|
|
23 ## is distorted, particularly in the stop band. |
|
|
24 ## |
|
|
25 ## In this version, I zero-pad the end of the signal to give the reverse |
|
|
26 ## filter time to ramp up to the level at the end of the signal. |
|
|
27 ## Unfortunately, the degree of padding required is dependent on the |
|
|
28 ## nature of the filter and not just its order, so this function needs |
|
|
29 ## some work yet. |
|
|
30 ## |
|
|
31 ## Example |
|
|
32 ## [b, a]=butter(3, 0.1); % 10 Hz low-pass filter |
|
|
33 ## t = 0:0.01:1.0; % 1 second sample |
|
|
34 ## x=sin(2*pi*t*2.3)+0.25*randn(size(t)); % 2.3 Hz sinusoid+noise |
|
|
35 ## y = filtfilt(b,a,x); z = filter(b,a,x); % apply filter |
|
|
36 ## plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;') |
|
|
37 |
|
|
38 ## Changelog: |
|
|
39 ## 2000 02 pkienzle@kienzle.powernet.co.uk |
|
|
40 ## - pad with zeros to load up the state vector on filter reverse. |
|
|
41 ## - add example |
|
|
42 |
|
|
43 ## TODO: In Matlab filtfilt `reduces filter startup transients by carefully |
|
|
44 ## TODO: choosing initial conditions, and by prepending onto the input |
|
|
45 ## TODO: sequence a short, reflected piece of the input sequence'. |
|
|
46 ## TODO: Once filtic is written, use that here. |
|
|
47 ## TODO: My version seems to have similar quality to matlab, but both are |
|
|
48 ## TODO: pretty bad. They do remove gross lag errors, though. |
|
|
49 ## TODO: Note that if x is really long, it might be worth doing |
|
|
50 ## TODO: the zero padding as a separate call to filter so that the |
|
|
51 ## TODO: vector never has to be copied. E.g., |
|
|
52 ## TODO: [y, state] = filter(b,a,x); |
|
|
53 ## TODO: tail = filter(b,a,zeros(1,max(length(b),length(a))),state); |
|
|
54 ## TODO: [tail, state] = filter(b,a,flipXX(tail)); |
|
|
55 ## TODO: y = flipXX(filter(b,a,flipXX(y), state)); |
|
|
56 ## TODO: Don't know for what n this would be faster, if any, but the |
|
|
57 ## TODO: memory saving might be nice. |
|
|
58 |
|
|
59 function y = filtfilt(b, a, x) |
|
|
60 if (nargin != 3) |
|
|
61 usage("y=filtfilt(b,a,x)"); |
|
|
62 end |
|
|
63 |
|
|
64 if (rows(x) == 1) |
|
|
65 y = filter(b,a,[x, zeros(1,2*max(length(a),length(b)))]); |
|
|
66 y = fliplr(filter(b,a,fliplr(y))); |
|
|
67 else |
|
|
68 y = filter(b,a,[x ; zeros(2*max(length(a),length(b)), 1)]); |
|
|
69 y = flipud(filter(b,a,flipud(y))); |
|
|
70 endif |
|
|
71 y = y(1:length(x)); |
|
|
72 endfunction |
