Savitzky-Golay Methods

Three options on the Spectrum Processing dialog window utilise a simple but effective means of extracting information from data containing noise. Two of these options are common place, where smoothing of experimental data and the determinations of derivatives are performed using the algorithm proposed by Savitzky and Golay (A. Savitzky and M. J. E. Golay, Anal. Chem., 36, 1627 (1964)). The third option utilises the same principle behind the smoothing and differentiation algorithm to perform integration on a set of data and represents an extension of the Savitzky-Golay approach.

Figure 1

Theory

Given a set of data containing both signal and noise, the initial objective of the Savitzky-Golay method is to replace the raw data by a smoother set of data representing the true signal responsible for the record intensities. For a constant underlying signal, the most natural means of estimating the true signal from a set of measurements would be to average the values. The act of averaging a set of values is in fact one example, and possibly, the simplest applications of the linear-least-squares principle. Spectral data, on the other hand, typically contain peaks superimposed on a background signal and therefore a more subtle use of averaging is required if the essential structure in the data set is to be retained. One way to use the averaging process, but to maintain information relating to the variation in the intensities, is to perform a local averaging for each bin within a spectrum, for example, each data bin could be replaced by the average of three bins, itself and the two bins either side of the bin. The averaging operation could be applied to a data set via a digital convolution of the data bins with a convolution kernel consisting of the values {1/3, 1/3, 1/3}. These simply, yet often used concepts are at the basis of the Savitzky-Golay method, which in essence applies the least-squares principle to determine an improved set of kernel coefficients for use in a digital convolution, where these improved coefficients are determined, in the least-squares sense, using polynomials rather than, for the case of averaging, simply assuming a constant value determined from a sub-range of data bins. Indeed, the Savitzky-Golay method could be seen as a generalisation of averaging data, since averaging a sub-range of data corresponds to using a Savitzky-Golay polynomial of degree zero.

To illustrate the Savitzky-Golay method, consider the specific example in which five data bins are used to approximate a quadratic polynomial. The polynomial can be expressed in the form:

.

Where the coefficients a_{0 ,} a₁ and a₂ are determined from the simultaneous equations in which the abscissa x is the index for the data bin; the origin is always placed at the central data bin and so the abscissa values corresponding to each of the data bins are {-2, -1, 0, 1, 2}:

Or

.

Where the evenly spaced data bins {d_-2, d_-1, d₀, d₁, d₂} are selected with the target of replacing the value for d₀ with the value for the polynomial at x = 0 or p(0) = a₀. Since there are five equations and only three unknowns, the coefficients to the polynomial must be determined in the least-squares sense, where the linearly independent basis functions are 1, x and x² , therefore the normal equations yield:

.

Since A^TA is a square symmetric matrix of rank three, the coefficient vector a is determined from [A^TA]^-1A^T, the top row of which yields the prescription for computing the value of a₀, namely:

Thus, for each set of five such data bins, the central bin can be replaced by the value determined for a₀, or in other words, a digital convolution using the five point kernel {s_i} and the raw data bins results in a smoothed set of data bins, where a linear least squares quadratic polynomial is used to model the data, five channels at a time.

Similarly, the derivative of a spectrum can be computed using the Savitzky-Golay polynomial. Again the intention is to approximate the derivative at a given point in the spectrum using the derivative of the polynomial at x = 0. Since dp(0)/dx = a₁, the second row of the matrix [A^TA]^-1A^T yields a second convolution kernel for computing the derivative of the spectrum and, apart from the difference in the kernel values, the computation of the derivative proceeds in an analogous fashion to that of the smoothing calculation.

The variation on a theme occurs with the integration algorithm. Here, the polynomial is again used in the approximation to the integral of the spectrum, however rather than evaluating the polynomial at the point x = 0 the integral between x = -1 and x = 0 is computed. Completing the integration over a range of data bin is then achieved by summing the set of sub-intervals computed using a convolution kernel determined from all three rows the matrix [A^TA]^-1A^T:

.

Therefore

Where the i_j in the convolution Kernel are the linear combination of the three terms from each of the rows of [A^TA]^-1A^T prescribed by the expression above. To compute the integral with respect to energy requires the multiplication of the convolution integrals by the energy step-size.

The use of the Savitzky-Golay polynomial to evaluate the derivative is essential for noisy data, since an interpolating polynomial would produce an inferior result to that of a least-squares approximation. On the other hand, there is an academic nature to using Savitzky-Golay integration, since simply summing the data bins includes the averaging mechanism required to reduce the influence of random noise on the computed values; however the symmetry of using the same algorithm to compute both the derivative and its inverse is somewhat pleasing.

The theory behind the Savitzky-Golay options, whilst illustrated using a five-point quadratic approximation, can be applied for any odd number of points applied to any degree of polynomial, provided the number of data bins over specifies the information required to compute the polynomial coefficients. CasaXPS offers quadratic, quartic and linear polynomials, while the number of data bins used to process the spectra is specified in a text-field on both the Smooth and also the Differentiate property pages on the Spectrum processing dialog window. Integration is always performed using a five-point quadratic Savitzky-Golay approximation. One further point regarding integration is that integration cannot recover constant offsets lost during differentiation, so there is an option (Zero Adjust) on the Integration property page, which when ticked, ensures the data is all positive following an integration calculation.