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 all650 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). 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