octave_packages/signal-1.1.3/sgolayfilt.m

   1 ## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
   2 ## Copyright (C) 2004 Pascal Dupuis <Pascal.Dupuis@esat.kuleuven.ac.be>
   3 ##
   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
   7 ## version.
   8 ##
   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
  12 ## details.
  13 ##
  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/>.
  16
  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.
  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: Patch filter.cc so that it accepts matrix arguments
  40
  41 function y = sgolayfilt (x, p = 3, n, m = 0, ts = 1)
  42
  43   if nargin < 1 || nargin > 5
  44     print_usage;
  45   endif
  46
  47   if (nargin >= 3)
  48     F = sgolay(p, n, m, ts);
  49   elseif (prod(size(p)) == 1)
  50     n = p+3-rem(p,2);
  51     F = sgolay(p, n);
  52   else
  53     F = p;
  54     n = size(F,1);
  55     if (size(F,1) != size(F,2))
  56       error("sgolayfilt(x,F): F is not a Savitzsky-Golay filter set");
  57     endif
  58   endif
  59
  60   transpose = (size(x,1) == 1);
  61   if (transpose) x = x.'; endif;
  62   len = size(x,1);
  63   if (len < n)
  64     error("sgolayfilt: insufficient data for filter");
  65   endif
  66
  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.
  75
  76   k = floor(n/2);
  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,:) ];
  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("");