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
Display message and abort function. | |
LU factorization. | |
Write formatted data to string. | |
Linear index from multiple subscripts. | |
Display warning message; disable or enable warning messages. |
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
Comments or corrections to the WAFO group