A Savitzky–Golay filter is a digital filter that can be applied to a set of digital data points for the purpose of smoothing the data, that is, to increase the precision of the data without distorting the signal tendency. This is achieved, in a process known as convolution, by fitting successive sub-sets of adjacent data points with a low-degree polynomial by the method of linear least squares. When the data points are equally spaced, an analytical solution to the least-squares equations can be found, in the form of a single set of "convolution coefficients" that can be applied to all data sub-sets, to give estimates of the smoothed signal, (or derivatives of the smoothed signal) at the central point of each sub-set. The method, based on established mathematical procedures,. "Graduation Formulae obtained by fitting a Polynomial." was popularized by Abraham Savitzky and Marcel J. E. Golay, who published tables of convolution coefficients for various polynomials and sub-set sizes in 1964. Some errors in the tables have been corrected. The method has been extended for the treatment of 2- and 3-dimensional data.
Savitzky and Golay's paper is one of the most widely cited papers in the journal Analytical Chemistry and is classed by that journal as one of its "10 seminal papers" saying "it can be argued that the dawn of the computer-controlled analytical instrument can be traced to this article".
Applications The data consists of a set of points , where is an independent variable and is an observed value. They are treated with a set of convolution coefficients, , according to the expression
Selected convolution coefficients are shown in the tables, below. For example, for smoothing by a 5-point quadratic polynomial, and the -th smoothed data point, , is given by
where, , etc. There are numerous applications of smoothing, such as avoiding the propagation of noise through an algorithm chain, or sometimes simply to make the data appear to be less noisy than it really is.
The following are applications of numerical differentiation of data. Note When calculating the nth derivative, an additional scaling factor of may be applied to all calculated data points to obtain absolute values (see expressions for , below, for details).
Location of maxima and minima in experimental data curves. This was the application that first motivated Savitzky. The first derivative of a function is zero at a maximum or minimum. The diagram shows data points belonging to a synthetic Lorentzian curve, with added noise (blue diamonds). Data are plotted on a scale of half width, relative to the peak maximum at zero. The smoothed curve (red line) and 1st derivative (green) were calculated with 7-point cubic Savitzky–Golay filters. Linear interpolation of the first derivative values at positions either side of the zero-crossing gives the position of the peak maximum. 3rd derivatives can also be used for this purpose. Location of an end-point in a titration curve. An end-point is an inflection point where the second derivative of the function is zero. The titration curve for malonic acid illustrates the power of the method. The first end-point at 4 ml is barely visible, but the second derivative allows its value to be easily determined by linear interpolation to find the zero crossing. Baseline flattening. In analytical chemistry it is sometimes necessary to measure the height of an absorption band against a curved baseline. Because the curvature of the baseline is much less than the curvature of the absorption band, the second derivative effectively flattens the baseline. Three measures of the derivative height, which is proportional to the absorption band height, are the "peak-to-valley" distances h1 and h2, and the height from baseline, h3. Resolution enhancement in spectroscopy. Bands in the second derivative of a spectroscopic curve are narrower than the bands in the spectrum: they have reduced half-width. This allows partially overlapping bands to be "resolved" into separate (negative) peaks. The diagram illustrates how this may be used also for chemical analysis, using measurement of "peak-to-valley" distances. In this case the valleys are a property of the 2nd derivative of a Lorentzian. (x-axis position is relative to the position of the peak maximum on a scale of half width at half height). Resolution enhancement with 4th derivative (positive peaks). The minima are a property of the 4th derivative of a Lorentzian.
Moving average
The "moving average filter" is a trivial example of a Savitzky–Golay filter that is commonly used with time series data to smooth out short-term fluctuations and highlight longer-term trends or cycles. Each subset of the data set is fit with a straight horizontal line as opposed to a higher order polynomial. An unweighted moving average filter is the simplest convolution filter.
The moving average is often used for a quick technical analysis of financial data, like stock prices, returns or trading volumes. It is also used in economics to examine gross domestic product, employment or other macroeconomic time series.
It was not included in some tables of Savitzky-Golay convolution coefficients as all the coefficient values are identical, with the value .
Derivation of convolution coefficients When the data points are equally spaced, an analytical solution to the least-squares equations can be found. This solution forms the basis of the convolution method of numerical smoothing and differentiation. Suppose that the data consists of a set of n points (xj, yj) (j = 1, ..., n), where xj is an independent variable and yj is a datum value. A polynomial will be fitted by linear least squares to a set of m (an odd number) adjacent data points, each separated by an interval h. Firstly, a change of variable is made
where is the value of the central point. z takes the values where (e.g. m = 5 → z = −2, −1, 0, 1, 2).With even values of m, z will run from 1 − m to m − 1 in steps of 2 The polynomial, of degree k is defined asThe simple moving average is a special case with k = 0, Y = a0. In this case all convolution coefficients are equal to 1/m.
The coefficients a0, a1 etc. are obtained by solving the normal equations (bold a represents a vector, bold J represents a matrix).
where is a Vandermonde matrix, that is row of has values .
For example, for a cubic polynomial fitted to 5 points, z = −2, −1, 0, 1, 2 the normal equations are solved as follows.
Now, the normal equations can be factored into two separate sets of equations, by rearranging rows and columns, with
Expressions for the inverse of each of these matrices can be obtained using Cramer's rule
The normal equations become
and
Multiplying out and removing common factors,
The coefficients of y in these expressions are known as convolution coefficients. They are elements of the matrix
In general,
In matrix notation this example is written as
Tables of convolution coefficients, calculated in the same way for m up to 25, were published for the Savitzky–Golay smoothing filter in 1964, The value of the central point, z = 0, is obtained from a single set of coefficients, a0 for smoothing, a1 for 1st derivative etc. The numerical derivatives are obtained by differentiating Y. This means that the derivatives are calculated for the smoothed data curve. For a cubic polynomial
In general, polynomials of degree (0 and 1),Smoothing using the moving average is equivalent, with equally spaced points, to local fitting with a (sloping) straight line (2 and 3), (4 and 5) etc. give the same coefficients for smoothing and even derivatives. Polynomials of degree (1 and 2), (3 and 4) etc. give the same coefficients for odd derivatives.
Algebraic expressions It is not necessary always to use the Savitzky–Golay tables. The summations in the matrix JTJ can be evaluated in closed form,
so that algebraic formulae can be derived for the convolution coefficients.The expressions given here are different from those of Madden, which are given in terms of the variable m' = (m − 1)/2. Functions that are suitable for use with a curve that has an inflection point are:
Smoothing, polynomial degree 2,3 : (the range of values for i also applies to the expressions below) 1st derivative: polynomial degree 3,4 2nd derivative: polynomial degree 2,3 3rd derivative: polynomial degree 3,4
Simpler expressions that can be used with curves that don't have an inflection point are: Smoothing, polynomial degree 0,1 (moving average): 1st derivative, polynomial degree 1,2:
Higher derivatives can be obtained. For example, a fourth derivative can be obtained by performing two passes of a second derivative function.
Use of orthogonal polynomials An alternative to fitting m data points by a simple polynomial in the subsidiary variable, z, is to use orthogonal polynomials.
where P0, ..., Pk is a set of mutually orthogonal polynomials of degree 0, ..., k. Full details on how to obtain expressions for the orthogonal polynomials and the relationship between the coefficients b and a are given by Guest. Expressions for the convolution coefficients are easily obtained because the normal equations matrix, JTJ, is a diagonal matrix as the product of any two orthogonal polynomials is zero by virtue of their mutual orthogonality. Therefore, each non-zero element of its inverse is simply the reciprocal of the corresponding element in the normal equation matrix. The calculation is further simplified by using recursion to build orthogonal Gram polynomials. The whole calculation can be coded in a few lines of PASCAL, a computer language well-adapted for calculations involving recursion.
Treatment of first and last points Savitzky–Golay filters are most commonly used to obtain the smoothed or derivative value at the central point, z = 0, using a single set of convolution coefficients. (m − 1)/2 points at the start and end of the series cannot be calculated using this process. Various strategies can be employed to avoid this inconvenience. The data could be artificially extended by adding, in reverse order, copies of the first (m − 1)/2 points at the beginning and copies of the last (m − 1)/2 points at the end. For instance, with m = 5, two points are added at the start and end of the data y1, ..., yn. y3,y2,y1, ... ,yn, yn−1, yn−2. Looking again at the fitting polynomial, it is obvious that data can be calculated for all values of z by using all sets of convolution coefficients for a single polynomial, a0 .. ak. For a cubic polynomial Convolution coefficients for the missing first and last points can also be easily obtained. This is also equivalent to fitting the first (m + 1)/2 points with the same polynomial, and similarly for the last points.
Weighting the data It is implicit in the above treatment that the data points are all given equal weight. Technically, the objective function
being minimized in the least-squares process has unit weights, wi = 1. When weights are not all the same the normal equations become
If the same set of diagonal weights is used for all data subsets, , an analytical solution to the normal equations can be written down. For example, with a quadratic polynomial,
An explicit expression for the inverse of this matrix can be obtained using Cramer's rule. A set of convolution coefficients may then be derived as
Alternatively the coefficients, C, could be calculated in a spreadsheet, employing a built-in matrix inversion routine to obtain the inverse of the normal equations matrix. This set of coefficients, once calculated and stored, can be used with all calculations in which the same weighting scheme applies. A different set of coefficients is needed for each different weighting scheme.
It was shown that Savitzky–Golay filter can be improved by introducing weights that decrease at the ends of the fitting interval.
Two-dimensional convolution coefficients