Anal. Chem. 1990, 62, 2749-2752
tation can be easily eliminated by reducing the size of the sample cell and by utilizing high-speed pumps for sample delivery/drain systems. The main point is that the proposed system can be easily constructed and is quite inexpensive.
ACKNOWLEDGMENT We thank W. B. Sisson, J. A. Easterling, and W. F. Syner of the Food and Drug Administration (FDA) for their assistance during this work.
LITERATURE CITED Inductively Coupled Plasm in Analytical Atomic Spectronmtry ; Montaser, A., Gollghtly, D. W., Eds., VCH: New York, 1987; 660 pages. Gustavsson, A. I n Inductively Coupled Plasma in Analytical Atomic Spectfometfy; Montaser, A., Gollghtly, D. W., Eds.; VCH: New York, 1987. Browner, R. F. I n Inductively Coupled "aEmission Spectroscopy; Boumans. P. W. J. M., Ed.; Wlley: New York, 1987; Part 11. West, C. D.; Hume, D. N. Anal. Chem. 1964, 36, 412-415. Wendt. R . H.; Fassel, V. A. Anal. Chem. 1965, 37, 920-922. Olsen, K. W.; Haas, W. J.; Fassel, V. A. Anal. Chem. 1977, 49, 632-637. Fassel, V. A.; Bear, B. R. Spectrochim. Acta 1986, 418, 1089-1113. Tang, Y. Q.; Trassy, C. Spectrochim. Acta 1986, 418, 143-150. Walters, P. E.; Barnardt. C. A. Spectrochim. Acta 1988, 438, 325-337. Petruccl, G. A.; Van Loon, J. D. Spectrochim. Acta 1990, 458. 959-968.
2749
(11) Instrument description from Applied Research Laboratorles, Inc., 24911 Stanford Ave, Valencla, CA 91355. (12) Instrument description from Baird Cap., Analytical Instruments Dlvision, 125 Mlddlesex Turnpike, Bedford, MA 01730-1468. (13) Instrument description from Cetac Technologies, Inc., P.O. Box 31118, 302 South 36th St, Omaha, NE 68131. (14) Clifford, R. H.; Montaser, A. A LowGost Ultrasonic Nebulizer for Plasma Spectrometry. 16th Annual Meeting of the Federation of Analytical Chemistry and Spectroscopy Societies, Chicago, IL. Oct 1989; paper no. 320. (15) Qlnhan, J.; Chu, 2 . ; Brushwyler, K.; Hleftje, G. M. Appl. Spectrosc. 1980. 4 4 , 183-186 and references therein. (16) Montaser, A.; Clifford, R. H.; Sinex, S.A.; Capar, S.G. Anal. Chem. 1989, 61, 2589-2592. (17) Montaser, A.; Clifford, R. H.; Sinex, S. A.; Capar, S. G. J . Anal. At. Spectrom. 1989, 4 , 499-503. (18) Montaser, A.; Ishil, I.; Tan, H.; Clifford, R. H.; Golightly, D. W. Spectrochim. Acta 1989, 448, 1163-1169. (19) Clifford, R. H.; Ishll, I.; Montaser, A.; Meyer, G. A. Anal. Chem. 1990, 62, 390-394. (20) Bachalo, W. D.; Houser, M. J. Opt. f n g . 1980, 23, 583-590. (21) . . Lana R. J. Ultfasonic Atomization of Liauhls. J . Acoust. SOC. Am. 196% 34, 6-8. (22) Clifford, R. H.; Montaser, A.; Sinex, S. A,; Capar, S.G. Anal. Chem. 1989, 61 2777-2784. I
RECEIVED for review August 8,1990. Accepted September 11, 1990. This research was sponsored in part by the US. Department of Energy under Grant No. DE-FG05-87-ER-13659. Partial support for A.M. and R.H.C. and use of facilities were provided by the FDA.
Predicting the Statistical Properties of Least-Squares Polynomial Filters G.R. Phillips and J. M. Harris* Department of Chemistry, University of Utah, Salt Lake City, Utah 84112 Experimental data generally consist of random noise superimposed on an unknown deterministic signal. Leastsquares polynomial filters, introduced by Savitzky and Golay (I),are frequently used for reducing such noise in time-dependent data and for numerical differentiation. Errors in the original article have been corrected by Steinier et al. (2) and Madden (3). In addition to hundreds of reported applications, several fundamental studies have addressed the selection of optimal smoothing parameters (4-8), such as the polynomial order, filter length, and recursive filter application. The reduction of random noise by digital filtering is independent of the form of the signal and depends on the design of the filter. Enke and Nieman (4) and Leach et al. (9) have studied noise reduction by using numerical simulation. In simulation studies, a deterministic signal is repeatedly combined with numerically generated random error and smoothed. Analysis of several hundred such records allows generalizations to be made about the average or expected behavior of the filter. This process must be repeated for each combination of filter parameters to completely characterize the filter. Such studies are time-consuming, and the results are limited to the conditions studied. The results of such studies are generally summarized in graphical form, making extraction of precise values more difficult than from algebraic expressions. Alternatively, a propagation-of-errors approach @ , l o ,l l ) can also be used, in which the ratio of the variance after smoothing divided by the variance of the raw data is equated to the sum of the squares of the smoothing coefficients. This approach does not provide covariance terms, and the accuracy of the results depends on the numerical methods used to estimate the convoluting coefficients. Steiner et al. (2),Leach et al. (9),and Bialkowski (10) have discussed the use of matrices to construct Savitzky-Golay 0003-2700/90/0362-2749$02.50/0
polynomial filters. Bialkowski (IO), in addition, has used a matrix approach with propagation-of-errors to determine the effects of polynomial order on the precision of filtered data. In the present work, symbolic programming software is used together with the matrix formulation of polynomial filters to predict the precision of filtered results without the need for time-consuming numerical simulation. The elements of the variance-covariance matrix provide both the variance of the amplitude and derivative estimates along with the covariance between these. The analytical expressions reveal the theoretical origins of the dependence of noise reduction on filter width and polynomial order. We demonstrate this approach by deriving expressions of the noise reduction by quadratic/cubic filters and compare their predictions to numerical simulations of the effects of filter width.
THEORY A matrix formulation for polynomial filters is introduced and used to predict the error structure of the filtered results. The unfiltered data must consist of a sequence of equally spaced data points with noise of constant variance
Y1 = s, + 5 (1) Here y , denotes an experimental observation, s, denotes the noise-free signal, and c, denotes random noise or error. The general expression for a digital filter used to estimate the midpoint of a ( 2 m + 1) point data segment is m
g1=0
=
c cy1
(2)
i=-m
where yl and 3, are the original and filtered data points, respectively, and the cI)s are the filter coefficients. The values of these coefficients are determined by the filter design (e.g. functional model, filter symmetry, etc.). 0 1990 American Chemical Society
2750
ANALYTICAL CHEMISTRY, VOL. 62, NO. 24, DECEMBER 15, 1990
Chart I
3(3m2+ 3m - 1) (2m - 1)(2m + 1)(2m + 3)
0
0
3 m(m + 1)(2m + 1) 0
v = Q2
-15 (2m - 1)(2m + 1)(2m + 3 )
If the underlying functional form of the signal, s, were known, regression techniques should be used in place of filtering (4). In the absence of such knowledge, segments of 2m + 1 data points can be modeled by an nth order polynomial (n C 2m + 1) n
fi
=
C&jil j=O
(3)
The 2m + 1 point data window is slid along the curve, and the abscissa reindexed for each data segment. Evaluation of eq 3 a t i = 0 reveals the following relationships between the coefficients P n j and the value of the filtered data point and its derivatives a t the midpoint fi=o
= Pn,o
(4)
Hence, Savitzky-Golay polynomial least-squares filtering is equivalent to the determination of the P,,&’s. Combining eqs 1 and 3 yields 2m + 1 simultaneous equations, one for each data point, y , n
yi =
E&jij
+ ti
j=O
Since n < 2m + 1, this is an overdetermined problem and can be easily solved by using matrix least squares (12). The set of simultaneous equations in eq 6 can be expressed in matrix notation, where lower-case letters capped by an arrow, -, denote column vectors and upper-case bold face letters denote matrices
Y’=xp+7 (7) In this equation, y’ contains the (2m + 1)observations of the dependent variable, Z is a vector of (2m + 1)random errors,
and the vector P contains the n + 1unknown parameters. The design matrix, X, has dimensions (2m + 1) by ( n + 1) and contains the (ij) factors from eq-6. The least-squares estimate of 6, denoted g, is obtained by minimizing the sum of squares of error and is given by (12)
6 = (XtX)-lXt$ The estimated error in 6 is contained
(8)
in the variance-co-
variance matrix (12)
where no2is the variance of the noise in the raw data. The first and second diagonal elements in V contain the expected value of the variance of LI,,,~ and bn,l, respectively. It is clear from eqs 2 and 8 that the matrix T = (XtX)-’Xt contains the elements C,,i of the filter function, C,. Convoluting coefficients for estimating the value and the first derivative at the midpoint of the filtered segment are contained in the first and second rows of T. Once these coefficients have been calculated and tabulated for a given polynomial order and filter length, no further matrix operations are required. The unsmoothed data are convoluted with these coefficients to yield the filtered data. The efficiency of this process often overshadows the underlying theoretical basis for computing
-15 (2m - 1)(2m + 1)(2m + 3) 0 45 --
(2m - l)m(m + 1)(2m + 1)(2m + 3)
I
(12)
the smoothing coefficients. However, the equivalence between this convolution and least-squares polynomial regression over a limited range of data yields estimates not only for the filtered data and its derivatives but also for the precision of these results.
EXPERIMENTAL SECTION Analytical expressions for the relative precision of symmetric quadratic filters were derived algebraically by hand. Formulas for higher order symmetric filters and initial point filters were derived by using the REDUCE algebraic programming system (Rand Corp., Santa Monica, CA) operating on a VAX 8300 computer. These formulas were shown to agree with the results of numerical matrix calculations performed by using a PC version of MATLAB (The Mathworks, Inc., South Natick, MA). The precision of data filtered by quadratic polynomial filters was evaluated by using simulated data consisting of a linear combination of Gaussian-distributed noise with linear and exponential decays. The magnitude of the noise was varied from 5 to 50% of the initial value of the test function, while the filter width was varied from 5 to 25 points. For each test function and noise magnitude, between 150 and 500 data sets were simulated and analyzed, using symmetric and initial point filters. Relative errors were calculated for each filter width as the standard deviation of the b,,;~ (for a series of the simulated data sets) divided by the standard deviation of the noise.
RESULTS AND DISCUSSION The results presented in this paper are for single pass quadratic and cubic filters. Multiple pass filters can be analyzed by convoluting and single pass filter with itself. For simulation studies, the standard deviation of the noise was varied over an order of magnitude without affecting the results, as predicted by eq 9. Symmetric (Midpoint) Filters. Symmetric filters are designed to estimate the middle point in a data window and are the most frequently studied and applied polynomial filters. The abscissa is indexed from -m to +m, with the filter width equal to 2m + 1. The elements of the design matrix, x i j , are given by
(10)
x Y. . = ii
where i = -m, ..., m and j = 0, ..., n. The awkward indexing results from the convention for numbering filter coefficients. This can be clarified by writing the elements in matrix form; for example, the design matrix for a quadratic polynomial filter is
I1
m
m 2
I
The diagonal elements of variance-covariance matrix contain the variance of the filtered amplitude and derivative estimates, while the off-diagonal elements predict the covariance between them. The matrix (XtX) for a quadratic filter is easily inverted, using Cramer’s rule, to yield the variance-covariance matrix given in Chart I. Symbolic manipulation produces exact results, in contrast to numerical methods which are sensitive to the numerical stability of algorithms used to evaluate the results (10). Equations for the relative standard deviation of filtered values and first
ANAL.YTICAL CHEMISTRY, VOL. 62, NO. 24, DECEMBER 15, 1990
Table I. Noise Reduction by a Symmetric Savitzky-Golay Filter of Width w
,
('/ao)2n,oa
= 4(w - 2)(w)(w
P)
a
('/u0)23J
v, U
-0 .-I
,-. .05
+ 3)
= 2 or 3,
,
,
,
,
,
, ,
,
,
,
,
,
,,
, ,,,
-
n N
bo
2
v
.02 1
.01
L
,
.
A
.
~
.
!
.
~
.
l
.
O
.
~
5
'
x
,
,,,,~,~ml
-20 1
.E .10 0 n
O n
,
-0
+ 2)
25(3w4 - 1 8 ~ + ' 31) = (w - 3)(w - 2)(w - l ) w ( w + l ) ( w + 2)(w
,
,
'\
Relative Error 3(3w2 - 7)
,
2751
.
~
.
~
1
. .
10
'
.
3
....'...LLLLul 20
Filter Width Flgure 2. Standard deviation ( U / U ~ ) ~of, , the first derivative relative to the raw data versus filter width, for a symmetric quadratic filter. Theoretical predictions (- - -) and simulation results (e)are shown on a log-log plot.
'
*'
0.9 1
.-0 a ..-
4
L
Y
C
0.8
0 .w
C -
R
Figwe 1. Standard deviation ( U / U ~ )of~ the , ~ smoothed mkipoint relative to the raw data versus filter width, for a symmetric quadratic filter (n = 2). Theoretical predictions (- - -) and simulation results from ref 4 ( * ) and this laboratory (e)are shown on a log-log plot.
derivatives versus w = 2m + 1 are given in Table I. Expressions for the convoluting coefficients have been derived by Madden (3). The relative standard deviation of the filtered point from a quadratic filter versus the filter width is shown in Figure 1 on a log-log scale. Theoretical predictions, shown as a dashed line, are in excellent agreement with simulation results from Enke and Nieman (4)and this laboratory. A linear least-squares fit of the log of relative standard deviation versus the log of filter width has a slope of -0.51, in agreement with the observation by Enke, that the noise reduction varies inversely with the square root of the filter width. This is also apparent in the net (w-'/~) dependence a t large w of the analytical expressions for ( U / U , ) ~ ,in~ Table I. The second diagonal element of V (eq 12) predicts the variance of the first derivative filter. Dividing this value by the noise variance, uo2,and taking the square root provides a theoretical prediction of the relative standard deviation of the first derivative filter, which is compared with numerical simulations in Figure 2. A linear least-squares fit of log,, ( u / ~ , ) versus ~ , ~ loglo (w)has a slope of -1.51, in agreement with the ( w - ~ /dependence ~) (at large w ) of (u/uo)z,l predicted in Table I. Initial Point Filters. A modified version of the Savitzky-Golay algorithm has been developed to estimate the initial point and initial slope of a data set(9, 13-15). Initial point filters have heen used to determine initial light intensities in thermal lens absorption measurements (16,17) and can be used to estimate initial reaction rates. Khan (14)and Wentzell et al. (11) have studied systematic errors that can be introduced if these filters are carelessly applied. Again, a matrix formalism enables prediction a priori of the statistical prop-
0.7
-
n
bo
\
s 0.6 : 5
10
20
Filter Width Figure 3. Standard deviation ( O / U ~ )in~ the , ~ smoothed value relative to the raw data versus filter width, for an initial point quadratic filter. Theoretical predictions (- - -) and simulation results (e)are shown on a log-log plot. erties of these filters. Initial point filters smooth or differentiate the first point in a data segment, in contrast to symmetric filters, which estimate the filtered value or derivative at the middle point. For initial point filters, the abscissa is numbered from 0 to 2m. Defining w = 2m + 1, the design matrix contains the elements x i j = ij, where i = 0, ..., w - 1and j = 0, ...,n. The XtX and (XtX)-l matrices are much more complicated for initial point fiiters than for midpoint filters because the design matrix is not symmetric about 0. Closed-form expressions have been obtained by using the REDUCE program for symbolic manipulation and are summarized in Table 11. The relative standard deviation in the initial point is shown in Figure 3 as a function of the filter width. The theoretical and simulated results are consistent and show a changing power dependence on w over the range shown. From the expression for ( ~ / C T , ) ~in, ,Table 11, one would expect a net (w-'I2) dependence at very large w. For the last seven points in Figure 3, a least-squares fit of log,, ( U / U ~ ) ~ versus ,O log,, (w) produces a slope of -0.58, which is close to the limiting value. Results for the precision of initial slope filter are shown in Figure 4. A linear least-squares fit of log (u/uo)z,l versus log w over the entire range of w in the plot has a slope of -1.46,
2752
ANALYTICAL CHEMISTRY, VOL. 62,
NO.24,
DECEMBER 15, 1990
Table 11. Properties of Initial Point Second and Third Order Polynomial Filter of Width w Convoluting Coefficients
3(3w2- 3w + 2) - 18(2w - 1)i + 30i2 C2,O.l = w(w + l ) ( w + 2)
c2,lJ. = - 1 8 ( 2 ~- l ) ( w - 2)(w - 1) + 1 2 ( 2 -~ 1 ) ( 8-~1l)i - 1 8 0 ( -~ l)i2 (w - 2)(u - l ) w ( w + l ) ( w + 2) 8(2w - 1)(w2- w + 3) - 20(6w2 - 6~ + 5)i + 1 2 0 ( 2 ~- l)i2 - 140i3 '3,0,1 = w(w + l ) ( w + 2)(w + 3) -2O(w - 3 ) ( -~ 2 ) ( -~ l)(6w2 - 6 w + 5) - 200(6w4 - 27w3 + 42w2 - 3 0 +~ 1l)i c3,1,,= -3OO(w - 1)(3w- 5)(3w - 2)i2 + 280((6w2 - 15w + ll)? (U!- 3 ) ( w - 2)(w - l ) w ( w + l ) ( w + 2)(w + 3) Relative Error
12(2w - 1 ) ( 8 -~ 11) = (w - 2)(w - l ) w ( w + l ) ( w
('/uO)2zJ
8 ( 2 -~ 1)(w2 - 10 200(6w4 - 27w3 ('/uO)23,1
.
=
\
8 \
\q
LITERATURE CITED
*\ r:
0.1
t.~.~.~.~.~.~.~.~.'.~ I
5
,
,
10
,
+ 42w2 - 3 0 +~ 11) + l ) ( w + 2)(w + 3)
prevention of signal distortion. The effect of smoothing on signal distortion can be determined by convolution of the functional form of the signal with the filter function. The effect of the filter on noise reduction has previously been studied by using numerical simulation, but the method is time-consuming and the accuracy of the results depends on the number of trials. In the present study, symbolic programming software has produced analytical expressions for the relative precision of polynomial filters obtained from a matrix formulation of the regression problem.
'.,
*\
+ 3)
(w - 3)(w - 2)(w - l ) w ( w
C
\
+ 2)
b
. , ( , . . I . . . .
20
Filter Width
Flgure 4. Standard deviation ( u / u & , in the first derivative relative to the raw data versus filter width, for an initial point quadratic filter.
Theoretical predictins (- - -) and simulation results (0)are shown on a log-log plot.
which is in agreement with the ( W - ~ I ~ dependence ) (at large w ) of (u/u& predicted in Table I and very similar to the dependence found for midpoint slope filters. It should be noted that these results measure the standard deviation of the derivative relative to the standard deviation of the noise. From numerical simulations, Leach et al. (9) reported the standard deviation of the filtered derivative relative to that of a simple two-point derivative, which cannot be directly compared to the results in Figure 4. Conclusion. The choice of optimal smoothing parameters represents a compromise between noise reduction and the
(1) Savitzky, A.; Golay, M. J. E. Anal. Chem. 1984, 36, 1627-1639. (2) Steiner, J.; Termonia, Y.: Dekour, J. Anal. Chem. 1972, 4 4 , 1906- 1909. (3) Madden, H. M. Anal. Chem. 1978, 50, 1383-1366. (4) Enke, C. G.:Nieman, T. A. Anal. Chem. 1978, 48, 705A-712A. (5) Betty, K. R.; Horlick, G. Anal. Chem. 1977, 48, 351-352. (6) Kaiser, J. F.; Reed, W. A. Rev. Scl. Insfrum. 1977, 48,1447-1457. (7) Edwards, T. H.; Wilson, P. D. Appl. Spectrosc. 1974, 28, 541-545. (8) Bromba, M. U. A.; Ziegler, H. Anal. Chem. 1981, 53, 1583-1586. (9) Leach, R. A.; Carter, C. A.; Harris, J. M. Anal. Chem. 1984, 56, 2304-2307. (10) Bialkowski, S. E. Anal. Chem. 1989, 61, 1308-1310. (11) Wentzell. P. D.; Doherty,T. P.; Crouch, S. R. Anal. Chem. 1987, 59, 367-37 1. (12) Deming, S. N.; Morgan, S. L. Clln. Chem. 1979, 2 5 , 840-855. Draper, N. R.; Smith, H. App/ied Regression Analysis, 2nd ed.: John Wiley 8 Sons: New York, 1981; Chapter 2. (13) Baedecker, P. A. Anal. Chem. 1985. 57, 1477-1479. (14) Khan, A. Anal. Chem. 1988, 6 0 , 369-371. (15) Proctor, A.; Sherwood, P. M. A. Anal. Chem. 1980, 52, 2315-2321. (16) Carter, C. A.; Harris, J. M. Anal. Chem. 1983, 55. 1256-1261. (17) Leach, R . A.: Harris, J. M. Anal. Chem. 1984, 56, 1481-1487.
RECEIVED for review July 24, 1990. Accepted September 13, 1990. This research was supported in part by the National Science Foundation under Grants CHE85-06667 and CHE90-10319.
