Anal. Chem. 1995, 67, 2758-2762
Smoothing and Differentiation by an Adaptive-Degree Polynomial Filter Phillip Barak Department of Soil Science, University of Wisconsin- Madison, Madison, Wisconsin 53706-7299
The Savitzky-Golay method for data smoothing and differentiation calculates convolutionweights using Gram polynomials that exactly reproduce the results of leastsquares polynomial regression. Use of the SavitzkyGolay method requires specification of both filter length and polynomial degree to calculate convolution weights. For maximum smoothing of statistical noise in data, polynomials with low degrees are desirable, while high polynomial degree is necessary for accurate reproduction of peaks in the data. Extension of the least-squares regression formalism with statistical testing of additional terms of polynomial degree to a heuristically chosen minimum for each data window leads to an adaptivedegree polynomial filter (ADPF). Based on noise reduction for data that consist of pure noise and on signal reproduction for data that is purely signal, ADPF performed nearly as well as the optimally chosen fixed-degree Savitzky-Golay filter and outperformed suboptimally chosen Savitzky-Golay filters. For synthetic data consisting of noise and signal, ADPF outperformed both optimally chosen and suboptimally chosen fixed-degree Savitzky-Golay filters.
original tabular data5 have been corrected7 and generalized to calculate the convolution weights at all positions, for all polynomial degrees, all filter lengths, and for any desired smoothed derivative,Q and the general treatment has been extended to unequally spaced datag and two dimensions.lO,ll The performance of polynomial filters using the SavitzkyGolay method in terms of noise reduction and signal enhancement has been determined numerically12 and calculated analytic all^.^^ Both filter length and polynomial degree together affect the smoothed values and the degree of noise reduction. Results of the smoothing operation are therefore dependent on a priori choices, and often several Savitzky-Golay smooths on the same data set using the various combinations of polynomial degree and filter lengths will be tested. Not only is it difficult to decide which combination of polynomial degree and filter length yields the optimum data smooth, but different combinations may be superior at various portions of the data set. If the true signal function were known, an objective figure, or figure of merit, could be constructed; in lieu of this information, use of the Savitzky-Golay method must rely on largely arbitrary decisions regarding smoothing parameters.
THEORY Data smoothing is an important part of the collection and interpretation of large sets of experimental data in which the signal is accompanied by noise. Although a number of techniques for data smoothing are known (e.g., Gaussian filtering,’ splines,2and Fourier transforms3), methods employing sliding, locally fit polynomial approximations have long been of interest4 A single publication by Savikky and Golay5 describing the use of piecewise polynomial fits to uniformly spaced data has been cited over 2000 times in the scientific literature since 1964; the technique is currently in wide use among chemists! both as a stand-alone routine and as part of proprietary software for data manipulation. The Savitzky-Golay method calculates a smoothed value of a polynomial function or its derivatives at a given point based on a series of convolution weights, or weighting factors, applied to neighboring points within a sliding data window. Convolution weights are derived from Gram polynomials, which exactly reproduce an equivalent least-squares polynomial regression. The E-mail:
[email protected]. (1) Marchand, P.; Marmet, L. Reo. Sci. instrum. 1983,54, 1034-1041. (2) Greville, T. N. E., Ed. Theory and Afifilications ofspline Functions; Academic Press: New York, 1969. (3) Press, W. H.; Flannery. B. P.; Teukolsky, S. A; Vetterling, W. T. Numerical Recipes; Cambridge University Press: Cambridge, U.K, 1986; Chapter 12. (4) Macaulay, F. R. The Smoothing of Time Series; National Bureau of Economic Research, Inc.: New York, 1931. (5) Savitzky, A; Golay, M. J. E. Anal. Chem. 1964,36, 1627-1639. (6) Gorry. P. A. Anal. Chem. 1990,62,570-573.
2758
Analytical Chemistry, Vol. 67, No. 77, September 7, 1995
The standard approach to Savitzky-Golay polynomial filters calculates the value of the smoothed function (s = 0) or desired derivatives of the function (s > 0) by the formula m
where h:,”.’is the convolution weight of the ith point to evaluate the sth derivative at point t using a polynomial of degree n on 2m 1 data points, y. The convolution weights may be calculated as6.8
+
where (a)(*)is a generalized factorial function (a) (a - 1) ... (a 6 + l),and (a)(o)= 0, and the Gram polynomials are defined as (7) Steiner, J.; Termonia. Y.; Deltour. J. Anal. Chem. 1972,44, 1906-1909. (8) Emst. R. R Adu. Magn. Reson. 1966,2, 1-135. (9) Gorry, P. A Anal. Chem. 1991,63, 534-536. (10) Ratzlaff, K L.; Johnson, J. T. Anal. Chem. 1989,61. 1303-1305. (11) Kuo, J. E.; Wang, H.: Pickup, S. Anal. Chem. 1991,63,630-645. (12) Enke, C. G; Nieman, T. A Anal. Chem. 1976,48, 705A-712A. (13) Phillips, G. R.; Harris, J. M. Anal. Chem. 1990,62,2749-2752.
0003-2700/95/0367-2758$9.00/0 0 1995 American Chemical Society
+
in the current data window and that j 3 polynomial degrees may be required in the next data window. The polynomial fit on the successive data window is therefore required to test additional terms from k = 0 to at least k =j 3 and to fail to add the highest two consecutive additional terms (up to either a user-defined maximum or N - 2) before concluding the search for significant polynomial degrees. For the initial data window, for which no value of j from a previous data window is available, j may default to 0 or be set to a user-defbed value. During testing for significance of additional polynomial degrees, the sum of squares of the residuals may be calculated using Gram polynomials based on previous values as
+
Filtering data by the Savitzky-Golay algorithm requires selection of a filter length 2m 1 and a polynomial degree n, determination of the corresponding convolution weights to evaluate the smoothed function or its derivatives at the central point of the filter length (t = 0), and application of the Savitzky-Golay filter in piecewise or sliding fashion for the length of the collected data. Interestingly, the Savitzky-Golay filter yields equivalent results for central point smoothing using even polynomial degree and the next highest odd degree, e.g., 0/1,2/3, etc., even though the convolution weights differ and the fit of the polynomial to the 2m 1 data points differs. Unlike the ked-degree Savitzky-Golay polynomial filter, standard polynomial regression analysis is based on changes in the sum of squares of residuals, x', upon changing the degree of the fitting polynomial,'* where
+
+
m
m
m
The sum of squares of residuals of two polynomials of degrees may be used to calculate a test statistic that follows the F distribution with degrees of freedom V I = n2 - n1 and vz = N - nz - 1, where N = 2m 1, as follows:15 HI and n2 (n2 > nl)
m
2 Xn
2 = Xn-1
m
-
The use of statistical testing of the significance of polynomial degree applied to each data window, together with a heuristic approach to the extent of additional terms tested, leads to a digital filter that varies the degree of the fitting polynomial as it slides down the data window. Such a filter may be termed an adaptivedegree polynomial filter (ADPF) and may be regarded as a modified or extended Savitzky-Golay filter, for which the polynomial degree is set a priori.
EXPERIMENTAL SECTION Noise Reduction. A set of lo00 normal deviates with a mean of zero and variance of unity were generated.17 For filter lengths of 5-23 data points and polynomial degrees from 0 to 9, convolution weights for the Savitzky-Golay method were calculated using eqs 1-3, and central values (t = 0) were calculated. Signal Fidelity. A Gaussian curve was generated of the form y(x) = expi- (z - p)2/(2uz)1, for x at unit intervals between -22 The resulting Fvalue is tested against critical values of F(V~,YZ), and $22, p = 0, and u = 15. Convolution weights for filter lengths available in standard statistical tables, at a preset probability level from 5 to 23 and polynomial degree from 0 to 9 were calculated, a, usually either 5 or 1%. The choice of polynomial degree may and the smoothed peak heights at x = t = 0 were calculated. be regarded as a sequence of hypothesis tests, with upper and Noise Reduction and S i a l Fidelity. A signal function was lower bounds for polynomial degrees to be tested. For standard, constructed of the form y(x) = C exp[-(x - pJ2/(2~?)I,x = 0, nonpiecewise polynomial regressions, the sequence of significance 0.05, 0.10, ..., where pi (I= 1-5) are 1, 1.62, 2.61, 4.22, and 6.82 tests commonly starts from the lowest degree and continues until and ui = ,dlO, appropriate to a constant chromatographic either one or two additional polynomial degrees fail the signifiseparation resolution of 2. Five sets of signal plus noise were cance test.149'6 created by adding normal deviates with a mean of zero and a For data sets with multiple peaks (e.g., chromatograms, IR standard deviation of 0.1 to the signal function. Overall fit and spectra, and X-ray diffractograms), it may be reasoned that each noise reduction is represented by the root mean square of the data window may have as many as one more signiicant peak than smoothed value minus the signal function without noise. the previous data window. Rather than fitting each piecewise Adaptive-Degree Polynomial Filter. For a given filter polynomial regression in a completely independent manner, length, convolution weights were calculated by using eqs 2 and 3 information from one data window can be used to heuristically fork = 0 to min{N - 2, ll}. For each data window, and for each determine the minimum polynomial degree to test for significance polynomial degree tested, F values were calculated using eqs 5 in the next data window. This information may be calculated from and 6 and were tested for significance at the a = 0.05 level. From a completed polynomial fit by determining the number of times j a completed polynomial fit on a given data window, the number that the first derivative of the fitted polynomial changes sign, of times j that the first derivative of the fitted polynomial changes indicating that j 1 polynomial degrees are required to reproduce sign was determined and the polynomial fit on the successive data the number of local maxima and minima, excluding end points, window is therefore required to test additional terms from k = 0 (14) D u m , B. S. Polynomial Regression. In Encyclopedia ofthe Statistical Sciences, to at least k = j 3 and until the two highest degrees Kotz, S.,Johnson, N. L., Eds.; Wiley: New York, 1986 Vol. 7, pp 700-703. tested failed the significance test or k = min{N - 2, ll}. A (15) Bevington, P. R Data Reduction and EworAnalysis for the Physical Sciences;
+
+
+
McGraw-Hill Book Co.: New York, 1969 Chapter 10. (16) Snedecor, G. W.; Cochran, W. G. Staristical Methods, 6th ed.; Iowa State University Press: Ames, IA,1967; Chapter 15.
(17) Press, W. H.; Flannery, B. P.;Teukolsky, S. A; Vetterling,W. T.Numerical Recipes; Cambridge University Press: Cambridge, U K , 1986; Chapter 7.
Analytical Chemistry, Vol. 67, No. 17, September 1, 1995 2759
0.9
I
0.8
n
0
0.7
f
n
0.6
5
0.5
.-
.-
0.4
E
Q, g
0.3
o0
.-c
0,
.-cc
E
I
1.0
0.9
I
0
Y CTI 0,
Z
0.2
a
Adaptive-Degree Polynomial Filter (ADPF)
5
7
9
11
15
19
23
Filter Length
\ sav;;itz; olayy
\
011
a,
.-%
Filter (Degree)
Adaptive Degree Polynomial Filter
0.7 -
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Filter Length / Peak Width
Figure 1. Noise reduction by ADPF and Savitzky-Golay filters applied to noise simulated by 1000 sets of normal deviates (u = 0, u = 1). Numbers adjacent to Savitzky-Golay results indicate equivalent polynomial degrees. Slope of -0.5 is indicated by the broken line.
Figure 2. Signal distortion, reported as peak height reduction, by ADPF and Savitzky-Golay filters applied to Gaussian curves. Numbers adjacent to Savitzky-Golay results indicate equivalent polynomial degrees.
VisualBasic-executable implementation is available from the author.
function,12 and of the degree of the polynomial used to fit the curve. When the filter length and polynomial degree are optimized with respect to the peak width of the signal, a polynomial of low degree may accurately reproduce the signal. If the polynomial degree is suboptimally low, the peak height returned by the polynomial filter may be significantly lower than the true value. Since many types of data analysis are based on peak height, reduction of peak height by data filtering is particularly undesirable. In the case of pure signal, adaptive-degree polynomial filtering at a = 0.05 reproduced the Gaussian signal height with high fidelity for all filter widths examined (Figure 2, filled circles), with less than 1%reduction in peak height even when the filter length was as much as 3 times greater than the peak width. At the same peak width, a standard Savitzky-Golay smooth of degree 2 or 3 would have led to a systematic peak height reduction of more than 30%. The fidelity of peak height returned by the adaptivedegree polynomial filter was achieved by the systematic selection of polynomial degree, ranging from 0 to 11,based on the results of an F test on the polynomial fits for the given filter length at the center of the peak. The argument for peak reproduction developed here using a Gaussian peak is believed to be generally true for any signal that can be approximated by polynomials. For example, sinusoidal functions may be expressed as series expansions, e.g., sin z = z - $ / 3 ! + $/5! - z7/7!+ .... Consequently, polynomial filters may fit sine functions very well if sufficient polynomial degrees are allocated (as does the adaptive-degree polynomial filter), with only truncation error left to consider. Using convolution weights hrso30, the transfer function and frequency response curve of polynomial filters may be calculated for any filter length and polynomial degree desired.18 The frequency response curve for adaptive-degree polynomial filters may be similarly calculated for the highest permitted polynomial degree, not to exceed 2m - 1 on a filter length of 2m + 1.
RESULTS AND DISCUSSION
Noise Reduction. For Savitzky-Golay filters applied to pure noise represented as a series of normal deviates, a log-log plot of relative standard deviation against filter length has a slope of approximately -0.5, indicating that the noise remaining after filtering varies inversely with the square root of the filter length (open circles, Figure 1). This relationship was previously shown both numerically12and analytically13for quadratic/cubic smooths, but is here seen to hold true for higher polynomial degrees up to 11. Clearly the choice of filter length and polynomial degree will determine the extent of removal of Gaussian noise, and for any filter length chosen, low polynomial degree is preferred for noise reduction. In the extreme case of pure noise, n = 0 or 1, corresponding to a running average and a sliding linear regression, respectively, provide the optimum smooths. Use of the adaptive-degree polynomial filter on normal deviates (Figure 1, filled circles) approaches the optimal level of noise reduction for pure noise with n = 0 or 1. The fact that the residual error using the adaptive degree polynomial filter is slightly higher than that for n = 0 or 1 is due to the occasional type I1 error, whereby the null hypothesis that the next higher term is not signi6cant is falsely accepted as true: this type of error must be accepted as an occasional occurrence so as to avoid an unacceptable level of type I error, whereby a true null hypothesis is rejected as false. Signal Fidelity. When Savitzky-Golay polynomial filters are applied to pure signal, the choice of inappropriately low polynomial degree can lead to failure to, accurately reproduce the signal, thereby introducing a systematic error quite different from suboptimal removal of statistical noise. The systematic distortion of a Gaussian curve by standard polynomial filters is shown in Figure 2 (open circles) as the relative peak height as a function of the ratio of the filter length to the peak width at half-height of a Gaussian peak (=2.35u), defined as the width of the smoothing 2760 Analytical Chemistry, Vo1. 67, No. 17, September 1, 7995
(18) Hamming, R. W. Digital W e n ,2nd ed.: Prentice-Hall: Englewood Cliffs, KJ, 1983: Chapter 3.
41
I
I
I
I
1.4
Savitzky-Golay
- Signal (composite 3
J
Gaussian curves) Signal t Gaussian noise (p=O, o=O.l)
‘I
0
h
4 B Y
‘
1.2
Adaptive Degree
0
1
.-
0.9
2
0.7
.E 0.8
.-8 0 z I
1
I
I
2
4
6
8
0.5
0.4
10
Time (arbitrary units)
Figure 3. Simulated chromatogram consisting of signal (solid line) and signal noise (points, one of five sets).
+
It may be noted that although Savitzky-Golay smoothing filters systematically distort peaks, the flattening of peaks is largely compensated by the broadening of peaks and the distortion of smoothed peak area is relatively small compared to the distortion of peak height.I2 Adaptive-degree polynomial filters will therefore offer little or no improvement to direct calculation of peak areas, although accurate peak shape is often useful in determining the beginning and end of peaks. Noise Reduction and Signal Fidelity. The results presented in Figures 1 and 2 for a standard polynomial filter on pure noise and pure signal demonstrate the opposing directions for optimal selection of polynomial degree for a given filter length: (1) reduction of polynomial degree to maximize noise reduction and (2) increase in the polynomial degree to maximize fidelity of signal reproduction. A poor match between filter length and polynomial degree leads to significant reduction in smoothed peak height if the polynomial degree is too low or suboptimal removal of statistical noise if the polynomial degree chosen is unnecessarily high. A compromise between these opposing trends is desirable. For purposes of evaluation, ADPF and Savitzky-Golay filters were applied to a synthetic function simulating a chromatographic separation in which the “true” signal and the magnitude of normally distributed, random noise were known (Figure 3). For Savitzky-Golay smooths, noise reduction (as measured by the root mean square of the smoothed value minus the true signal value) was highly dependent upon the choice of polynomial degree for each filter length applied to signal noise data sets (Figure 4,open circles). Multiple local minima in selection of filter length and filter degree are apparent for this data set, all leaving approximately 60%of the noise. For these five data sets, residual noise was 65% of the original value when a second or third degree polynomial was used with a filter length of five data points, decreasing to a minimum of 58% at longer filter windows for optimal polynomial degrees. Particularly troubling with the Savitzky-Golay filters is the observation that suboptimal choices of polynomial degree can systematically generate distortion greater than the amount of noise originally present ( U / Q 7 1) due to distortion of the signals
+
0.6
5
7
9
11
15
19
23
FiIter Length Figure 4. Noise reduction and signal fidelity by ADPF and Savitzky-Golay filters applied to five sets of synthetic chromatogram data in Figure 3. Error bars indicate standard deviation of the mean. Numbers adjacent to Savitzky-Golay results indicate equivalent polynomial degrees.
(Figure 4). Without prior knowledge of the synthetic function, selection of the optimum values of the filter length and filter degree would be essentially hit or miss. On the synthetic signal and noise data in Figure 3, the adaptivedegree polynomial filter was superior in noise reduction to any single combination of fixed polynomial degree and filter length (Figure 4). Improvement in noise reduction by adaptive degree polynomial filtering over the best matched fixed-degree filter ranged from 6% for a five-point smooth to 28% for a 23-point smooth, due to the larger amount of information about statistical noise in the longer filter lengths. Improvements in noise reduction by adaptive-degree filtering over suboptimal fixed polynomial degrees were even greater. From these results, it would appear that the use of ADPF both removes the guesswork in the choice of polynomial degree of the smoothing function and outperforms the best selection of a fixed polynomial degree. The adaptivedegree polynomial filter provided noise reduction at least equal to the optimal Savitzky-Golay smooth for short filter length and noise reduction beyond that of the optimal fixed-degree SavitzkyGolay smooth with all longer filter lengths. The improvement in noise reduction by ADPF compared to Savitzky-Golay filters with optimally-chosen polynomial degrees is due to the piecewise choice of polynomial degrees that are not identical to the fixed-degree choices of Savitzky-Golay optima. In Figure 5, a cumulative frequency plot of ADPF degree choices across the data set in Figure 3 has optimal Savitzky-Golay degrees highlighted for several filter lengths. In all cases, a signifcant percentage (60-90%) of the ADPF degree choices were lower than the optimum Savitzky-Golay filter, permitting greater removal of noise. In the case of 11-pointfilter lengths, additional improvement in a/uo was also obtained by permitting occasional higher degrees to better reproduce signal shape. Using ADPF, the variable polynomial degrees used throughout the data set cause variable levels of noise reduction, as previously shown in Analytical Chemistry, Vol. 67, No. 17, September 1, 1995
2761
x
1 -
-
Filter Length o 0.2 - /
/ ./
L
I
,
I
,
0
0
I1
”
A A
23
I’
I
-
5 points
,
-
/
,
I
Figure 1. As a consequence, if raw data contained uniform noise levels, the amounts of residual noise throughout the smoothed
data would be variable, depending upon the polynomial degree used to fit a particular window. Compared to other methods of implementing piecewise polynomial regressions, the main advantage claimed for the Savitzky-Golay method is its computational However, the polynomial regression, about which Savitzky and Golay5wrote, “(e)ven with a high-speed computer, it is a tedious proposition at best”, can now be executed on desktop computers at essentially no cost beyond the initial price of the computer. The use of orthogonal polynomials, such as Gram polynomials, remains advantageous for such calculations because it avoids many of the computational errors resulting from limits of computer precision due to small differences in numbers taken to high powers that otherwise easily occur using Gaussian polynomial regression techniques.lg Additionally, Gram polynomials are used by ADPF to speed calculation of the sum of squares of residuals (eq 5), permitting ,y2 to be evaluated as the sum of N multiplications rather than N(N + 1) multiplications (eq 4). The advantages of the computational speed and accuracy of the Savitzky-Golay method of data smoothing are retained in this adaptive-degree modification, while the piecewise and stepwise polynomial regression reduces the need to specify a priori polynomial degree and reduces noise to values lower than those achieved by the best SavitzkyGolay smooth. Received for review March 20, 1995. Accepted June 9,
1995.B AC950273Q
(19) Ralston, A. A First Course in Numerical Analysis; McGraw-Hill: New York, 1965; Chapter 6.
2762 Analytical Chemistry, Vol. 67, No. 17, September 1, 1995
* Abstract published
in Advance ACS Abstracts, August 1, 1995.
