Home > wafo > misc > savgol.m

savgol

PURPOSE ^

Savitzky-Golay filter coefficients.

SYNOPSIS ^

c = savgol(nl,nr,m,ld,np)

DESCRIPTION ^

 SAVGOL  Savitzky-Golay filter coefficients. 
  
  CALL:  c = savgol(nl,nr,m,ld,np); 
  
   c     = vector with Savitzky-Golay filter coefficients. 
   nl,nr = number of leftward (past) and rightward (future) data points 
           used, respectively, making the total number of data points 
           used nl+nr+1. (nl+nr>=m) 
   m     = order of the smoothing polynomial, also equal to the highest 
           conserved moment; usual values are m = 2 or m = 4. 
   ld    = order of the derivative desired. (default ld=0 for smoothed 
           function) (ld<=m) 
   np    = length of c  (np>=nr+nl+1) (defalt 2*max(nr,nl)+1)  
   
  The idea of Savitzky-Golay filtering is to find filter coefficients Cn 
  that preserves higher moments, i.e., to approximate the underlying 
  function within the moving window not by a constant (whose estimate is 
  the average), but by a polynomial of higher order, typically quadratic 
  or quartic. For each point Y(i), we least-squares fit a polynomial of 
  order M to all NL+NR+1 points in the moving window and then set YF(i) 
  to the value of this polynomial at position I, i.e.,  
  
              nr 
     YF(i) = sum Cn(n)*Y(i+n) 
            n=-nl 
   
  SAVGOL Returns the filter coefficients C(1:np), in wrap-around order, 
  i.e., C(1) = Cn(0), C(2)=Cn(-1),...C(Nl+1)=Cn(-Nl), and C(Np) = Cn(1), 
  C(Np-1) = Cn(2),...C(Np-Nr) = Cn(Nr), which is consistent for use with 
  fft to perform the convolution. 
  
  Note the restrictions:  np >= nr+nl+1 > m >= ld 
  
  Example: 
    x = linspace(0,1); 
    y = exp(x)+1e-1*randn(size(x)); 
    nl=2;nr=2;m=2;ld=0; 
    c = savgol(nl,nr,m,ld,length(x)); 
    yf = real(ifft(fft(y).*fft(c))); % convolution with pronounced end effects 
    yf2 = convlv(y,c);               % convolution with less end effects 
    plot(x,y,x,yf,x,yf2) 
    ld =1; m =4; 
    c = savgol(nl,nr,m,ld);          % differentiation filter 
    dyf = convlv(y,c)*gamma(ld+1)/(x(2)-x(1))^ld; % Derivative 
    ix = nl+1:length(x)-nr;                       % for all these ix 
    semilogy(x(ix),abs(y(ix)-dyf(ix))) % Error of the derivative 
  
  See also  smooth

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function c = savgol(nl,nr,m,ld,np) 
002 %SAVGOL  Savitzky-Golay filter coefficients. 
003 % 
004 % CALL:  c = savgol(nl,nr,m,ld,np); 
005 % 
006 %  c     = vector with Savitzky-Golay filter coefficients. 
007 %  nl,nr = number of leftward (past) and rightward (future) data points 
008 %          used, respectively, making the total number of data points 
009 %          used nl+nr+1. (nl+nr>=m) 
010 %  m     = order of the smoothing polynomial, also equal to the highest 
011 %          conserved moment; usual values are m = 2 or m = 4. 
012 %  ld    = order of the derivative desired. (default ld=0 for smoothed 
013 %          function) (ld<=m) 
014 %  np    = length of c  (np>=nr+nl+1) (defalt 2*max(nr,nl)+1)  
015 %  
016 % The idea of Savitzky-Golay filtering is to find filter coefficients Cn 
017 % that preserves higher moments, i.e., to approximate the underlying 
018 % function within the moving window not by a constant (whose estimate is 
019 % the average), but by a polynomial of higher order, typically quadratic 
020 % or quartic. For each point Y(i), we least-squares fit a polynomial of 
021 % order M to all NL+NR+1 points in the moving window and then set YF(i) 
022 % to the value of this polynomial at position I, i.e.,  
023 % 
024 %             nr 
025 %    YF(i) = sum Cn(n)*Y(i+n) 
026 %           n=-nl 
027 %  
028 % SAVGOL Returns the filter coefficients C(1:np), in wrap-around order, 
029 % i.e., C(1) = Cn(0), C(2)=Cn(-1),...C(Nl+1)=Cn(-Nl), and C(Np) = Cn(1), 
030 % C(Np-1) = Cn(2),...C(Np-Nr) = Cn(Nr), which is consistent for use with 
031 % fft to perform the convolution. 
032 % 
033 % Note the restrictions:  np >= nr+nl+1 > m >= ld 
034 % 
035 % Example: 
036 %   x = linspace(0,1); 
037 %   y = exp(x)+1e-1*randn(size(x)); 
038 %   nl=2;nr=2;m=2;ld=0; 
039 %   c = savgol(nl,nr,m,ld,length(x)); 
040 %   yf = real(ifft(fft(y).*fft(c))); % convolution with pronounced end effects 
041 %   yf2 = convlv(y,c);               % convolution with less end effects 
042 %   plot(x,y,x,yf,x,yf2) 
043 %   ld =1; m =4; 
044 %   c = savgol(nl,nr,m,ld);          % differentiation filter 
045 %   dyf = convlv(y,c)*gamma(ld+1)/(x(2)-x(1))^ld; % Derivative 
046 %   ix = nl+1:length(x)-nr;                       % for all these ix 
047 %   semilogy(x(ix),abs(y(ix)-dyf(ix))) % Error of the derivative 
048 % 
049 % See also  smooth 
050  
051 % Reference  
052 % William H. Press, Saul Teukolsky,  
053 % William T. Wetterling and Brian P. Flannery (1997) 
054 % "Numerical recipes in Fortran 77", Vol. 1, pp 644-649 
055  
056 % History 
057 % by pab 2000 
058  
059  
060 error(nargchk(3,5,nargin)) 
061 if nargin<5|isempty(np),np=2*max(nl,nr)+1;end 
062 if nargin<4|isempty(ld),ld=0;end 
063  
064 if (np<nl+nr+1) 
065   error(sprintf('np must be larger than nl+nr = %d or nl+nr less than np = %d',nr+nl,np)) 
066 end 
067 if nl<0,  
068   warning('nl must be positive or zero') 
069   nl = max(nl,0); 
070 end 
071 if nr<0,  
072   warning('nr must be positive or zero') 
073   nr = max(nr,0); 
074 end 
075 if ld>m  
076   error(sprintf('m must be larger than ld = %d ',ld)) 
077 end 
078 if nl+nr<m, 
079   error(sprintf('m must be smaller than nl+nr = %d',nl+nr)) 
080 end 
081  
082 asiz = [m+1,m+1]; 
083 a    = zeros(asiz); 
084 for ipj=0:2*m   
085   %Set up the normal equations of the desired least-squares fit.    
086   tmp = sum([1:nr].^ipj) + (ipj==0) + sum([-1:-1:-nl].^ipj);  
087    
088   mm = min(ipj,2*m-ipj); 
089   imj=-mm:2:mm; 
090   ind = sub2ind(asiz, 1+(ipj+imj)/2,1+(ipj-imj)/2); 
091   a(ind)=tmp; 
092   %for imj=-mm:2:mm  
093   %  a(1+(ipj+imj)/2,1+(ipj-imj)/2)=tmp; 
094   %end  
095 end  
096  
097 % Solve them by LU decomposition. 
098 [L,U] = lu(a); 
099  
100 % Right-hand side vector is unit vector,  
101 % depending on which derivative we want. 
102 b       = zeros(m+1,1); 
103 b(ld+1) = 1; 
104  
105 % Backsubstitute, giving one row of the inverse matrix. 
106  
107 d = (U\(L\b)).'; 
108  
109  
110 %Zero the output array (it may be bigger than number of coefficients).  
111 c = zeros(1,np); 
112  
113 if 0, 
114   for k=-nl:nr  
115     %Each Savitzky-Golay coefficient is the dot product 
116     %of powers of an integer with the inverse matrix row. 
117     tmp = sum(d(1:m+1).*[1 k.^(1:m)]); 
118     kk = mod(np-k,np)+1; %Store in wrap-around order. 
119     c(kk) = tmp; 
120   end  
121 else 
122   %Each Savitzky-Golay coefficient is the dot product 
123   %of powers of an integer with the inverse matrix row. 
124   k    = (-nl:nr ).'; 
125   Nk   = length(k); 
126   tmp0 = repmat(k,1,m+1).^repmat(0:m,Nk,1); 
127   tmp  = sum(d(ones(Nk,1),1:m+1).*tmp0,2)'; 
128   kk   = mod(np-k,np)+1; %Store in wrap-around order. 
129   c(kk) = tmp; 
130 end 
131 return 
132  
133

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003