c14-8 (779583)
Текст из файла
650Chapter 14.Statistical Description of Data14.8 Savitzky-Golay Smoothing Filtersgi =nRXcnfi+n(14.8.1)n=−nLHere nL is the number of points used “to the left” of a data point i, i.e., earlier than it, whilenR 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 simplestpossible averaging procedure: For some fixed nL = nR , compute each gi as the average ofthe data points from fi−nL to fi+nR . This is sometimes called moving window averagingand corresponds to equation (14.8.1) with constant cn = 1/(nL + nR + 1).
If the underlyingfunction 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 onthe average balanced by lower points at the other end. A bias is introduced, however, ifthe 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 parametersare themselves of physical interest, the bias introduced is distinctly undesirable.Note, however, that moving window averaging does preserve the area under a spectralline, which is its zeroth moment, and also (if the window is symmetric with nL = nR ) itsmean 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 highermoments. 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, typically quadratic or quartic: For each point fi , we least-squares fit a polynomial to allSample 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).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, dependson what you want to accomplish by filtering. One obvious use for low-pass filters is tosmooth noisy data.The premise of data smoothing is that one is measuring a variable that is both slowlyvarying 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 nearbypoints measure very nearly the same underlying value, averaging can reduce the level of noisewithout (much) biasing the value obtained.We must comment editorially that the smoothing of data lies in a murky area, beyondthe fringe of some better posed, and therefore more highly recommended, techniques thatare discussed elsewhere in this book.
If you are fitting data to a parametric model, forexample (see Chapter 15), it is almost always better to use raw data than to use data thathas been pre-processed by a smoothing procedure. Another alternative to blind smoothing isso-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, toguide the eye through a forest of data points all with large error bars; or as a means of makinginitial rough estimates of simple parameters from a graph.In this section we discuss a particular type of low-pass filter, well-adapted for datasmoothing, and termed variously Savitzky-Golay [1], least-squares [2], or DISPO (DigitalSmoothing Polynomial) [3] filters.
Rather than having their properties defined in the Fourierdomain, 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-Golay filters were initially (and are still often) used to render visible the relativewidths 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 impulseresponse filter) replaces each data value fi by a linear combination gi of itself and somenumber of nearby neighbors,65114.8 Savitzky-Golay Smoothing FiltersM nL nRSample Savitzky-Golay Coefficients−0.086 0.343 0.486 0.343 −0.0862222312402554440.035 −0.128 0.070 0.315 0.417 0.315 0.070 −0.128 0.0354550.042 −0.105 −0.023 0.140 0.280 0.333 0.280 0.140 −0.023 −0.105 0.042−0.143 0.171 0.343 0.371 0.2570.086 −0.143 −0.086 0.257 0.8860.1030.161 0.196 0.207 0.196 0.1610.1030.021 −0.084nL + nR + 1 points in the moving window, and then set gi to be the value of that polynomialat position i. (If you are not familiar with least-squares fitting, you might want to look aheadto Chapter 15.) We make no use of the value of the polynomial at any other point.
When wemove 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 theprocess 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 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 apolynomial of degree M in i, namely a0 + a1 i + · · · + aM iM to the values f−nL , . . . , fnR .Then g0 will be the value of that polynomial at i = 0, namely a0 .
The design matrix forthis problem (§15.4) isi = −nL , . . . , nR ,Aij = ijj = 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 · fora = (AT · A)−1 · (AT · f)We also have the specific formsnRnRnoXXAT · A=Aki Akj =ki+jijk=−nL(14.8.3)(14.8.4)k=−nLandnAT · fo=jnRXAkj fk =k=−nLnRXk j fk(14.8.5)k=−nLSince the coefficient cn is the component a0 when f is replaced by the unit vector en ,−nL ≤ n < nR , we haveM nonoX(AT · A)−1cn = (AT · A)−1 · (AT · en ) =0m=0nm(14.8.6)0mNote that equation (14.8.6) says that we need only one row of the inverse matrix. (Numericallywe can get this by LU decomposition with only a single backsubstitution.)The function savgol, below, implements equation (14.8.6). As input, it takes theparameters nl = nL , nr = nR , and m = M (the desired order).
Also input is np, thephysical length of the output array c, and a parameter ld which for data fitting should bezero. In fact, ld specifies which coefficient among the ai ’s should be returned, and we arehere 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, thefiltered first derivative is the convolution (14.8.1) divided by the stepsize ∆.
For derivatives,one usually wants m = 4 or larger.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).−0.084 0.021652Chapter 14.Statistical Description of Data#include <math.h>#include "nrutil.h"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 desiredsum=(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 thanfor (k = -nl;k<=nr;k++) {number of coefficients).sum=b[1];Each Savitzky-Golay coefficient is the dotfac=1.0;product of powers of an integer with thefor (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;}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 .
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
