c14-8 (779583), страница 2
Текст из файла (страница 2)
These are stored inc in “wrap-around order”; that is, c0 is in c[1], c−1 is in c[2], and so on for further negativeindices. The value c1 is stored in c[np], c2 in c[np-1], and so on for positive indices. Thisorder may seem arcane, but it is the natural one where causal filters have nonzero coefficientsin 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 and4, 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 coefficientsalways 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 correctionsof both positive and negative sign. In practice, the Savitzky-Golay filters are most usefulfor much larger values of nL and nR , since these few-point formulas can accomplish onlya relatively small amount of smoothing.Figure 14.8.1 shows a numerical experiment using a 33 point smoothing filter, that is,Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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 inroutine 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 totalnumber of data points used nl + nr + 1. ld is the order of the derivative desired (e.g., ld = 0for smoothed function). m is the order of the smoothing polynomial, also equal to the highestconserved 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;65314.8 Savitzky-Golay Smoothing Filters86420before86420100200300400500600700800900300400500600700800900300400500600700800900after square (16,16,0)086420100200after S – G (16,16,4)0100200Figure 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 windowaverage. The window extends 16 points leftward and rightward, for a total of 33 points. 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-Golaysmoothing filter (of degree 4) using the same 33 points. While there is less smoothing of the broadestfeature, narrower features have their heights and widths preserved.nL = nR = 16.
The upper panel shows a test function, constructed to have six “bumps” ofvarying 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 windowaverage. One sees that the window of width 33 does quite a nice job of smoothing the broadestbump, 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 theidentical width, and degree M = 4. One sees that the heights and widths of the bumps arequite extraordinarily preserved.
A trade-off is that the broadest bump is less smoothed. Thatis because the central positive lobe of the Savitzky-Golay filter coefficients fills only a fractionof the full 33 point width. As a rough guideline, best results are obtained when the full widthof the degree 4 Savitzky-Golay filter is between 1 and 2 times the FWHM of desired featuresin the data. (References [3] and [4] give additional practical hints.)Figure 14.8.2 shows the result of smoothing the same noisy “data” with broaderSavitzky-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 filtersize, then even the Savitzky-Golay filter must at some point give out. The higher order filtermanages 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 smoothingSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).065486420Chapter 14.Statistical Description of Dataafter S – G (32,32,2)86420100200300400500600700800900300400500600700800900300400500600700800900after S – G (32,32,4)086420100200after S – G (32,32,6)0100200Figure 14.8.2. Result of applying wider 65 point Savitzky-Golay filters to the same data set as in Figure14.8.1. Top: degree 2.
Center: degree 4. Bottom: degree 6. All of these filters are inoptimally broadfor the resolution of the narrow features. Higher-order filters do best at preserving feature heights andwidths, but do less smoothing on broader features.without loss of resolution. It does this by assuming that relatively distant data points havesome significant redundancy that can be used to reduce the level of noise. The specific natureof the assumed redundancy is that the underlying function should be locally well-fitted by apolynomial. When this is true, as it is for smooth line profiles not too much narrower thanthe filter width, then the performance of Savitzky-Golay filters can be spectacular. When itis not true, then these filters have no compelling advantage over other classes of smoothingfilter coefficients.A last remark concerns irregularly sampled data, where the values fi are not uniformlyspaced in time.
The obvious generalization of Savitzky-Golay filtering would be to do aleast-squares fit within a moving window around each data point, one containing a fixednumber 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 onedata point.
One must instead do the actual least-squares fits for each data point. This becomescomputationally 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 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. Specifically, if the location of the points isapproximately random within the window, then a rough criterion is this: Ifpthe change in facross the full width of the N = nL + nR + 1 point window is less than N/2 times themeasurement noise on a single point, then the cheap method can be used.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
