Savitzky–Golay filter

The data consists of a set of points {x_j, y_j}, j = 1, ..., n, where x_j is an independent variable and y_j is an observed value. They are treated with a set of m convolution coefficients, C_i, according to the expression

${\displaystyle Y_{j}=\sum _{i={\tfrac {1-m}{2}}}^{\tfrac {m-1}{2}}C_{i}\,y_{j+i},\qquad {\frac {m+1}{2}}\leq j\leq n-{\frac {m-1}{2}}}$

Selected convolution coefficients are shown in the tables, below. For example, for smoothing by a 5-point quadratic polynomial, m = 5, i = −2, −1, 0, 1, 2 and the jth smoothed data point, Y_j, is given by

${\displaystyle Y_{j}={\frac {1}{35}}(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+2})}$,

where, C₋₂ = −3/35, C₋₁ = 12 / 35, etc. There are numerous applications of smoothing, which is performed primarily to make the data appear to be less noisy than it really is. The following are applications of numerical differentiation of data.[8] Note When calculating the nth derivative, an additional scaling factor of ${\displaystyle {\frac {n!}{h^{n}}}}$ may be applied to all calculated data points to obtain absolute values (see expressions for ${\displaystyle {\frac {d^{n}Y}{dx^{n}}}}$, below, for details).

(1) Synthetic Lorentzian + noise (blue) and 1st derivative (green)

(2) Titration curve (blue) for malonic acid and 2nd derivative (green). The part in the light blue box is magnified 10 times

(3) Lorentzian on exponential baseline (blue) and 2nd derivative (green)

(4) Sum of two Lorentzians (blue) and 2nd derivative (green)

(5) 4th derivative of the sum of two Lorentzians

1. Location of maxima and minima in experimental data curves. This was the application that first motivated Savitzky.[4] 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.
2. 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.[9] 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.
3. Baseline flattening. In analytical chemistry it is sometimes necessary to measure the height of an absorption band against a curved baseline.[10] 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.[11]
4. 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.[12] 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).
5. Resolution enhancement with 4th derivative (positive peaks). The minima are a property of the 4th derivative of a Lorentzian.

When the data points are equally spaced, an analytical solution to the least-squares equations can be found.[2] 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 (x_j, y_j) (j = 1, ..., n), where x_j is an independent variable and y_j 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

${\displaystyle z={{x-{\bar {x}}} \over h}}$

where ${\displaystyle {\bar {x}}}$ is the value of the central point. z takes the values ${\displaystyle {\tfrac {1-m}{2}},\cdots ,0,\cdots ,{\tfrac {m-1}{2}}}$ (e.g. m = 5 → z = −2, −1, 0, 1, 2).[note 1] The polynomial, of degree k is defined as

${\displaystyle Y=a_{0}+a_{1}z+a_{2}z^{2}\cdots +a_{k}z^{k}.}$[note 2]

The coefficients a₀, a₁ etc. are obtained by solving the normal equations (bold a represents a vector, bold J represents a matrix).

${\displaystyle {\mathbf {a} }=\left({{\mathbf {J} }^{\mathbf {T} }{\mathbf {J} }}\right)^{-{\mathbf {1} }}{\mathbf {J} }^{\mathbf {T} }{\mathbf {y} },}$

where ${\displaystyle \mathbf {J} }$ is a Vandermonde matrix, that is ${\displaystyle i}$-th row of ${\displaystyle \mathbf {J} }$ has values ${\displaystyle 1,z_{i},z_{i}^{2},\dots }$.

For example, for a cubic polynomial fitted to 5 points, z= −2, −1, 0, 1, 2 the normal equations are solved as follows.

${\displaystyle \mathbf {J} ={\begin{pmatrix}1&-2&4&-8\\1&-1&1&-1\\1&0&0&0\\1&1&1&1\\1&2&4&8\end{pmatrix}}}$

${\displaystyle \mathbf {J^{T}J} ={\begin{pmatrix}m&\sum z&\sum z^{2}&\sum z^{3}\\\sum z&\sum z^{2}&\sum z^{3}&\sum z^{4}\\\sum z^{2}&\sum z^{3}&\sum z^{4}&\sum z^{5}\\\sum z^{3}&\sum z^{4}&\sum z^{5}&\sum z^{6}\\\end{pmatrix}}={\begin{pmatrix}m&0&\sum z^{2}&0\\0&\sum z^{2}&0&\sum z^{4}\\\sum z^{2}&0&\sum z^{4}&0\\0&\sum z^{4}&0&\sum z^{6}\\\end{pmatrix}}={\begin{pmatrix}5&0&10&0\\0&10&0&34\\10&0&34&0\\0&34&0&130\\\end{pmatrix}}}$

Now, the normal equations can be factored into two separate sets of equations, by rearranging rows and columns, with

${\displaystyle \mathbf {J^{T}J} _{\text{even}}={\begin{pmatrix}5&10\\10&34\\\end{pmatrix}}\quad \mathrm {and} \quad \mathbf {J^{T}J} _{\text{odd}}={\begin{pmatrix}10&34\\34&130\\\end{pmatrix}}}$

Expressions for the inverse of each of these matrices can be obtained using Cramer's rule

${\displaystyle (\mathbf {J^{T}J} )_{\text{even}}^{-1}={1 \over 70}{\begin{pmatrix}34&-10\\-10&5\\\end{pmatrix}}\quad \mathrm {and} \quad (\mathbf {J^{T}J} )_{\text{odd}}^{-1}={1 \over 144}{\begin{pmatrix}130&-34\\-34&10\\\end{pmatrix}}}$

The normal equations become

${\displaystyle {\begin{pmatrix}{a_{0}}\\{a_{2}}\\\end{pmatrix}}_{j}={1 \over 70}{\begin{pmatrix}34&-10\\-10&5\end{pmatrix}}{\begin{pmatrix}1&1&1&1&1\\4&1&0&1&4\\\end{pmatrix}}{\begin{pmatrix}y_{j-2}\\y_{j-1}\\y_{j}\\y_{j+1}\\y_{j+2}\end{pmatrix}}}$

and

${\displaystyle {\begin{pmatrix}a_{1}\\a_{3}\\\end{pmatrix}}_{j}={1 \over 144}{\begin{pmatrix}130&-34\\-34&10\\\end{pmatrix}}{\begin{pmatrix}-2&-1&0&1&2\\-8&-1&0&1&8\\\end{pmatrix}}{\begin{pmatrix}y_{j-2}\\y_{j-1}\\y_{j}\\y_{j+1}\\y_{j+2}\\\end{pmatrix}}}$

Multiplying out and removing common factors,

${\displaystyle {\begin{aligned}a_{0,j}&={1 \over 35}(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+2})\\a_{1,j}&={1 \over 12}(y_{j-2}-8y_{j-1}+8y_{j+1}-y_{j+2})\\a_{2,j}&={1 \over 14}(2y_{j-2}-y_{j-1}-2y_{j}-y_{j+1}+2y_{j+2})\\a_{3,j}&={1 \over 12}(-y_{j-2}+2y_{j-1}-2y_{j+1}+y_{j+2})\end{aligned}}}$

The coefficients of y in these expressions are known as convolution coefficients. They are elements of the matrix

${\displaystyle \mathbf {C=(J^{T}J)^{-1}J^{T}} }$

In general,

${\displaystyle (C\otimes y)_{j}\ =Y_{j}=\sum _{i=-{\tfrac {m-1}{2}}}^{\tfrac {m-1}{2}}C_{i}\,y_{j+i},\qquad {\frac {m+1}{2}}\leq j\leq n-{\frac {m-1}{2}}}$

In matrix notation this example is written as

${\displaystyle {\begin{pmatrix}Y_{3}\\Y_{4}\\Y_{5}\\\vdots \end{pmatrix}}={1 \over 35}{\begin{pmatrix}-3&12&17&12&-3&0&0&\cdots \\0&-3&12&17&12&-3&0&\cdots \\0&0&-3&12&17&12&-3&\cdots \\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\ddots \end{pmatrix}}{\begin{pmatrix}y_{1}\\y_{2}\\y_{3}\\y_{4}\\y_{5}\\y_{6}\\y_{7}\\\vdots \end{pmatrix}}}$

Tables of convolution coefficients, calculated in the same way for m up to 25, were published for the Savitzky–Golay smoothing filter in 1964,[3][5] The value of the central point, z = 0, is obtained from a single set of coefficients, a₀ for smoothing, a₁ 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

${\displaystyle {\begin{aligned}Y&=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}=a_{0}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {dY}{dx}}&={\frac {1}{h}}\left({a_{1}+2a_{2}z+3a_{3}z^{2}}\right)={\frac {1}{h}}a_{1}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {d^{2}Y}{dx^{2}}}&={\frac {1}{h^{2}}}\left({2a_{2}+6a_{3}z}\right)={\frac {2}{h^{2}}}a_{2}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {d^{3}Y}{dx^{3}}}&={\frac {6}{h^{3}}}a_{3}\end{aligned}}}$

In general, polynomials of degree (0 and 1),[note 3] (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.

Two-dimensional smoothing and differentiation can also be applied to tables of data values, such as intensity values in a photographic image which is composed of a rectangular grid of pixels.[17] [18] Such a grid is referred as a kernel, and the data points that constitute the kernel are referred as nodes. The trick is to transform the rectangular kernel into a single row by a simple ordering of the indices of the nodes. Whereas the one-dimensional filter coefficients are found by fitting a polynomial in the subsidiary variable z to a set of m data points, the two-dimensional coefficients are found by fitting a polynomial in subsidiary variables v and w to a set of the values at the m × n kernel nodes. The following example, for a bivariate polynomial of total degree 3, m = 7, and n = 5, illustrates the process, which parallels the process for the one dimensional case, above.[19]

${\displaystyle v={\frac {x-{\bar {x}}}{h(x)}};w={\frac {y-{\bar {y}}}{h(y)}}}$

${\displaystyle Y=a_{00}+a_{10}v+a_{01}w+a_{20}v^{2}+a_{11}vw+a_{02}w^{2}+a_{30}v^{3}+a_{21}v^{2}w+a_{12}vw^{2}+a_{03}w^{3}}$

The rectangular kernel of 35 data values, d₁ − d₃₅

v w	−3	−2	−1	0	1	2	3
−2	d1	d2	d3	d4	d5	d6	d7
−1	d8	d9	d10	d11	d12	d13	d14
0	d15	d16	d17	d18	d19	d20	d21
1	d22	d23	d24	d25	d26	d27	d28
2	d29	d30	d31	d32	d33	d34	d35

becomes a vector when the rows are placed one after another.

d = (d₁ ... d₃₅)^T

The Jacobian has 10 columns, one for each of the parameters a₀₀ − a₀₃, and 35 rows, one for each pair of v and w values. Each row has the form

${\displaystyle J_{\text{row}}=1\quad v\quad w\quad v^{2}\quad vw\quad w^{2}\quad v^{3}\quad v^{2}w\quad vw^{2}\quad w^{3}}$

The convolution coefficients are calculated as

${\displaystyle \mathbf {C} =\left(\mathbf {J} ^{T}\mathbf {J} \right)^{-1}\mathbf {J} ^{T}}$

The first row of C contains 35 convolution coefficients, which can be multiplied with the 35 data values, respectively, to obtain the polynomial coefficient ${\displaystyle a_{00}}$, which is the smoothed value at the central node of the kernel (i.e. at the 18th node of the above table). Similarly, other rows of C can be multiplied with the 35 values to obtain other polynomial coefficients, which, in turn, can be used to obtain smoothed values and different smoothed partial derivatives at different nodes.

Nikitas and Pappa-Louisi showed that depending on the format of the used polynomial, the quality of smoothing may vary significantly.[20] They recommend using the polynomial of the form

${\displaystyle Y=\sum _{i=0}^{p}\sum _{j=0}^{q}a_{ij}v^{i}w^{j}}$

because such polynomials can achieve good smoothing both in the central and in the near-boundary regions of a kernel, and therefore they can be confidently used in smoothing both at the internal and at the near-boundary data points of a sampled domain. In order to avoid ill-conditioning when solving the least-squares problem, p < m and q < n. For software that calculates the two-dimensional coefficients and for a database of such C's, see the section on multi-dimensional convolution coefficients, below.

The idea of two-dimensional convolution coefficients can be extended to the higher spatial dimensions as well, in a straightforward manner,[17][21] by arranging multidimensional distribution of the kernel nodes in a single row. Following the aforementioned finding by Nikitas and Pappa-Louisi[20] in two-dimensional cases, usage of the following form of the polynomial is recommended in multidimensional cases:

${\displaystyle Y=\sum _{i_{1}=0}^{p_{1}}\sum _{i_{2}=0}^{p_{2}}\cdots \sum _{i_{D}=0}^{p_{D}}(a_{i_{1}i_{2}\cdots i_{D}}\times u_{1}^{i_{1}}u_{2}^{i_{2}}\cdots u_{D}^{i_{D}})}$

where D is the dimension of the space, ${\displaystyle a}$'s are the polynomial coefficients, and u's are the coordinates in the different spatial directions. Algebraic expressions for partial derivatives of any order, be it mixed or otherwise, can be easily derived from the above expression.[21] Note that C depends on the manner in which the kernel nodes are arranged in a row and on the manner in which the different terms of the expanded form of the above polynomial is arranged, when preparing the Jacobian.

Accurate computation of C in multidimensional cases becomes challenging, as precision of standard floating point numbers available in computer programming languages no longer remain sufficient. The insufficient precision causes the floating point truncation errors to become comparable to the magnitudes of some C elements, which, in turn, severely degrades its accuracy and renders it useless. Chandra Shekhar has brought forth two open source softwares, Advanced Convolution Coefficient Calculator (ACCC) and Precise Convolution Coefficient Calculator (PCCC), which handle these accuracy issues adequately. ACCC performs the computation by using floating point numbers, in an iterative manner.[22] The precision of the floating-point numbers is gradually increased in each iteration, by using GNU MPFR. Once the obtained C's in two consecutive iterations start having same significant digits until a pre-specified distance, the convergence is assumed to have reached. If the distance is sufficiently large, the computation yields a highly accurate C. PCCC employs rational number calculations, by using GNU Multiple Precision Arithmetic Library, and yields a fully accurate C, in the rational number format.[23] In the end, these rational numbers are converted into floating point numbers, until a pre-specified number of significant digits.

A database of C's that are calculated by using ACCC, for symmetric kernels and both symmetric and asymmetric polynomials, on unity-spaced kernel nodes, in the 1, 2, 3, and 4 dimensional spaces, is made available.[24] Chandra Shekhar has also laid out a mathematical framework that describes usage of C calculated on unity-spaced kernel nodes to perform filtering and partial differentiations (of various orders) on non-uniformly spaced kernel nodes,[21] allowing usage of C provided in the aforementioned database. Although this method yields approximate results only, they are acceptable in most engineering applications, provided that non-uniformity of the kernel nodes is weak.

1. The sum of convolution coefficients for smoothing is equal to one. The sum of coefficients for odd derivatives is zero.[25]
2. The sum of squared convolution coefficients for smoothing is equal to the value of the central coefficient.[26]
3. Smoothing of a function leaves the area under the function unchanged.[25]
4. Convolution of a symmetric function with even-derivative coefficients conserves the centre of symmetry.[25]
5. Properties of derivative filters.[27]

Signal distortion and noise reduction

It is inevitable that the signal will be distorted in the convolution process. From property 3 above, when data which has a peak is smoothed the peak height will be reduced and the half-width will be increased. Both the extent of the distortion and S/N (signal-to-noise ratio) improvement:

• decrease as the degree of the polynomial increases
• increase as the width, m of the convolution function increases

Effect of smoothing on data points with uncorrelated noise of unit standard deviation

For example, If the noise in all data points is uncorrelated and has a constant standard deviation, σ, the standard deviation on the noise will be decreased by convolution with an m-point smoothing function to[26][note 5]

polynomial degree 0 or 1:${\displaystyle {\sqrt {1 \over m}}\sigma }$ (moving average)

polynomial degree 2 or 3: ${\displaystyle {\sqrt {\frac {3(3m^{2}-7)}{4m(m^{2}-4)}}}\sigma }$.

These functions are shown in the plot at the right. For example, with a 9-point linear function (moving average) two thirds of the noise is removed and with a 9-point quadratic/cubic smoothing function only about half the noise is removed. Most of the noise remaining is low-frequency noise(see Frequency characteristics of convolution filters, below).

Although the moving average function gives better noise reduction it is unsuitable for smoothing data which has curvature over m points. A quadratic filter function is unsuitable for getting a derivative of a data curve with an inflection point because a quadratic polynomial does not have one. The optimal choice of polynomial order and number of convolution coefficients will be a compromise between noise reduction and distortion.[28]

Frequency characteristics of convolution filters

Fourier transform of the 9-point quadratic/cubic smoothing function

Convolution maps to multiplication in the Fourier co-domain. The discrete Fourier transform of a convolution filter is a real-valued function which can be represented as

${\displaystyle FT(\theta )=\sum _{j={\tfrac {1-m}{2}}}^{\tfrac {m-1}{2}}C_{j}\cos(j\theta )}$

θ runs from 0 to 180 degrees, after which the function merely repeats itself. The plot for a 9-point quadratic/cubic smoothing function is typical. At very low angle, the plot is almost flat, meaning that low-frequency components of the data will be virtually unchanged by the smoothing operation. As the angle increases the value decreases so that higher frequency components are more and more attenuated. This shows that the convolution filter can be described as a low-pass filter: the noise that is removed is primarily high-frequency noise and low-frequency noise passes through the filter.[31] Some high-frequency noise components are attenuated more than others, as shown by undulations in the Fourier transform at large angles. This can give rise to small oscillations in the smoothed data[32] and phase reversal, i.e., high-frequency oscillations in the data get inverted by Savitzky–Golay filtering.[33]

Compared with other smoothing filters, e.g. convolution with a Gaussian or multi-pass moving-average filtering, Savitzky–Golay filters have an initially flatter response and sharper cutoff in the frequency domain, especially for high orders of the fit polynomial (see frequency characteristics). For data with limited signal bandwidth, this means that Savitzky–Golay filtering can provide better signal-to-noise ratio than many other filters; e.g., peak heights of spectra are better preserved than for other filters with similar noise suppression. Disadvantages of the Savitzky–Golay filters are comparably poor suppression of some high frequencies (poor stopband suppression) and artifacts when using polynomial fits for the first and last points.[16]

Alternative smoothing methods that share the advantages of Savitzky–Golay filters and mitigate at least some of their disadvantages are Savitzky–Golay filters with properly chosen fitting weights, Whittaker–Henderson smoothing (a method closely related to smoothing splines), and convolution with a windowed sinc function.[16]

• Kernel smoother – Different terminology for many of the same processes, used in statistics
• Local regression — the LOESS and LOWESS methods
• Numerical differentiation – Application to differentiation of functions
• Smoothing spline
• Stencil (numerical analysis) – Application to the solution of differential equations

• Gans, Peter (1992). Data fitting in the chemical sciences: By the method of least squares. ISBN 9780471934127.

Wikimedia Commons has media related to Savitzky–Golay filter.

• Advanced Convolution Coefficient Calculator (ACCC) for multidimensional least-squares filters
• Savitzky–Golay filter in Fundamentals of Statistics
• A wider range of coefficients for a range of data set sizes, orders of fit, and offsets from the centre point