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

## CROSS-REFERENCE INFORMATION

This function calls:
 error Display message and abort function. lu LU factorization. sprintf Write formatted data to string. sub2ind Linear index from multiple subscripts. warning Display warning message; disable or enable warning messages.
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 %
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