Mercurial > forge
comparison main/signal/sgolayfilt.m @ 0:6b33357c7561 octave-forge
Initial revision
author | pkienzle |
---|---|
date | Wed, 10 Oct 2001 19:54:49 +0000 |
parents | |
children | 1ef8e424b6dc |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:6b33357c7561 |
---|---|
1 ## Copyright (C) 2001 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 ## y = sgolayfilt (x, p, n) | |
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. | |
21 ## | |
22 ## y = sgolayfilt (x, F) | |
23 ## Smooth the data in x with smoothing filter F computed by sgolay. | |
24 ## | |
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: | |
29 ## | |
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;"); | |
36 ## | |
37 ## See also: sgolay | |
38 | |
39 ## TODO: Doesn't accept "weight vector" as fourth parameter. | |
40 ## TODO: Patch filter.cc so that it accepts matrix arguments | |
41 | |
42 function y = sgolayfilt (x, p, n) | |
43 | |
44 if nargin < 1 || nargin > 3 | |
45 usage("y = sgolayfilt(x,p,n) or y = sgolayfilt(x,F)"); | |
46 endif | |
47 | |
48 if (nargin < 2) | |
49 p = 3; | |
50 endif | |
51 if (nargin == 3) | |
52 F = sgolay(p, n); | |
53 elseif (prod(size(p)) == 1) | |
54 n = p+3-rem(p,2); | |
55 F = sgolay(p, n); | |
56 else | |
57 F = p; | |
58 n = size(F,1); | |
59 if (size(F,1) != size(F,2)) | |
60 error("sgolayfilt(x,F): F is not a Savitzsky-Golay filter set"); | |
61 endif | |
62 endif | |
63 | |
64 transpose = (size(x,1) == 1); | |
65 if (transpose) x = x.'; endif; | |
66 len = size(x,1); | |
67 if (len < n) | |
68 error("sgolayfilt: insufficient data for filter"); | |
69 endif | |
70 | |
71 ## The first k rows of F are used to filter the first k points | |
72 ## of the data set based on the first n points of the data set. | |
73 ## The last k rows of F are used to filter the last k points | |
74 ## of the data set based on the last n points of the dataset. | |
75 ## The remaining data is filtered using the central row of F. | |
76 k = floor(n/2); | |
77 z = filter(F(k+1,:), 1, x); | |
78 y = [ F(1:k,:)*x(1:n,:) ; z(n:len,:) ; F(k+2:n,:)*x(len-n+1:len,:) ]; | |
79 | |
80 if (transpose) y = y.'; endif | |
81 | |
82 endfunction | |
83 | |
84 %!demo | |
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(""); | |
93 %! | |
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(""); | |
100 | |
101 %!demo | |
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(""); | |
111 %! | |
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(""); | |
118 %! | |
119 %! oneplot(); |