652 Chapter 14.Statistical Description of Data #include <math.h> #include "nrutil.h" void savgol(float c[],int np,int nl,int nr,int ld,int m) Returns in c[1..np],in wrap-around order (N.B.!)consistent with the argument respns in routine convlv,a set of Savitzky-Golay filter coefficients.nl is the number of leftward (past) data points used,while nr is the number of rightward (future)data points,making the total number of data points used nl+nr+1.Id is the order of the derivative desired (e.g.,ld =0 for smoothed function).m is the order of the smoothing polynomial,also equal to the highest conserved moment;usual values are m =2 or m =4. void lubksb(float **a,int n,int *indx,float b[]); void ludcmp(float **a,int n,int *indx,float *d); int imj,ipj,j,k,kk,mm,*indx; float d,fac,sum,**a,*b; (includi 18881892 if (np nl+nr+1 II nl 0 I nr 0 II ld m Il nl+nr m) nrerror("bad args in savgol"); 1-800 indx=ivector(1,m+1); a=matr1x(1,m+1,1,m+1); b=vector(1,m+1); NUMERICAL RECIPES for (ipj=0;ipj<=(m <1);ipj++){ Set up the normal equations of the desired sum=(1pj？0.0:1.0); least-squares fit. for (k=1;k<=nr;k++)sum +pow((double)k,(double)ipj); (North for (k=1;k<=nl;k++)sum +pow((double)-k,(double)ipj); mm=IMIN(ipj,2*m-ipj); for (imj =-mm;imj<=mm;imj+=2)a[1+(ipj+imj)/2][1+(ipj-imj)/2]=sum; Americ computer Press. ART ludcmp(a,m+1,indx,&d); Solve them:LU decomposition. for(j=1;j<=m+1;j++)b[j]=0.0; b[1d+1]=1.0; 9 Programs Right-hand side vector is unit vector,depending on which derivative we want. lubksb(a,m+1,indx,b); Get one row of the inverse matrix. for (kk=1;kk<=np;kk++)c[kk]=0.0; Zero the output array (it may be bigger than SCIENTIFIC for (k =-nl;k<=nr;k++){ number of coefficients). sum=b[1]; Each Savitzky-Golay coefficient is the dot fac=1.0; product of powers of an integer with the for (mm=1;mm<=m;mm++)sum +b[mm+1]*(fac *k);inverse matrix row. kk=((np-k)np)+1; Store in wrap-around order. c[kk]=sum; COMPUTING(ISBN 198918920 free_vector(b,1,m+1); free_matrix(a,1,m+1,1,m+1); free_ivector(indx,1,m+1); Numerical Recipes -43108 As output,savgol returns the coefficients cn,for -nL<n nR.These are stored in c in"wrap-around order";that is,co is in c[1],c-1 is in c[2],and so on for further negative (outside indices.The value ci is stored in c [np],c2 in c [np-1],and so on for positive indices.This order may seem arcane,but it is the natural one where causal filters have nonzero coefficients North Software. in low array elements of c.It is also the order required by the function convlv in 813.1, which can be used to apply the digital filter to a data set. The accompanying table shows some typical output from savgol.For orders 2 and 4,the coefficients of Savitzky-Golay filters with several choices of n and nR are shown. The central column is the coefficient applied to the data fa in obtaining the smoothed ga Coefficients to the left are applied to earlier data;to the right,to later.The coefficients always add(within roundoff error)to unity.One sees that,as befits a smoothing operator, the coefficients always have a central positive lobe,but with smaller,outlying corrections of both positive and negative sign.In practice,the Savitzky-Golay filters are most useful for much larger values of nL and nR,since these few-point formulas can accomplish only a relatively small amount of smoothing. Figure 14.8.1 shows a numerical experiment using a 33 point smoothing filter,that is,652 Chapter 14. Statistical Description of Data Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include <math.h> #include "nrutil.h" void savgol(float c[], int np, int nl, int nr, int ld, int m) Returns in c[1..np], in wrap-around order (N.B.!) consistent with the argument respns in routine convlv, a set of Savitzky-Golay filter coefficients. nl is the number of leftward (past) data points used, while nr is the number of rightward (future) data points, making the total number of data points used nl + nr + 1. ld is the order of the derivative desired (e.g., ld = 0 for smoothed function). m is the order of the smoothing polynomial, also equal to the highest conserved moment; usual values are m = 2 or m = 4. { void lubksb(float **a, int n, int *indx, float b[]); void ludcmp(float **a, int n, int *indx, float *d); int imj,ipj,j,k,kk,mm,*indx; float d,fac,sum,**a,*b; if (np < nl+nr+1 || nl < 0 || nr < 0 || ld > m || nl+nr < m) nrerror("bad args in savgol"); indx=ivector(1,m+1); a=matrix(1,m+1,1,m+1); b=vector(1,m+1); for (ipj=0;ipj<=(m << 1);ipj++) { Set up the normal equations of the desired sum=(ipj ? 0.0 : 1.0); least-squares fit. for (k=1;k<=nr;k++) sum += pow((double)k,(double)ipj); for (k=1;k<=nl;k++) sum += pow((double)-k,(double)ipj); mm=IMIN(ipj,2*m-ipj); for (imj = -mm;imj<=mm;imj+=2) a[1+(ipj+imj)/2][1+(ipj-imj)/2]=sum; } ludcmp(a,m+1,indx,&d); Solve them: LU decomposition. for (j=1;j<=m+1;j++) b[j]=0.0; b[ld+1]=1.0; Right-hand side vector is unit vector, depending on which derivative we want. lubksb(a,m+1,indx,b); Get one row of the inverse matrix. for (kk=1;kk<=np;kk++) c[kk]=0.0; Zero the output array (it may be bigger than for (k = -nl;k<=nr;k++) { number of coefficients). sum=b[1]; Each Savitzky-Golay coefficient is the dot product of powers of an integer with the inverse matrix row. fac=1.0; for (mm=1;mm<=m;mm++) sum += b[mm+1]*(fac *= k); kk=((np-k) % np)+1; Store in wrap-around order. c[kk]=sum; } free_vector(b,1,m+1); free_matrix(a,1,m+1,1,m+1); free_ivector(indx,1,m+1); } As output, savgol returns the coefficients cn, for −nL ≤ n ≤ nR. These are stored in c in “wrap-around order”; that is, c0 is in c[1], c−1 is in c[2], and so on for further negative indices. The value c1 is stored in c[np], c2 in c[np-1], and so on for positive indices. This order may seem arcane, but it is the natural one where causal filters have nonzero coefficients in low array elements of c. It is also the order required by the function convlv in §13.1, which can be used to apply the digital filter to a data set. The accompanying table shows some typical output from savgol. For orders 2 and 4, the coefficients of Savitzky-Golay filters with several choices of nL and nR are shown. The central column is the coefficient applied to the data fi in obtaining the smoothed gi. Coefficients to the left are applied to earlier data; to the right, to later. The coefficients always add (within roundoff error) to unity. One sees that, as befits a smoothing operator, the coefficients always have a central positive lobe, but with smaller, outlying corrections of both positive and negative sign. In practice, the Savitzky-Golay filters are most useful for much larger values of nL and nR, since these few-point formulas can accomplish only a relatively small amount of smoothing. Figure 14.8.1 shows a numerical experiment using a 33 point smoothing filter, that is