f14-8

644 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).CITED REFERENCES AND FURTHER READING: Fasano, G. and Franceschini, A. 1987, Monthly Notices of the Royal Astronomical Society , vol. 225, pp. 155–170. [1] Peacock,J.A.1983, MonthlyNoticesoftheRoyalAstronomicalSociety ,vol.202,pp.615–627.[2] Spergel, D.N., Piran, T., Loeb, A., Goodman, J., and Bahcall, J.N. 1987, Science, vol. 237, pp. 1471–1473. [3] 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 replaceeach data point by some kind of local average of surrounding data points. Since nearbypointsmeasurevery nearlythesameunderlying value, averagingcanreduce thelevelofnoisewithout (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 arediscussed 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 beenpre-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 isprobably most justifiedwhen itisused simply asa graphical technique, to guidetheeyethrough aforestofdatapointsallwithlargeerrorbars;orasameans ofmaking initialroughestimates 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],o rDISPO(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 froma particular formulation of the data smoothing problem in the time domain, as we will nowsee. Savitzky-Golayfilterswereinitially(andarestilloften)usedtorendervisibletherelativewidths and heights of spectral lines in noisy spectrometric data. Recall thata digital filteris applied to a series of equally spaced data values f i≡f(ti), where ti≡t0+i∆for some constant sample spacing ∆and i=...−2,−1,0,1,2,.... Wehave seen ( §13.5) that the simplesttype of digital filter(thenonrecursive or finiteimpulse response filter) replaces each data value fiby a linear combination giof itself and some number of nearby neighbors, gi=nR/summationdisplay n=−nLcnfi+n (14.8.1 ) Here nLis the number of points used “to the left” of a data point i, i.e., earlier than it, while nRis the number used to the right, i.e.,later. A so-called causalfilter 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 gias the average of the data points from fi−nLtofi+nR. This is sometimes called moving window averaging andcorresponds toequation(14.8.1)withconstant cn=1 /(nL+nR+1 ). Iftheunderlying function is constant, or is changing linearly with time (increasing or decreasing), then nobias is introduced into the result. Higher points at one end of the averaging interval are on 14.8Savitzky-GolaySmoothingFilters 645Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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,movingwindowaveragingalwaysreducesthefunctionvalue. Inthespectrometricapplication,a narrow spectral line has its height reduced and its width increased. Since these parametersare 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 L=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 c nthat preserve higher moments. Equivalently, the idea is to approximate the underlying function within the movingwindow not by a constant (whose estimate is the average), but by a polynomial of higherorder,typicallyquadraticorquartic: Foreachpoint f i,weleast-squares fitapolynomialtoall nL+nR+1points inthe moving window, and then set gitobe thevalue ofthatpolynomial 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 fitusing 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 afitted polynomial are themselves linear in the values of the data. That means that we can doall the fitting in advance, for fictitious data consisting of all zeros except for a single 1, andthen do the fitson the realdata justby taking linearcombinations. Thisisthe key point, then:There are particular sets of filter coefficients c nfor which equation (14.8.1) “automatically” accomplishes the process of polynomial least-squares fitting inside a moving window. To derive such coefficients, consider how g0might be obtained: We want to fit a polynomial of degree Mini, namely a0+a1i+···+aMiMto the values f−nL,...,f nR. Then g0will be the value of that polynomial at i=0, namely a0. The design matrix for this problem ( §15.4) is Aij=iji=−nL,...,n R,j =0 ,...,M (14.8.2 ) andthenormalequationsforthevectorof aj’sintermsofthevectorof fi’sisinmatrixnotation (AT·A)·a=AT·fora=(AT·A)−1·(AT·f)( 14.8.3 ) We also have the specific forms /braceleftBig AT·A/bracerightBig ij=nR/summationdisplay k=−nLAkiAkj=nR/summationdisplay k=−nLki+j(14.8.4 ) and /braceleftBig AT·f/bracerightBig j=nR/summationdisplay k=−nLAkjfk=nR/summationdisplay k=−nLkjfk (14.8.5 ) Since the coefficient cnis the component a0whenfis replaced by the unit vector en, −nL≤n<n R, we have cn=/braceleftBig (AT·A)−1·(AT·en)/bracerightBig 0=M/summationdisplay m=0/braceleftBig (AT·A)−1/bracerightBig 0mnm(14.8.6 ) Notethatequation(14.8.6)saysthatweneedonlyonerowoftheinversematrix. (Numerically we can get this by LUdecomposition with only a single backsubstitution.) The subroutine savgol, below, implements equation (14.8.6). As input, it takes the parameters nl=nL,nr=nR, and m=M(the desired order). Also input is np, the physical length of the output array c, and a parameter ldwhich for data fitting should be zero. In fact, ldspecifies which coefficient among the ai’s should be returned, and we are here interested in a0. For another purpose, namely the computation of numerical derivatives (already mentioned in §5.7) the useful choice is ld≥1. With ld=1, for example, the filteredfirstderivative isthe convolution (14.8.1)divided by the stepsize ∆.F o r ld=k>1, the array cmust be multiplied by k!to give derivative coefficients. For derivatives, one usually wants m=4or larger. 646 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).M nLnR Sample Savitzky-Golay Coefficients 222 −0.086 0.343 0.4860.343−0.086 231 −0.143 0.171 0.343 0.3710.257 240 0.086 −0.143−0.086 0.257 0.886 255−0.084 0.021 0.103 0.161 0.196 0.2070.196 0.161 0.103 0.021 −0.084 444 0.035 −0.128 0.070 0.315 0.4170.315 0.070 −0.128 0.035 4550.042−0.105−0.023 0.140 0.280 0.3330.280 0.140 −0.023−0.105 0.042 SUBROUTINE savgol(c,np,nl,nr,ld,m) INTEGER ld,m,nl,np,nr,MMAXREAL c(np) PARAMETER (MMAX=6) C USES lubksb,ludcmp Returns in c(1:np), in wrap-around order (N.B.!) consistent with the argument respns inroutine convlv, aset ofSavitzky-Golayfilter coefficients. nlis the number of leftward (past) data points used, while nris the number of rightward (future) data points, making the total number ofdata pointsused nl+nr+1.ldistheorderofthederivativedesired (e.g., ld =0for smoothed function). mis the order of the smoothing polynomial, also equal to the highest conserved moment; usual values are m=2orm=4. INTEGER imj,ipj,j,k,kk,mm,indx(MMAX+1)REAL d,fac,sum,a(MMAX+1,MMAX+1),b(MMAX+1) if(np.lt.nl+nr+1.or.nl.lt.0.or.nr.lt.0.or.ld.gt.m.or.m.gt.MMAX * .or.nl+nr.lt.m) pause ’bad args in savgol’ do 14ipj=0,2*m Set up the normal equations of the desired least- squares fit. sum=0. if(ipj.eq.0)sum=1. do11k=1,nr sum=sum+float(k)**ipj enddo 11 do12k=1,nl sum=sum+float(-k)**ipj enddo 12 mm=min(ipj,2*m-ipj)do 13imj=-mm,mm,2 a(1+(ipj+imj)/2,1+(ipj-imj)/2)=sum enddo 13 enddo 14 call ludcmp(a,m+1,MMAX+1,indx,d) Solve them: LUdecomposition. do15j=1,m+1 b(j)=0. enddo 15 b(ld+1)=1. Right-handsidevectorisunitvector,dependingonwhichderivativewewant. call lubksb(a,m+1,MMAX+1,indx,b) Backsubstitute,givingonerowoftheinversematrix. do16kk=1,np Zerotheoutputarray(itmaybebiggerthannumber of coefficients). c(kk)=0. enddo 16 do18k=-nl,nr Each Savitzky-Golay coefficient is the dot product of powers of an integer with the inverse matrixrow.sum=b(1) fac=1.do 17mm=1,m fac=fac*k sum=sum+b(mm+1)*fac enddo 17 kk=mod(np-k,np)+1 Store in wrap-around order. c(kk)=sum enddo 18 return END 14.8Savitzky-GolaySmoothingFilters 647Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).8 6420after square (16,16,0) 0 100 200 300 400 500 600 700 800 900 8 6420after S–G (16,16,4) 0 100 200 300 400 500 600 700 800 9008 6420before 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. Thewindow extends 16points leftward and rightward, for a total of 33points. Note that narrowfeatures are broadened and suffer corresponding loss of amplitude. The dotted curve is the underlyingfunction 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. As output, savgolreturns the coef ficients cn, for−nL≤n≤nR. These are stored in cin“wrap-around order ”; thatis, c0isinc(1),c−1isinc(2),and soon forfurther negative indices. The value c1is stored in c(np),c2inc(np-1), and so on for positive indices. This order may seem arcane, but itis the natural one where causal filtershave nonzero coef ficients in low array elements of c. It is also the order required by the subroutine convlvin§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 coef ficients of Savitzky-Golay filters with several choices of nLand nRare shown. The central column is the coef ficient applied to the data fiin obtaining the smoothed gi. Coefficients to the left are applied to earlier data; to the right, to later. The coef ficients always add (within roundoff error) to unity. One sees that, as be fits 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 nLand 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, nL=nR=1 6. 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 variancehas been added. (The test function without noise is shown as the dotted curves in the centerand lower panels.) The widths of the bumps (full width at half of maximum, or FWHM) are140, 43, 24, 17, 13, and 10, respectively. The middle panel of Figure 14.8.1 shows the result of smoothing by a moving window 648 Chapter14. StatisticalDescriptionofDataSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).after S–G (32,32,4)after S–G (32,32,2) 8 6420 0 100 200 300 400 500 600 700 800 900 86420after S–G (32,32,6) 0 100 200 300 400 500 600 700 800 9008 6420 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. average. Oneseesthatthewindowofwidth33doesquiteanicejobofsmoothingthebroadest 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. ThatisbecausethecentralpositivelobeoftheSavitzky-Golay filtercoefficientsfillsonlyafraction of the full 33 point width. As a rough guideline, best results are obtained when the full widthof the degree 4 Savitzky-Golay filteris 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=3 2(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. Tosummarize: Withinlimits,Savitzky-Golay filteringdoesmanagetoprovidesmoothing without loss of resolution. It does this by assuming that relatively distant data points havesome signi ficant redundancy thatcan be used to reduce the levelof noise. The speci fic 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 pro files not too much narrower than thefilter 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 f iare not uniformly spaced in time. The obvious generalization of Savitzky-Golay filtering would be to do a 14.8Savitzky-GolaySmoothingFilters 649Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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 datapoint. Onemustinsteaddotheactualleast-squares fitsforeachdatapoint. Thisbecomes computationally burdensome for larger nL,nR, and M. As a cheap alternative, one can simply pretend that the data points areequally spaced. This amounts to virtually shifting, within each moving window, the data points to equallyspaced positions. Such a shift introduces the equivalent of an additional source of noiseinto the function values. In those cases where smoothing is useful, this noise will often bemuch smaller than the noise already present. Speci fically, 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=n L+nR+1point window is less than/radicalbig N/2times the measurement noise on a single point, then the cheap method can be used. 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]