octave_packages/signal-1.1.3/sgolay.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 ## F = sgolay (p, n [, m [, ts]])
  18 ##   Computes the filter coefficients for all Savitzsky-Golay smoothing
  19 ##   filters of order p for length n (odd). m can be used in order to
  20 ##   get directly the mth derivative. In this case, ts is a scaling factor.
  21 ##
  22 ## The early rows of F smooth based on future values and later rows
  23 ## smooth based on past values, with the middle row using half future
  24 ## and half past.  In particular, you can use row i to estimate x(k)
  25 ## based on the i-1 preceding values and the n-i following values of x
  26 ## values as y(k) = F(i,:) * x(k-i+1:k+n-i).
  27 ##
  28 ## Normally, you would apply the first (n-1)/2 rows to the first k
  29 ## points of the vector, the last k rows to the last k points of the
  30 ## vector and middle row to the remainder, but for example if you were
  31 ## running on a realtime system where you wanted to smooth based on the
  32 ## all the data collected up to the current time, with a lag of five
  33 ## samples, you could apply just the filter on row n-5 to your window
  34 ## of length n each time you added a new sample.
  35 ##
  36 ## Reference: Numerical recipes in C. p 650
  37 ##
  38 ## See also: sgolayfilt
  39
  40 ## Based on smooth.m by E. Farhi <manuf@ldv.univ-montp2.fr>
  41
  42 function F = sgolay (p, n, m = 0, ts = 1)
  43
  44   if (nargin < 2 || nargin > 4)
  45     print_usage;
  46   elseif rem(n,2) != 1
  47     error ("sgolay needs an odd filter length n");
  48   elseif p >= n
  49     error ("sgolay needs filter length n larger than polynomial order p");
  50   else
  51     if length(m) > 1, error("weight vector unimplemented"); endif
  52
  53     ## Construct a set of filters from complete causal to completely
  54     ## noncausal, one filter per row.  For the bulk of your data you
  55     ## will use the central filter, but towards the ends you will need
  56     ## a filter that doesn't go beyond the end points.
  57     F = zeros (n, n);
  58     k = floor (n/2);
  59     for row = 1:k+1
  60       ## Construct a matrix of weights Cij = xi ^ j.  The points xi are
  61       ## equally spaced on the unit grid, with past points using negative
  62       ## values and future points using positive values.
  63       C = ( [(1:n)-row]'*ones(1,p+1) ) .^ ( ones(n,1)*[0:p] );
  64       ## A = pseudo-inverse (C), so C*A = I; this is constructed from the SVD
  65       A = pinv(C);
  66       ## Take the row of the matrix corresponding to the derivative
  67       ## you want to compute.
  68       F(row,:) = A(1+m,:);
  69     end
  70     ## The filters shifted to the right are symmetric with those to the left.
  71     F(k+2:n,:) = (-1)^m*F(k:-1:1,n:-1:1);
  72
  73   endif
  74   F =  F * ( prod(1:m) / (ts^m) );
  75 endfunction
  76
  77 %!test
  78 %! N=2^12;
  79 %! t=[0:N-1]'/N;
  80 %! dt=t(2)-t(1);
  81 %! w = 2*pi*50;
  82 %! offset = 0.5; # 50 Hz carrier
  83 %! # exponential modulation and its derivatives
  84 %! d = 1+exp(-3*(t-offset));
  85 %! dd = -3*exp(-3*(t-offset));
  86 %! d2d = 9*exp(-3*(t-offset));
  87 %! d3d = -27*exp(-3*(t-offset));
  88 %! # modulated carrier and its derivatives
  89 %! x = d.*sin(w*t);
  90 %! dx = dd.*sin(w*t) + w*d.*cos(w*t);
  91 %! d2x = (d2d-w^2*d).*sin(w*t) + 2*w*dd.*cos(w*t);
  92 %! d3x = (d3d-3*w^2*dd).*sin(w*t) + (3*w*d2d-w^3*d).*cos(w*t);
  93 %!
  94 %! y = sgolayfilt(x,sgolay(8,41,0,dt));
  95 %! assert(norm(y-x)/norm(x),0,5e-6);
  96 %!
  97 %! y = sgolayfilt(x,sgolay(8,41,1,dt));
  98 %! assert(norm(y-dx)/norm(dx),0,5e-6);
  99 %!
 100 %! y = sgolayfilt(x,sgolay(8,41,2,dt));
 101 %! assert(norm(y-d2x)/norm(d2x),0,1e-5);
 102 %!
 103 %! y = sgolayfilt(x,sgolay(8,41,3,dt));
 104 %! assert(norm(y-d3x)/norm(d3x),0,1e-4);