650 Chapter 14.Statistical Description of Data 14.8 Savitzky-Golay Smoothing Filters In $13.5 we learned something about the construction and application of digital filters, but little guidance was given on which particular filter to use.That,of course,depends on what you want to accomplish by filtering.One obvious use for low-pass filters is to smooth noisy data. The premise of data smoothing is that one is measuring a variable that is both slowly varying and also corrupted by random noise.Then it can sometimes be useful to replace each data point by some kind of local average of surrounding data points.Since nearby points measure very nearly the same underlying value,averaging can reduce the level of noise without (much)biasing the value obtained. 8 We must comment editorially that the smoothing of data lies in a murky area,beyond the fringe of some better posed,and therefore more highly recommended,techniques that are discussed elsewhere in this book.If you are fitting data to a parametric model,for example (see Chapter 15),it is almost always better to use raw data than to use data that has been pre-processed by a smoothing procedure.Another alternative to blind smoothing is so-called optimal"or Wiener filtering,as discussed in $13.3 and more generally in $13.6.Data smoothing is probably most justified when it is used simply as a graphical technique,to guide the eye through a forest of data points all with large error bars;or as a means of making initial rough estimates of simple parameters from a graph. 令 In this section we discuss a particular type of low-pass filter,well-adapted for data smoothing,and termed variously Savitzky-Golay[1],least-squares [2],or DISPO (Digital Smoothing Polynomial)[3]filters.Rather than having their properties defined in the Fourier domain,and then translated to the time domain,Savitzky-Golay filters derive directly from 、 Press. a particular formulation of the data smoothing problem in the time domain,as we will now see.Savitzky-Golay filters were initially (and are still often)used to render visible the relative widths and heights of spectral lines in noisy spectrometric data. Recall that a digital filter is applied to a series of equally spaced data values fa=f(t:), whereti=to +i for some constant sample spacing A and i =...-2,-1,0,1,2,.... SCIENTIFIC We have seen (813.5)that the simplest type of digital filter (the nonrecursive or finite impulse response filter)replaces each data value fi by a linear combination gi of itself and some number of nearby neighbors. nR 9= Cnfitn (14.8.1) n=-nL Here n is the number of points used"to the left"of a data point i,i.e.,earlier than it,while nR is the number used to the right,i.e.,later.A so-called causal filter would have n =0. 10621 As a starting point for understanding Savitzky-Golay filters,consider the simplest possible averaging procedure:For some fixed nL=nR,compute each g:as the average of 43106 the data points from fto fThis is sometimes called moving window averaging Numerical Recipes and corresponds to equation (14.8.1)with constant c 1/(nL +nR+1).If the underlying function is constant,or is changing linearly with time (increasing or decreasing),then no (outside bias is introduced into the result.Higher points at one end of the averaging interval are on the average balanced by lower points at the other end.A bias is introduced,however,if 首 Software. the underlying function has a nonzero second derivative.At a local maximum,for example, moving window averaging always reduces the function value.In the spectrometric application, a narrow spectral line has its height reduced and its width increased.Since these parameters are themselves of physical interest,the bias introduced is distinctly undesirable. Note,however,that moving window averaging does preserve the area under a spectral line,which is its zeroth moment,and also (if the window is symmetric with n=nR)its mean position in time,which is its first moment.What is violated is the second moment. equivalent to the line width. The idea of Savitzky-Golay filtering is to find filter coefficients cn that preserve higher moments.Equivalently,the idea is 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 fi,we least-squares fit a polynomial to all
650 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 machinereadable 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). 14.8 Savitzky-Golay Smoothing Filters In §13.5 we learned something about the construction and application of digital filters, but little guidance was given on which particular filter to use. That, of course, depends on what you want to accomplish by filtering. One obvious use for low-pass filters is to smooth noisy data. The premise of data smoothing is that one is measuring a variable that is both slowly varying and also corrupted by random noise. Then it can sometimes be useful to replace each data point by some kind of local average of surrounding data points. Since nearby points measure very nearly the same underlying value, averaging can reduce the level of noise without (much) biasing the value obtained. We must comment editorially that the smoothing of data lies in a murky area, beyond the fringe of some better posed, and therefore more highly recommended, techniques that are discussed elsewhere in this book. If you are fitting data to a parametric model, for example (see Chapter 15), it is almost always better to use raw data than to use data that has been pre-processed by a smoothing procedure. Another alternative to blind smoothing is so-called “optimal” or Wiener filtering, as discussed in §13.3 and more generally in §13.6. Data smoothing is probably most justified when it is used simply as a graphical technique, to guide the eye through a forest of data points all with large error bars; or as a means of making initial rough estimates of simple parameters from a graph. In this section we discuss a particular type of low-pass filter, well-adapted for data smoothing, and termed variously Savitzky-Golay [1], least-squares [2], or DISPO (Digital Smoothing Polynomial) [3] filters. Rather than having their properties defined in the Fourier domain, and then translated to the time domain, Savitzky-Golay filters derive directly from a particular formulation of the data smoothing problem in the time domain, as we will now see. Savitzky-Golay filters were initially (and are still often) used to render visible the relative widths and heights of spectral lines in noisy spectrometric data. Recall that a digital filter is applied to a series of equally spaced data values fi ≡ f(ti), where ti ≡ t0 + i∆ for some constant sample spacing ∆ and i = ... − 2, −1, 0, 1, 2,... . We have seen (§13.5) that the simplest type of digital filter (the nonrecursive or finite impulse response filter) replaces each data value fi by a linear combination gi of itself and some number of nearby neighbors, gi = nR n=−nL cnfi+n (14.8.1) Here nL is the number of points used “to the left” of a data point i, i.e., earlier than it, while nR is the number used to the right, i.e., later. A so-called causal filter would have nR = 0. As a starting point for understanding Savitzky-Golay filters, consider the simplest possible averaging procedure: For some fixed nL = nR, compute each gi as the average of the data points from fi−nL to fi+nR . This is sometimes called moving window averaging and corresponds to equation (14.8.1) with constant cn = 1/(nL + nR + 1). If the underlying function is constant, or is changing linearly with time (increasing or decreasing), then no bias is introduced into the result. Higher points at one end of the averaging interval are on the average balanced by lower points at the other end. A bias is introduced, however, if the underlying function has a nonzero second derivative. At a local maximum, for example, moving window averaging always reduces the function value. In the spectrometric application, a narrow spectral line has its height reduced and its width increased. Since these parameters are themselves of physical interest, the bias introduced is distinctly undesirable. Note, however, that moving window averaging does preserve the area under a spectral line, which is its zeroth moment, and also (if the window is symmetric with nL = nR) its mean position in time, which is its first moment. What is violated is the second moment, equivalent to the line width. The idea of Savitzky-Golay filtering is to find filter coefficients cn that preserve higher moments. Equivalently, the idea is 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 fi, we least-squares fit a polynomial to all
14.8 Savitzky-Golay Smoothing Filters 651 nL nR Sample Savitzky-Golay Coefficients -0.0860.3430.4860.343-0.086 -0.1430.1710.3430.3710.257 2 0 0.086-0.143-0.0860.2570.886 2 -0.0840.0210.1030.1610.1960.2070.1960.1610.1030.021 -0.084 0.035-0.1280.0700.3150.4170.3150.070-0.1280.035 5 0.042-0.105-0.0230.1400.2800.3330.2800.140-0.023-0.1050.042 nL+nR+1 points in the moving window,and then set gi to be the value of that polynomial at position i.(If you are not familiar with least-squares fitting,you might want to look ahead to Chapter 15.)We make no use of the value of the polynomial at any other point.When we nted for move on to the next point f+,we do a whole new least-squares fit using a shifted window. All these least-squares fits would be laborious if done as described.Luckily,since the process of least-squares fitting involves only a linear matrix inversion,the coefficients of a fitted polynomial are themselves linear in the values of the data.That means that we can do all the fitting in advance,for fictitious data consisting of all zeros except for a single 1,and then do the fits on the real data just by taking linear combinations.This is the key point,then: 高A 令 There are particular sets of filter coefficients cn for which equation (14.8.1)"automatically" accomplishes the process of polynomial least-squares fitting inside a moving window. Toderive such cs consider how might obtained:We want tofit a polynomial of degree M in i,namely ao+ai+..+ai" to the values f.-nL,··,JnR Press. Then go will be the value of that polynomial at i 0,namely ao.The design matrix for this problem (815.4)is Programs 9 Aij=ii=-nL;...:nR:j=0,...:M (14.8.2) and the normal equations for the vector of a;'s in terms of the vector of fi's is in matrix notation (AT.A).a=AT.f or a=(AT.A)-1.(AT.f) (14.8.3) 6 We also have the specific forms AA= (14.8.4) k=一nL k三一卫 and {ar.,=定h=2 (14.8.5) Numerical -43126 Since the coefficient c is the component ao when f is replaced by the unit vector en, -n≤n1, the array c must be multiplied by k!to give derivative coefficients.For derivatives,one usually wants m 4 or larger
14.8 Savitzky-Golay Smoothing Filters 651 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 machinereadable 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). M nL nR Sample Savitzky-Golay Coefficients 2 2 2 −0.086 0.343 0.486 0.343 −0.086 2 3 1 −0.143 0.171 0.343 0.371 0.257 2 4 0 0.086 −0.143 −0.086 0.257 0.886 2 5 5 −0.084 0.021 0.103 0.161 0.196 0.207 0.196 0.161 0.103 0.021 −0.084 4 4 4 0.035 −0.128 0.070 0.315 0.417 0.315 0.070 −0.128 0.035 4 5 5 0.042 −0.105 −0.023 0.140 0.280 0.333 0.280 0.140 −0.023 −0.105 0.042 nL + nR + 1 points in the moving window, and then set gi to be the value of that polynomial at position i. (If you are not familiar with least-squares fitting, you might want to look ahead to Chapter 15.) We make no use of the value of the polynomial at any other point. When we move on to the next point fi+1, we do a whole new least-squares fit using a shifted window. All these least-squares fits would be laborious if done as described. Luckily, since the process of least-squares fitting involves only a linear matrix inversion, the coefficients of a fitted polynomial are themselves linear in the values of the data. That means that we can do all the fitting in advance, for fictitious data consisting of all zeros except for a single 1, and then do the fits on the real data just by taking linear combinations. This is the key point, then: There are particular sets of filter coefficients cn for which equation (14.8.1) “automatically” accomplishes the process of polynomial least-squares fitting inside a moving window. To derive such coefficients, consider how g0 might be obtained: We want to fit a polynomial of degree M in i, namely a0 + a1i + ··· + aMi M to the values f−nL ,...,fnR . Then g0 will be the value of that polynomial at i = 0, namely a0. The design matrix for this problem (§15.4) is Aij = i j i = −nL,...,nR, j = 0,...,M (14.8.2) and the normal equations for the vector of aj ’s in terms of the vector of fi’s is in matrix notation (AT · A) · a = AT · f or a = (AT · A) −1 · (AT · f) (14.8.3) We also have the specific forms AT · A ij = nR k=−nL AkiAkj = nR k=−nL ki+j (14.8.4) and AT · f j = nR k=−nL Akj fk = nR k=−nL kj fk (14.8.5) Since the coefficient cn is the component a0 when f is replaced by the unit vector en, −nL ≤ n 1, the array c must be multiplied by k! to give derivative coefficients. For derivatives, one usually wants m = 4 or larger.
652 Chapter 14.Statistical Description of Data #include #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 machinereadable 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 #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 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
14.8 Savitzky-Golay Smoothing Filters 653 8 6 before 4 0 0 100 200 300 400 500 600 700 800 900 6 after square (16,16,0) 4 83g nted for 19881992 0 Ennhh1 1-.200 0 100 200 300400500 600700800 900 to any 0T11+ server 6 after S-G(16.16.4) (Nort America computer, University Press. THE ART 0 100 200300400500600700 800 900 9 Progra Figure 14.8.1.Top:Synthetic noisy data consisting of a sequence of progressively narrower bumps, and additive Gaussian white noise.Center:Result of smoothing the data by a simple moving window average.The window extends 16 points leftward and rightward,for a total of 33 points.Note that narrow features are broadened and suffer corresponding loss of amplitude.The dotted curve is the underlying function used to generate the synthetic data.Bottom:Result of smoothing the data by a Savitzky-Golay to dir smoothing filter (of degree 4)using the same 33 points.While there is less smoothing of the broadest feature,narrower features have their heights and widths preserved. OF SCIENTIFIC COMPUTING (ISBN nL=nR =16.The upper panel shows a test function,constructed to have six"bumps"of varying widths,all of height 8 units.To this function Gaussian white noise of unit variance has been added.(The test function without noise is shown as the dotted curves in the center 10621 and lower panels.)The widths of the bumps(full width at half of maximum,or FWHM)are Numerica 140.43.24.17.13.and 10.respectively. 43106 The middle panel of Figure 14.8.1 shows the result of smoothing by a moving window average.One sees that the window of width 33 does quite a nice job of smoothing the broadest Recipes bump,but that the narrower bumps suffer considerable loss of height and increase of width The underlying signal (dotted)is very badly represented. The lower panel shows the result of smoothing with a Savitzky-Golay filter of the North Software. identical width,and degree M=4.One sees that the heights and widths of the bumps are quite extraordinarily preserved.A trade-off is that the broadest bump is less smoothed.That is because the central positive lobe of the Savitzky-Golay filter coefficients fills only a fraction of the full 33 point width.As a rough guideline,best results are obtained when the full width of the degree 4 Savitzky-Golay filter is between 1 and 2 times the FWHM of desired features in the data.(References [3]and [4]give additional practical hints.) Figure 14.8.2 shows the result of smoothing the same noisy "data"with broader Savitzky-Golay filters of 3 different orders.Here we have nL=nR=32(65 point filter) and M =2,4,6.One sees that,when the bumps are too narrow with respect to the filter size,then even the Savitzky-Golay filter must at some point give out.The higher order filter manages to track narrower features,but at the cost of less smoothing on broad features To summarize:Within limits,Savitzky-Golay filtering does manage to provide smoothing
14.8 Savitzky-Golay Smoothing Filters 653 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 machinereadable 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). 8 6 4 2 0 after square (16,16,0) 0 100 200 300 400 500 600 700 800 900 8 6 4 2 0 after S–G (16,16,4) 0 100 200 300 400 500 600 700 800 900 8 6 4 2 0 before 0 100 200 300 400 500 600 700 800 900 Figure 14.8.1. Top: Synthetic noisy data consisting of a sequence of progressively narrower bumps, and additive Gaussian white noise. Center: Result of smoothing the data by a simple moving window average. The window extends 16 points leftward and rightward, for a total of 33 points. Note that narrow features are broadened and suffer corresponding loss of amplitude. The dotted curve is the underlying function used to generate the synthetic data. Bottom: Result of smoothing the data by a Savitzky-Golay smoothing filter (of degree 4) using the same 33 points. While there is less smoothing of the broadest feature, narrower features have their heights and widths preserved. nL = nR = 16. The upper panel shows a test function, constructed to have six “bumps” of varying widths, all of height 8 units. To this function Gaussian white noise of unit variance has been added. (The test function without noise is shown as the dotted curves in the center and lower panels.) The widths of the bumps (full width at half of maximum, or FWHM) are 140, 43, 24, 17, 13, and 10, respectively. The middle panel of Figure 14.8.1 shows the result of smoothing by a moving window average. One sees that the window of width 33 does quite a nice job of smoothing the broadest bump, but that the narrower bumps suffer considerable loss of height and increase of width. The underlying signal (dotted) is very badly represented. The lower panel shows the result of smoothing with a Savitzky-Golay filter of the identical width, and degree M = 4. One sees that the heights and widths of the bumps are quite extraordinarily preserved. A trade-off is that the broadest bump is less smoothed. That is because the central positive lobe of the Savitzky-Golay filter coefficients fills only a fraction of the full 33 point width. As a rough guideline, best results are obtained when the full width of the degree 4 Savitzky-Golay filter is between 1 and 2 times the FWHM of desired features in the data. (References [3] and [4] give additional practical hints.) Figure 14.8.2 shows the result of smoothing the same noisy “data” with broader Savitzky-Golay filters of 3 different orders. Here we have nL = nR = 32 (65 point filter) and M = 2, 4, 6. One sees that, when the bumps are too narrow with respect to the filter size, then even the Savitzky-Golay filter must at some point give out. The higher order filter manages to track narrower features, but at the cost of less smoothing on broad features. To summarize: Within limits, Savitzky-Golay filtering does manage to provide smoothing
654 Chapter 14.Statistical Description of Data 6 after S-G(32.32.2) 2 0 0 100 200 300 400500600 700800 900 after S-G(32.32,4) E 83g 23 granted for 19881992 0 1-.200 0 100 200 300 400500 600 700 800 900 from NUMERICAL RECIPESI 8TT+T++T 6 after S-G(32.32,6) server (Nor to make 4 America computer, THE mw,,1,,,,1,,,,1, ART 0 100.200300400500600700800 900 9 Programs Figure 14.8.2.Result of applying wider 65 point Savitzky-Golay filters to the same data set as in Figure 14.8.1.Top:degree 2.Center:degree 4.Bottom:degree 6.All of these filters are inoptimally broad for the resolution of the narrow features.Higher-order filters do best at preserving feature heights and widths,but do less smoothing on broader features. 兰云侯 to dir OF SCIENTIFIC COMPUTING(ISBN without loss of resolution.It does this by assuming that relatively distant data points have 1988-19920 some significant redundancy that can be used to reduce the level of noise.The specific nature of the assumed redundancy is that the underlying function should be locally well-fitted by a 10-621 polynomial.When this is true,as it is for smooth line profiles not too much narrower than the filter width,then the performance of Savitzky-Golay filters can be spectacular.When it idge.org Numerical Recipes 43108 is not true,then these filters have no compelling advantage over other classes of smoothing filter coefficients. A last remark concerns irregularly sampled data,where the values fi are not uniformly (outside spaced in time.The obvious generalization of Savitzky-Golay filtering would be to do a least-squares fit within a moving window around each data point,one containing a fixed North Software. number of data points to the left (nL)and right(nR).Because of the irregular spacing. however,there is no way to obtain universal filter coefficients applicable to more than one data point.One must instead do the actual least-squares fits for each data point.This becomes computationally burdensome for larger nL,nR,and M. As a cheap alternative,one can simply pretend that the data points are equally spaced. This amounts to virtually shifting,within each moving window,the data points to equally spaced positions.Such a shift introduces the equivalent of an additional source of noise into the function values.In those cases where smoothing is useful,this noise will often be much smaller than the noise already present.Specifically,if the location of the points is approximately random within the window,then a rough criterion is this:If the change in f across the full width of the N =nL+nR+1 point window is less than VN/2 times the measurement noise on a single point,then the cheap method can be used
654 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 machinereadable 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). after S–G (32,32,4) after S–G (32,32,2) 8 6 4 2 0 0 100 200 300 400 500 600 700 800 900 8 6 4 2 0 after S–G (32,32,6) 0 100 200 300 400 500 600 700 800 900 8 6 4 2 0 0 100 200 300 400 500 600 700 800 900 Figure 14.8.2. Result of applying wider 65 point Savitzky-Golay filters to the same data set as in Figure 14.8.1. Top: degree 2. Center: degree 4. Bottom: degree 6. All of these filters are inoptimally broad for the resolution of the narrow features. Higher-order filters do best at preserving feature heights and widths, but do less smoothing on broader features. without loss of resolution. It does this by assuming that relatively distant data points have some significant redundancy that can be used to reduce the level of noise. The specific nature of the assumed redundancy is that the underlying function should be locally well-fitted by a polynomial. When this is true, as it is for smooth line profiles not too much narrower than the filter width, then the performance of Savitzky-Golay filters can be spectacular. When it is not true, then these filters have no compelling advantage over other classes of smoothing filter coefficients. A last remark concerns irregularly sampled data, where the values fi are not uniformly spaced in time. The obvious generalization of Savitzky-Golay filtering would be to do a least-squares fit within a moving window around each data point, one containing a fixed number of data points to the left (nL) and right (nR). Because of the irregular spacing, however, there is no way to obtain universal filter coefficients applicable to more than one data point. One must instead do the actual least-squares fits for each data point. This becomes computationally burdensome for larger nL, nR, and M. As a cheap alternative, one can simply pretend that the data points are equally spaced. This amounts to virtually shifting, within each moving window, the data points to equally spaced positions. Such a shift introduces the equivalent of an additional source of noise into the function values. In those cases where smoothing is useful, this noise will often be much smaller than the noise already present. Specifically, if the location of the points is approximately random within the window, then a rough criterion is this: If the change in f across the full width of the N = nL + nR + 1 point window is less than N/2 times the measurement noise on a single point, then the cheap method can be used
14.8 Savitzky-Golay Smoothing Filters 655 CITED REFERENCES AND FURTHER READING: Savitzky A.,and Golay,M.J.E.1964,Analytica/Chemistry,vol.36,pp.1627-1639.[1] Hamming,R.W.1983,Digita/Filters,2nd ed.(Englewood Cliffs,NJ:Prentice-Hall).[2] Ziegler,H.1981,Applied Spectroscopy,vol.35,pp.88-92.[3] Bromba,M.U.A.,and Ziegler,H.1981,Analytica/Chemistry,vol.53,pp.1583-1586.[4] Permission is http://w.nr.com or call 1-800-872-7423(North America only),orsend email to directcustserv@cambridge.org (outside North America). granted for internet users to make one paper copy for their 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)
14.8 Savitzky-Golay Smoothing Filters 655 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 machinereadable 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). CITED REFERENCES AND FURTHER READING: Savitzky A., and Golay, M.J.E. 1964, Analytical Chemistry, vol. 36, pp. 1627–1639. [1] Hamming, R.W. 1983, Digital Filters, 2nd ed. (Englewood Cliffs, NJ: Prentice-Hall). [2] Ziegler, H. 1981, Applied Spectroscopy, vol. 35, pp. 88–92. [3] Bromba, M.U.A., and Ziegler, H. 1981, Analytical Chemistry, vol. 53, pp. 1583–1586. [4]