1 ## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
2 ## Copyright (C) 2004 Pascal Dupuis <Pascal.Dupuis@esat.kuleuven.ac.be>
4 ## This program is free software; you can redistribute it and/or modify it under
5 ## the terms of the GNU General Public License as published by the Free Software
6 ## Foundation; either version 3 of the License, or (at your option) any later
9 ## This program is distributed in the hope that it will be useful, but WITHOUT
10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 ## You should have received a copy of the GNU General Public License along with
15 ## this program; if not, see <http://www.gnu.org/licenses/>.
17 ## y = sgolayfilt (x, p, n [, m [, ts]])
18 ## Smooth the data in x with a Savitsky-Golay smoothing filter of
19 ## polynomial order p and length n, n odd, n > p. By default, p=3
20 ## and n=p+2 or n=p+3 if p is even.
22 ## y = sgolayfilt (x, F)
23 ## Smooth the data in x with smoothing filter F computed by sgolay.
25 ## These filters are particularly good at preserving lineshape while
26 ## removing high frequency squiggles. Particularly, compare a 5 sample
27 ## averager, an order 5 butterworth lowpass filter (cutoff 1/3) and
28 ## sgolayfilt(x, 3, 5), the best cubic estimated from 5 points:
30 ## [b, a] = butter(5,1/3);
31 ## x=[zeros(1,15), 10*ones(1,10), zeros(1,15)];
32 ## plot(sgolayfilt(x),"r;sgolayfilt;",...
33 ## filtfilt(ones(1,5)/5,1,x),"g;5 sample average;",...
34 ## filtfilt(b,a,x),"c;order 5 butterworth;",...
35 ## x,"+b;original data;");
39 ## TODO: Patch filter.cc so that it accepts matrix arguments
41 function y = sgolayfilt (x, p = 3, n, m = 0, ts = 1)
43 if nargin < 1 || nargin > 5
48 F = sgolay(p, n, m, ts);
49 elseif (prod(size(p)) == 1)
55 if (size(F,1) != size(F,2))
56 error("sgolayfilt(x,F): F is not a Savitzsky-Golay filter set");
60 transpose = (size(x,1) == 1);
61 if (transpose) x = x.'; endif;
64 error("sgolayfilt: insufficient data for filter");
67 ## The first k rows of F are used to filter the first k points
68 ## of the data set based on the first n points of the data set.
69 ## The last k rows of F are used to filter the last k points
70 ## of the data set based on the last n points of the dataset.
71 ## The remaining data is filtered using the central row of F.
72 ## As the filter coefficients are used in the reverse order of what
73 ## seems the logical notation, reverse F(k+1, :) so that antisymmetric
74 ## sequences are used with the right sign.
77 z = filter(F(k+1,n:-1:1), 1, x);
78 y = [ F(1:k,:)*x(1:n,:) ; z(n:len,:) ; F(k+2:n,:)*x(len-n+1:len,:) ];
80 if (transpose) y = y.'; endif
85 %! [b, a] = butter(5,1/3);
86 %! x=[zeros(1,15), 10*ones(1,10), zeros(1,15)];
87 %! subplot(121); title("boxcar");
88 %! axis([1 40 -2 15]);
89 %! plot(sgolayfilt(x),"r;sgolay(3,5);",...
90 %! filtfilt(ones(1,5)/5,1,x),"g;5 sample average;",...
91 %! filtfilt(b,a,x),"c;order 5 butterworth;",...
92 %! x,"+b;original data;"); title("");
94 %! x=x+randn(size(x))/2;
95 %! subplot(122); title("boxcar+noise");
96 %! plot(sgolayfilt(x,3,5),"r;sgolay(3,5);",...
97 %! filtfilt(ones(1,5)/5,1,x),"g;5 sample average;",...
98 %! filtfilt(b,a,x),"c;order 5 butterworth;",...
99 %! x,"+b;original data;"); title("");
102 %! [b, a] = butter(5,1/3);
103 %! t = 0:0.01:1.0; % 1 second sample
104 %! x=cos(2*pi*t*3); % 3 Hz sinusoid
105 %! subplot(121); title("sinusoid");
106 %! axis([0 1 -1.5 2.5]);
107 %! plot(t,sgolayfilt(x,3,5),"r;sgolay(3,5);",...
108 %! t,filtfilt(ones(1,5)/5,1,x),"g;5 sample average;",...
109 %! t,filtfilt(b,a,x),"c;order 5 butterworth;",...
110 %! t,x,"+b;original data;"); title("");
112 %! x=x+0.2*randn(size(x)); % signal+noise
113 %! subplot(122); title("sinusoid+noise");
114 %! plot(t,sgolayfilt(x',3,5),"r;sgolay(3,5);",...
115 %! t,filtfilt(ones(1,5)/5,1,x),"g;5 sample average;",...
116 %! t,filtfilt(b,a,x),"c;order 5 butterworth;",...
117 %! t,x,"+b;original data;"); title("");