Signal-to-Noise Ratio Enhancement by Least-Squares Polynomial

Instrumentation

Signal-to-Noise Ratio Enhancement by Least-Squares Polynomial Smoothing Every real-world experimenter must cope with noise and the complexities it adds to the measurement step. To reduce the noise content of our mea­ sured signals, we employ (at acquisi­ tion time) such techniques as shielding, electrical filters, signal aver­ aging, or integration. But as comput­ ing power permeates our laboratories, there is an increasing tendency also to employ numerical methods to "mas­ sage" the data for additional noise re­ duction, "after-the-fact". For this pur­ pose, least-squares polynomial smoothing has become a popular tech­ nique, as evidenced by frequent refer­ ence in the literature (1-17). An experimenter must confront himself with two questions regarding smoothing: "Do I want to smooth my d a t a ? " "How should I adjust the pa­ rameters of the smoothing process?" By the time one has gotten far enough to ask the first question, he has usual­ ly already decided that the answer is yes. This decision is reached on the basis of its ease of implementation once the data and the smoothing pro­ grams are already in the computer rather than on the basis of what is to be gained (or lost) through smoothing. The question of smoothing parameters is all too often answered by intuition

C. G. Enke Department of Chemistry Michigan State University East Lansing, Mich. 48824

Timothy A. Nieman School of Chemical Sciences University of Illinois Urbana, III. 61801

"Least-squares smoothing is of cosmetic value only. Smoothed data do not contain any additional information" or blind guess rather than by consider­ ation of firm guidelines. This paper has been written in an attempt to clar­ ify the capabilities and limitations of the smoothing process and to aid in the selection of the optimum smooth­ ing parameters. Smoothing Process To smooth raw data, a polynomial can be used to approximate local sec­ tions of data. The fitted polynomial is

Figure 1 . S m o o t h i n g n o r m a l v s . u n i f o r m noise

then used to calculate a "better" value for the central point. Savitzky and Golay (18) and Steiner et al. (19) pop­ ularized this technique by publishing tables of weighting coefficients to be convoluted with the raw data to yield smoothed values. This approach is given by Equation 1

where Y and Y* are the raw and smoothed data, respectively, the C, are the weighting coefficients, Ν is a normalization coefficient, and Yj-m, . . . , Yj+m is the range of raw data points used in providing a smoothed estimate of the central point, Yj. Their method assumes that the data points are evenly spaced on the abscis­ sa axis, the data curves are continu­ ous, and that the uncertainty in the ordinate values is much greater than the uncertainty in the abscissa values. This type of smooth reduces noise (high-frequency fluctuations) by act­ ing as a low pass filter (16). Informa­ tion carried by the signal, within the remaining bandwidth, will be re­ tained. Since the smoothing function will only be an approximate represen­ tation of each local section of data, the true signal will suffer some distortion

Figure 2 . Noise r e d u c t i o n w i t h m u l t i p l e pass s m o o t h i n g

Solid line: normal distribution; dashed line: uniform distribution. Single pass smoothing used

ANALYTICAL CHEMISTRY, VOL. 48, NO. 8, JULY 1976 · 705 A

which will depend on the details of the smoothing process and on the proper­ ties of the data being smoothed: All techniques for improvement of the signal-to-noise ratio (S/N) of a set of data produce this improvement at the expense of some other property of the data. Analog filters cause band­ width reduction and asymmetry in time. Signal averaging results in a loss of time resolution and/or bandwidth but does not cause signal distortion. The loss due to these signal condition­ ing types of procedures is determined at the time of data acquisition. The smoothing technique of data condi­ tioning results in bandwidth reduction and in signal distortion due to loss in resolution along the abscissa (wave­ length axis of a spectrum). However, the type and degree of compromise (between noise reduction and signal distortion) depend on parameters of the smoothing process and may be de­ cided after data acquisition is com­ pleted. For a Savitzky and Golay smooth, the parameters which may be adjusted include the type of smooth­ ing function (quadratic-cubic, quarticquintic, derivative, etc.), the width of the smoothing interval (the number of points in the smoothing function), and the number of repetitions of the smoothing process. The effect of these parameters on the degree of noise re­ duction and signal distortion has been determined by an empirical study (20), the major results of which are described below. Discussion will be limited to quadratic-cubic smooths containing 5-23 points and from 1 to 128 repetitions. Noise Reduction Study

Noise Source Considerations and Single-Pass Smoothing. In any real situation, signal and noise would exist together and be smoothed together. To facilitate study, it was advanta­ geous to separately generate noise-free signals and random noise and to smooth each separately. Since the smoothing method studied is linear, this approach yielded the same results as when signal and noise are added to­ gether before smoothing. We can assume that the major noise component of a system is white noise. White noise is a mixture of signals of all frequencies with random ampli­ tudes and phase (21). The probability distribution of the amplitudes is nor­ mal (Gaussian) if the observed noise is the sum of several independent noise sources (22). In this case, our noise generator can be any algorithm that yields a random sequence of numbers with a normal distribution of magni­ tudes. Since least-squares smoothing makes no assumption concerning the probability distribution of the uncer­ tainties (noise) other than it has a 706 A ·

mean of zero (23), one may well won­ der why the concern for a normal dis­ tribution in the model noise source. The reason is shown in Figure 1, which illustrates a log-log plot of noise re­ duction (for both uniform and normal distributions) vs. the width of the smoothing function. The degree of noise reduction is measured as the ratio of the standard deviation after (σ) and before (σ0) a single pass of the smoothing function. For normally dis­ tributed noise, the amount of noise re­ maining after a single pass smooth is inversely proportional to the square root of the width of the smoothing function, as indicated by Savitzky and Golay (18). However, the results for uniformly distributed noise show sig­ nificant deviation from the squareroot relationship and decrease at a faster rate than predicted. Although the least-squares smoothing process makes no assumption concerning the form of the noise distribution, the out­ come of the smooth is seen to depend upon the distribution. Studies employing synthetic noise with a nonnormal distribution may yield results inconsistent with real world situa­ tions. Creation of a normally distributed random sequence begins with a uni­ formly distributed random sequence as generated by a congruential algo­ rithm (Equation 2) contained in many computer software systems (24). Ui = Uc-^dnod

D)

(2)

The Ui are the elements of the ran­ dom sequence, and M is a multiplier. The modulus D is usually the comput­ er word length. These uniform values can be converted into standard normal values (μ = 0, σ = 1) either by using a look-up table constructed from a stan­ dard normal table or by using the Central Limit theorem approximation (Equation 3), where Ε and V are the expected value and variance, respec­ tively

(25). A sequence of numbers with specified mean and standard deviation can be generated as in Equation 4. Effect of Multiple Smooths. Mul­ tiple passes of a smoothing function can be used to reduce data irregulari­ ties to a further extent than with sin­ gle-pass smoothing. Figure 2 illus­ trates σ/σ0 of white noise as a function of the number of times the smooth was performed for several different width smooths. Although a log-log plot is not quite linear, the fraction of noise remaining after multiple smooths varies approximately with the inverse eighth root of the number

A N A L Y T I C A L C H E M I S T R Y , V O L . 4 8 , N O . 8, J U L Y

1976

Figure 3. Peak distortion due to smoothing a: Single pass smoothing (subscript G: Gaussian peaks; subscript L: Lorentzian peaks), b: Height response factor for Gaussian peaks (number by curves: smoothing ratio)

of smooths. As a result, the first pass of a particular smoothing function produces the most noticeable noise re­ duction. Given enough passes, a short smooth can reduce noise to the same extent as a longer smooth. But for a limited number of passes, and consid­ ering only noise reduction and not sig­ nal distortion, choice of as long a smoothing function as feasible is more important than the number of times the smooth is repeated. For a 2 η + 1 point Savitzky and Golay smooth, there are η points on each end of the data array which can­ not be replaced by smoothed values. For this reason, some workers studying repeated smoothing have dropped η points from each end of the data array after each smooth (10). This leads to rapid loss of data points, as five passes of a 23-point smooth would drop 110 points. The opposite extreme, to retain all smoothed and unsmoothed points, leads to unsatis­ factory distortion at the ends of the smoothed sections. The best proce-

dure to eliminate end distortion and yet maintain the maximum number of data points is to use a shorter smooth on the ends than on the bulk of the data. To reduce programming com­ plexity, a satisfactory compromise can be reached in which the unsmoothed end points are set equal to the end smoothed point after each pass of the smoothing function. Signal Distortion Study

The ideas in the previous section concerning noise reduction were de­ veloped without any consideration of the fate of the analytical signal as a re­ sult of the smoothing. Since the smoothing function will reduce the bandwidth and also will not be an exact fit for the noise-free signal, dis­ tortion will occur. Isolated peak-shaped signals (Gaussian and Lorentzian) were useful models for this study due to their widespread occurrence in signals from many types of data sources (spectros­ copy, chromatography, and electro­ chemistry). A Gaussian-shaped peak can be generated by y = (h) exp [-(χ - μ)2/(2 σ2))

(5)

where μ is the mean, σ is the standard deviation, and h is the peak height. For a spectroscopist or chromatographer, it is more convenient to express peak width in terms of the full-widthat-half maximum (FWHM) instead of σ. Since σ2 = (FWHM)2/8(\n

2)

(6)

y = (h) exp [-4(ln 2) Χ (χ - ^/(FWHM)'1]

(7)

for a Gaussian peak, then

Similarly, a Lorentzian peak may be generated by y =h

(FWHM)2/{(FWHM)2 + 4(χ-μ)2]

(8)

[For the purposes of the discussion in this paper, peak widths are defined in terms of the number of times the sig­ nal is sampled across the peak FWHM. The time variable of sam­ pling rate is thus replaced by the spa­ tial variable of data point number (or core memory address).] As peaks of these types are smoothed, they become lower and broader. The degree of this distortion is determined not by the width of the smoothing function but by the ratio of the width of the smoothing function and the peak FWHM (hereafter called the "smoothing ratio"). Figure 3a il­ lustrates the percentage change in the height and area of Gaussian and Lo­ rentzian peaks as a function of the smoothing ratio for a single pass smooth. Peak height is reduced, and area increased by smoothing. The

Gaussian peak suffers less peak height distortion but more area distortion than the Lorentzian peak. Results for a Voigt profile would be intermediate between the Lorentzian and Gaussian. To keep peak height distortion below 1%, the smoothing ratio must be less than 0.9 for a Gaussian peak and less than 0.7 for a Lorentzian peak. As long as the smoothing ratio is less than unity, the peak is not noticeably broadened. Figure 3b illustrates the percentage change in peak height for Gaussian peaks due to multiple pass smoothing. The variable on the ordinate is called the smoothing response factor (/) and is the value of the peak parameter after smoothing divided by the value before smoothing. On the graph, the family of curves represents different values of the smoothing ratio as indi­ cated next to each curve. As would be expected, the same response can be achieved by several combinations of smoothing ratio and number of repeti­ tions. Curves for height, width, and area are similar, and, in general, any combinations of smoothing ratio and number of repetitions which result in approximately the same value of the peak height response factor (//, ) also yield the same values for the peak area (fa ) and peak width (fw ) response fac­ tors. As the degree of smoothing is in­ creased to eliminate more noise, the parameters of peaks composing the signal suffer increased distortion. Since these parameters contain infor­ mation concerning the chemical sys­ tem being studied, smoothing results in a change in the relationship be­ tween the peak parameter and its as­ sociated physical property. Peak height before smoothing (h„) might be directly proportional to some physical property (p) as h„ = kp. Since height after smoothing (h) is related to h0 through the height response factor (//, ), the new relationship becomes h = kpfi,- The factor //, can be taken from graphs such as in Figure 3b and has been seen to be a function of the width and number of repetitions of the

708 A · ANALYTICAL CHEMISTRY, VOL. 48, NO. 8, JULY 1976

smoothing function and of the FWHM of the peak. If all of the peaks in the signal being smoothed are of equal FWHM, then the peak heights will be reduced proportionally. However, if the peaks are of different widths (as in a chromatogram or IR spectrum), and if the same smooth is used over the entire data array, each peak will be distorted to a different extent. To cor­ rect this disproportionate reduction, each peak must be treated separately according to its FWHM. For a com­ plex spectrum this could be a practical impossibility due to overlapping cor­ rections from adjacent peaks. The al­ ternatives are either to smooth in such a way that the maximum observed error is within acceptable limits or to use different smoothing functions on different peaks so that the smoothing ratio is a constant (4). Signal-to-Noise Ratio Enhancement

The two elements of the study de­ scribed above considered first pure noise and then noise-free signals. In real measurements the signal and noise will be present together, and the signal-to-noise ratio (S/N) is then used to describe the quality of the measurement. The larger the value of S/N, the more readily is information extracted from the data when human interpretation is involved at any step. The results from smoothing pure noise and noise-free peaks can be combined to yield information concerning the ef­ fect smoothing has on the S/N of a signal and to learn how to use smooth­ ing to maximize S/N. Since peak height is often used to quantitate a chemical measurement, results from peak height reduction of Gaussian peaks will be used to develop the S/N relationships. The data in Figure 3b are the ratios of peak height after and before smoothing (h/ha), whereasthe data in Figure 2 are the ratios of the root-mean-square (rms) noise after and before smoothing (σ/ σ„). Dividing h/h0 by σ/συ yields (hi σ)Ι(Η0/σ0), or (S/N)/(S/N)0, which is the factor of S/N enhancement result-

Figure 4. S/N en­ hancement for 16-point wide peak (number by curves: repetition)

Figure 5. Maximum S/N enhancement factor vs. peak width

ing from smoothing. In this case, S/N is defined as the ratio of the observed (not true) signal and the rms noise. Due to smoothing, the observed signal may be altered from its true form. This discussion does not consider dis­ tortion to be noise since the distortion is not random and is predictable. Figure 4 illustrates the S/N en­ hancement factor as a function of the width of the smoothing function for a typical peak. In this case, the peak FWHM is 16 points. T h e graph con­ tains a family of curves with different numbers of repetitions of smoothing. There is a definite limit to the degree of S/N enhancement through smooth­ ing, and t h a t limit can be reached by several combinations of width of smoothing function and number of repetitions. If this type of study is re­ peated on peaks of different widths, the limit of S/N enhancement in­ creases (Figure 5) as the FWHM of the peak increases (approximately as a function of FWHM06). It is important to notice t h a t the degree of S/N en­ hancement is generally small. For a peak sampled by even as many as 32 points at its FWHM, only a fourfold factor of S/N enhancement is possi­ ble. In addition, smoothing to the maximum enhancement produces ap­ proximately 12% distortion in both the peak height and width; distortion in the area is limited to less than 1%. It is an awkward procedure to have to consult graphs like t h a t in Figure 4 for each different peak width to deter­ mine the proper smoothing procedure. In addition, the need to vary both the smoothing width and number of repe­ titions increases the complexity of cal­ culation. If we limit ourselves to single pass smoothing, a simple guideline arises. Figure 6 illustrates a plot of the fraction of the maximum S/N en­ hancement obtained vs. the smoothing ratio using single pass smoothing. T h e

Figure 6. Fraction of maximum S/N enhancement for single pass smoothing

definite maximum at a smoothing ratio of about 2 is noteworthy. For sin­ gle pass smoothing, the maximum S/N enhancement is obtained by using a smoothing function twice as wide as the peak FWHM. T h e S/N enhancement obtained by this method is within 5% of the maximum S/N en­ hancement obtained if one allows the possibility of more than one pass of a smoothing function. If this study is re­ peated with Lorentzian-shaped peaks, the same result (choice of smoothing ratio equal to 2) is reached for qua­ dratic-cubic smoothing. It is possible to predict this empiri­ cally obtained guideline from the theory of matched filters and the gen­ eralized mathematics for the construc­ tion of the polynomial smoothing functions. For measurements of peak height, the best filter is the matched filter (26-28) in which the FWHM of the filter matches the FWHM of the peak to be smoothed. T h e problem now reduces to determining the rela­ tionship between the number of points in the smoothing function and its FWHM. For a quadratic-cubic smooth, the general formula for the weighting coefficients, C,, is given by Equation 9 (29)

where the nomenclature is the same as in Equation 1 (2 m + 1 is the number of points in the smooth). From this equation, one can derive the following expression for the desired ratio:

T h e value of this ratio varies from 1.92 for a 5-point smooth (m = 2) to 1.83

710 A · ANALYTICAL CHEMISTRY, VOL. 48, NO. 8, JULY 1976

for any smoothing function of 17 points (m = 8) or wider. (The limiting value for infinite m is ν Ί θ / 3 or 1.826.) In addition to the quadratic-cubic smoothing function analysis described above, it is possible to apply the same approach to higher order smoothing functions. For example, consideration of quartic-quintic smoothing func­ tions results in the guideline t h a t S/N is maximized by using a smoothing ratio of about 2.9 to 3. Application of the criterion for qua­ dratic-cubic is shown in Figure 7. T h e raw data show a pair of lines (each sampled by 8 points at FWHM) out of the spectrum of a mercury lamp as re­ corded by a vidicon spectrometer (30). T h e source was attenuated to produce a signal of low S/N. T h e data after a 5-point smooth (smoothing ratio = 0.6) have suffered negligible distortion of the peaks b u t still contain signifi­ cant noise. A 23-point smooth (smoothing ratio = 2.9) produces a sig­ nal considerably free of noise b u t at the expense of significant peak distor­ tion. T h e intermediate case, using a 15-point smooth (smoothing ratio = 1.9) as would be chosen using the smoothing ratio of 2 (or 1.83) guide­ line, produces the desired compromise between noise reduction and signal distortion. For comparison, a 15-point smoothing function (Figure 8) has a FWHM of 8.17, confirming t h a t the appropriate function was one having the same FWHM as the data peak. Fourier transform approaches to digital filtering are gaining increasing popularity (28). T h e type of smooth­ ing function considered in this paper, in addition to many other types of smoothing functions, can be applied in the Fourier domain. Guidelines devel­ oped in this paper are relevant to Fou­ rier filtering (when the frequency ver­ sion of a polynomial smoothing func­ tion is used); however, it is first neces-

of s m o o t h i n g . As long as t h e u s e r m a i n t a i n s p e r s p e c t i v e on his r e a s o n s for d a t a collection a n d analysis, leastsquares polynomial smoothing can be a useful tool. Acknowledgment W e t h a n k B r i a n H a h n for a s s i s t a n c e in t h e p r e l i m i n a r y stages of t h i s work a n d S a n d r a N i e m a n for helpful advice a n d discussions d u r i n g t h e e n t i r e course of t h e project. References

Figure 7. Example of smoothing to maximize S/N enhancement

s a r y t o t r a n s f o r m t h e g u i d e l i n e s from t h e t i m e d o m a i n to t h e frequency d o ­ main. Conclusions A t t h i s p o i n t it is n e c e s s a r y t o real­ ize t h a t least-squares smoothing is of cosmetic value only. If t h e informa­ t i o n c o n t a i n e d in d a t a is t o be e x t r a c t ­ ed using l e a s t - s q u a r e s fitting b y a n automated procedure with no h u m a n i n t e r a c t i o n , t h e n p r i o r u s e of leasts q u a r e s s m o o t h i n g does n o t r e s u l t in i m p r o v e d precision. S m o o t h e d d a t a d o n o t c o n t a i n a n y a d d i t i o n a l informa­ tion; in fact, s o m e i n f o r m a t i o n h a s b e e n lost d u e t o b a n d w i d t h r e d u c t i o n . S m o o t h i n g " c l e a n s u p " noisy d a t a t o m a k e i n f o r m a t i o n in t h e d a t a m o r e easily accessible t o h u m a n i n t e r p r e t a ­ t i o n . F o r m e a s u r e m e n t s involving h u m a n interaction with data interpre­ tation, the optimum smooth produces t h e m a x i m u m S/N e n h a n c e m e n t of t h e raw d a t a . W h e n s m o o t h i n g t o t h i s

e x t e n t , t h e r e is a n o b v i o u s d e g r e e of d i s t o r t i o n i n t r o d u c e d i n t o t h e signal. T h e u s e r m u s t d e c i d e if, for his p u r ­ poses, t h e loss of a c c u r a c y is t o l e r a b l e . If s m o o t h i n g is u s e d t o p r e p a r e d a t a for d i s p l a y for h u m a n i n t e r p r e t a t i o n , t h e n S/N e n h a n c e m e n t m o r e t h a n c o m p e n s a t e s for signal d i s t o r t i o n . Similarly, if p a r a m e t e r s are t o b e m a n u a l l y e x t r a c t e d , a n d if t h e analysis involves only relative m e a s u r e m e n t s against standards t h a t have been t r e a t e d identically, t h e n signal distor­ t i o n factors will cancel (provided t h e p e a k s are e q u a l in w i d t h ) a n d t h e in­ creased S/N will yield a h i g h e r preci­ sion m e a s u r e m e n t . If a b s o l u t e m e a ­ s u r e m e n t s a r e r e q u i r e d or if p e a k b r o a d e n i n g causes r e s o l u t i o n t o d e ­ crease below a n a c c e p t a b l e level, t h e n it is n o t d e s i r a b l e t o s m o o t h t o t h e m a x i m u m S/N. In t h i s case, a b e t t e r a p p r o a c h w o u l d b e t o use d e g r e e of d i s t o r t i o n a n d n o t S/N e n h a n c e m e n t as t h e c r i t e r i o n for j u d g i n g t h e e x t e n t

Figure 8. Weighting coefficients for 15point quadraticcubic smoothing function

712 A · ANALYTICAL CHEMISTRY, VOL. 48, NO. 8, JULY 1976

(1) R. R. Ernst, "Sensitivity Enhancement in Magnetic Resonance", in "Advances in Magnetic Resonance", Vol 2, pp 3450, J. S. Waugh, Ed., Academic Press, New York, N.Y., 1966. (2) H. P. Yule, Anal. Chem., 38, 103 (1966). (3) H. P. Yule, Trans. Am. Nucl. Soc, 10, 38 (1967). (4) H. P. Yule, Nucl. Instrum. Methods, 54,61 (1967). (5) R. N. Jones, R. Venkataraghavan, and J. W. Hopkins, Spectrochim. Acta, 23A, 941 (1967). (6) K. Yamashita and S. Minami, Jpn. J. Appl. Phys., 8, 1505 (1969). (7) J. P. Porchet and Hs. H. Gunthard, J. Phys., E, 3, 261 (1970). (8) M. Caprini, S. Cohn-Sfetcu, and A. Manof, IEEE Trans., AU-18, 389 (1970). (9) K. Yamashita and S. Minami, Jpn. J. Appl. Phys., 10, 1097 (1971). (10) H. P. Yule, Anal. Chem., 44, 1245 (1972). (11) M.J.E. Golay, IEEE Trans., C21, 299 (1972). (12) H. Tominaga, M. Dojyo, and M. Tanaka, Nucl. Instrum. Methods, 98, 69 (1972). (13) T. H. Edwards and P. D. Willson, Appl. Spectrosc, 28, 541 (1974). (14) G. Dulaney, Anal. Chem., 47, 24A (1975). (15) M. W. Overton, L. L. Alber, and D. E. Smith, ibid., ρ 363Α. (16) P. D. Willson and T. H. Edwards, Appl. Spectrosc. Rev., in press. (17) G. Giannelli and O. Altimura, Rev. Sci. Instrum., 47, 27 (1976). (18) A. Savitzky and M.J.E. Golay, Anal. Chem., 36, 1627 (1964). (19) J. Steiner, Y. Termonia, and J. Deltour, ibid., 44, 1906 (1972). (20) T. A. Nieman, PhD thesis, Michigan State University, East Lansing, Mich., 1975. (21) H. V. Malmstadt, C. G. Enke, S. R. Crouch, and G. Horlick, "Optimization of Electronic Measurements", Module 4, pp 8, 9, Benjamin, Menlo Park, Calif., 1974. (22) W. C. Hamilton, "Statistics in Physi­ cal Sciences", ρ 67, Ronald Press, New York, N.Y., 1964. (23) W. C. Hamilton, ibid., ρ 130. (24) Τ. Ε. Hull and D.D.F. Day, "Comput­ ers and Problem Solving", pp 137-39, Addison-Wesley, Ontario, Canada, 1970. (25) M. Dwass, "Probability and Statis­ tics", pp 297, 318, 329-30, Benjamin, New York, N.Y., 1970. (26) G. L. Turin, IRE Trans. Inf. Theory, IT-6, 311 (1960). (27) G. Horlick, Anal. Chem., 44, 943 (1972). (28) K. R. Betty and G. Horlick, Appl. Spectrosc, 30, 23 (1976). (29) E. Whittaker and G. Robinson, "The Calculus of Observation", pp 291-96, 4th éd., Dover, New York, N.Y., 1967. (30) T. A. Nieman and C. G. Enke, Anal. Chem., 48,619(1976).