Savitzky-Golay Smoothing

In 1964 two scientists at the Perkin-Elmer Corporation, one of the leading instrument makers of the time, published a paper that became one of the most cited in the history of analytical chemistry. Abraham Savitzky and Marcel J. E. Golay, a Swiss-born physicist also known for the Golay error-correcting codes, the Golay infrared detector and, for chemists, the open-tubular capillary column for gas chromatography, wrote Smoothing and Differentiation of Data by Simplified Least Squares Procedures in Analytical Chemistry [1]. The timing was not accidental. Instruments were starting to produce digital output that could be fed to the first laboratory computers, and practitioners needed practical ways to clean up noisy spectra without distorting the bands that carry the chemistry.

The paper's insight was this: instead of averaging neighbouring points, fit a polynomial through a sliding window of points by least squares, and take the value of that polynomial at the centre. Band shapes, heights and widths survive far better than with simple averaging. Even better, Savitzky and Golay showed that the whole fit collapses into a fixed set of convolution weights: integers that could be printed in a table and applied with a handful of multiplications, a decisive advantage on 1960s hardware.

Some of the published tables contained errors, which Steinier, Termonia and Deltour corrected in 1972 [2], and in 1978 Madden gave general formulas for the weights, so that nobody needs to copy them from a table [7]. Neither the errata nor the six decades since have dented the method's popularity: Savitzky-Golay filtering remains a standard tool in spectroscopy and chromatography whenever band fidelity matters [3].

The moving average and Gaussian smoothing both reduce noise by averaging, and both share the same weakness: they flatten and broaden bands. A sharp maximum is averaged with the lower points around it, so it comes out lower. Two bands that were barely resolved can merge into one. In spectroscopy, where band positions, heights and widths carry chemical information, that is a serious problem.

Savitzky-Golay smoothing keeps the idea of a sliding window but changes what we do inside it.

Fit a low-order polynomial to the points in the window, by least squares, and use its value at the centre as the smoothed value.

Take a window of 5 points around point $i$:

• A moving average adds the five values and divides by 5. If the window sits on the top of a band, the average is necessarily lower than the top.
• Savitzky-Golay fits, for example, a parabola $a_0+a_{1}t+a_2{t}^2$ through the five points and reads it at the centre ($t=0$). A parabola can bend, so on the top of a band it follows the curvature instead of cutting it off.

A walking analogy: you want your altitude on a foggy hillside, and you ask the people around you for theirs. Averaging their answers puts you too low if you stand on a summit, because everyone around you is lower. Fitting a smooth curve through their answers and reading it at your position recognises that you are on a summit.

The panel below shows the procedure. The red curve is the polynomial fitted to the points inside the shaded window; the ring marks its value at the centre. Slide the window, or let it slide, and change the window size and polynomial order.

Some things to try:

• With order 0, the "polynomial" is a constant, and the least squares constant through a set of points is their mean: Savitzky-Golay of order 0 is the moving average.
• With order 2 on the top of the sharp band, the parabola follows the curvature and the centre value stays close to the true maximum, while the moving average value (shown in the metrics) falls short.
• With a high order and a small window (for example w = 7 with order 5 or 6), the polynomial starts to follow the noise. You are fitting almost as many coefficients as there are points, and with order 6 the curve passes through all seven points and smooths nothing.

This is ordinary linear least squares, the same machinery as in the normal equations. For a window of $w=2m+1$ points we place the centre at $t=0$, so the positions are $t=-m,\dots ,m$, and we fit

The design matrix has one row per point in the window and one column per power of $t$:

and the normal equations give the coefficients:

The smoothed value is the polynomial at the centre, $\hat{y}(0)={\hat{a}}_0$, the first element of $\hat{\mathbf{a}}$.

Here comes the key observation. The matrix $({\mathbf{A}}^T\mathbf{A}{)}^{-1}{\mathbf{A}}^T$ depends only on the window size and the polynomial order, not on the data. Its first row is therefore a fixed set of numbers $c_{-m},\dots ,c_m$, and the smoothed value is just a weighted sum:

The same weights are used at every position, so the filter is a convolution, exactly like a moving average, only with cleverer weights. The weights can be computed once and reused for every point and every spectrum.

For a window of 5 points and a quadratic (order 2), the weights are

so

Notice that:

• The centre point gets the largest weight (17/35, about 0.49).
• The immediate neighbours get 12/35 (about 0.34) each.
• The outer points get a negative weight, -3/35.

The negative weights are what let the filter follow curvature. On the top of a band the outer points are low; subtracting a little of them pushes the estimate back up towards the true maximum. The weights still sum to 1, so a constant signal is left untouched.

The second panel puts the weights side by side with the flat weights of a moving average of the same length, and applies both filters to a spectrum with a narrow band.

Two numbers in the metrics deserve attention.

Band height. With the same window, the Savitzky-Golay filter keeps far more of the narrow band than the moving average: with the starting settings (w = 15, p = 2) it keeps about 85 % of the band height, against about 50 % for the moving average. A polynomial of order $p$ is reproduced exactly by the filter, so any band shape that looks locally like a parabola passes through almost untouched.

Noise. For independent random noise of standard deviation $\sigma$, the noise left after smoothing is $\sigma \sqrt{{\sum }_j{c}_j^2}$. For least squares weights this sum equals the centre weight (check it in the panel: the noise factor is the square root of c₀), so

For the 5-point quadratic, $\sqrt{17/35}\approx 0.70$, while a 5-point moving average gives $1/\sqrt{5}\approx 0.45$. For the same window, Savitzky-Golay removes less noise. That is the honest trade-off: the polynomial is flexible enough to follow bands, and the same flexibility lets some noise through. In practice you compensate by using a longer window than you would with a moving average. Because the filter keeps band shapes so much better, the longer window is usually affordable, and for a given noise reduction Savitzky-Golay distorts sharp bands much less [4],[5].

The same trade-off can be read in terms of frequencies. Think of the spectrum as a sum of waves: broad bands are slow waves, noise is mostly fast waves. A Savitzky-Golay filter is a low-pass filter, one that lets slow waves through and damps fast ones. Compared with a moving average of the same length, it keeps a wider range of slow waves almost untouched (a flatter pass band) and moves the cut-off to faster waves, which is why it removes less noise. Like the moving average, it lets part of the fast waves leak through in small ripples called side lobes [4].

Savitzky-Golay smoothing has two parameters.

As with the moving average, the window decides how many points take part:

• Small windows (5 to 7 points): light smoothing, fine details kept.
• Medium windows (9 to 21 points): moderate smoothing, typical for spectra.
• Large windows (25 points and more): heavy smoothing for very noisy data or densely sampled spectra.

What really matters is the window compared with the width of your narrowest band. A quadratic filter tolerates windows up to roughly the full width at half maximum (FWHM) of a band before the band height drops noticeably; a moving average starts to hurt well before that. For a Gaussian band whose FWHM equals the window, a quadratic filter keeps about 98 % of the height and a moving average about 80 %; with a window of 1.5 times the FWHM the figures are about 93 % and 64 %.

The order sets the flexibility of the local fit:

• Order 0: a constant, which is the moving average.
• Order 2 (equivalently 3, for smoothing): the usual choice; follows the curvature of band tops.
• Order 4 (equivalently 5): more flexible, keeps sharper bands with a given window, but removes less noise, so it needs a longer window.
• Higher orders: rarely useful for smoothing.

A polynomial of order $p$ has $p+1$ coefficients, so the window must contain more than $p+1$ points (if $w=p+1$, the polynomial passes through every point and nothing is smoothed). A comfortable rule is $w\ge 2p+3$ (for example $w\ge 7$ for $p=2$), and preferably more.

The first and last $m$ points do not have a full centred window. A natural solution, used by scipy.signal.savgol_filter (default mode='interp'), MATLAB's sgolayfilt and R's signal::sgolayfilt, is to fit the polynomial to the first (or last) complete window and evaluate it at the edge positions instead of at the centre. Those edge estimates are noisier than the interior ones, because a polynomial is less well determined near the end of its fitting range. This edge treatment, with its own tables of weights for each edge position, was worked out by Gorry [6]. Other implementations pad the data by mirroring or repeating values. Whatever the method, treat the first and last half-window of points with some caution.

Because Savitzky-Golay fits a polynomial, it can also return the derivative of that polynomial at the centre. The first derivative at $t=0$ is ${\hat{a}}_1$, the second is $2{\hat{a}}_2$, and so on. Each is again a fixed set of weights, so smoothed derivatives cost no more than smoothing [1]. The same odd-even pairing as for smoothing appears here, shifted by one: for the first derivative, order 2 gives exactly the same weights as order 1. For example, the 5-point quadratic first-derivative weights are

(per data point; divide by the spacing $\Delta x$ to get the derivative per wavelength unit). Differentiating raw noisy data amplifies the noise badly; differentiating a local least squares fit does not. This is why Savitzky-Golay is the standard way to compute spectral derivatives in chemometrics. It is covered in detail in Spectral Derivatives.

import numpy as np
from scipy.signal import savgol_filter, savgol_coeffs
import matplotlib.pyplot as plt

# Noisy synthetic spectrum
rng = np.random.default_rng(42)
x = np.linspace(0, 100, 500)
dx = x[1] - x[0]
true_signal = (0.8 * np.exp(-((x - 25)**2) / 20) +
               0.6 * np.exp(-((x - 50)**2) / 15) +
               0.9 * np.exp(-((x - 75)**2) / 25))
noisy_signal = true_signal + rng.normal(0, 0.08, len(x))

# The weights themselves: window 5, order 2  ->  [-3, 12, 17, 12, -3] / 35
print(savgol_coeffs(5, 2) * 35)

# savgol_filter(data, window_length, polyorder)
settings = [(5, 2), (11, 2), (11, 4), (21, 2)]
fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True, sharey=True)
for ax, (w, p) in zip(axes.flat, settings):
    smoothed = savgol_filter(noisy_signal, window_length=w, polyorder=p)
    ax.plot(x, noisy_signal, color='gray', alpha=0.5, label='Noisy')
    ax.plot(x, smoothed, linewidth=2, label=f'w = {w}, p = {p}')
    ax.plot(x, true_signal, 'k--', alpha=0.7, label='True')
    ax.set_title(f'Savitzky-Golay, w = {w}, p = {p}')
    ax.legend()
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Smoothing and derivatives in one call.
# delta = sample spacing, so the derivatives are per x unit, not per point.
smooth = savgol_filter(noisy_signal, 11, 2)
first_derivative = savgol_filter(noisy_signal, 11, 2, deriv=1, delta=dx)
second_derivative = savgol_filter(noisy_signal, 11, 2, deriv=2, delta=dx)

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
axes[0].plot(x, noisy_signal, color='gray', alpha=0.5)
axes[0].plot(x, smooth, linewidth=2)
axes[0].set_ylabel('Signal')
axes[1].plot(x, first_derivative, color='red')
axes[1].set_ylabel('1st derivative')
axes[2].plot(x, second_derivative, color='green')
axes[2].set_ylabel('2nd derivative')
axes[2].set_xlabel('x')
plt.tight_layout()
plt.show()

% Savitzky-Golay smoothing with the Signal Processing Toolbox
rng(42);
x = linspace(0, 100, 500);
dx = x(2) - x(1);
true_signal = 0.8 * exp(-((x - 25).^2) / 20) + ...
              0.6 * exp(-((x - 50).^2) / 15) + ...
              0.9 * exp(-((x - 75).^2) / 25);
noisy_signal = true_signal + 0.08 * randn(size(x));

% sgolayfilt(data, polynomial_order, window_length)  <- note the order!
settings = [5 2; 11 2; 11 4; 21 2];
figure('Position', [100 100 1200 800]);
for k = 1:size(settings, 1)
    w = settings(k, 1); p = settings(k, 2);
    subplot(2, 2, k);
    plot(x, noisy_signal, 'Color', [0.6 0.6 0.6]); hold on;
    plot(x, sgolayfilt(noisy_signal, p, w), 'LineWidth', 2);
    plot(x, true_signal, '--k', 'LineWidth', 1.5);
    title(sprintf('Savitzky-Golay, w = %d, p = %d', w, p));
    grid on;
end

% Derivatives: sgolay returns the differentiation filters in g.
% Column d+1 of g estimates the d-th derivative (per point) divided by d!.
% conv flips the filter, hence the (-1)^d factor, as in the MATLAB documentation.
order = 2;
framelen = 11;
[~, g] = sgolay(order, framelen);
smoothed     = conv(noisy_signal,  g(:, 1), 'same');
first_deriv  = conv(noisy_signal, -g(:, 2), 'same') / dx;          % 1!/(-dx)^1
second_deriv = conv(noisy_signal, 2 * g(:, 3), 'same') / dx^2;     % 2!/(-dx)^2
% conv(...,'same') pads with zeros, so ignore the first and last (framelen-1)/2 points.

figure;
subplot(3, 1, 1); plot(x, noisy_signal, 'Color', [0.6 0.6 0.6]); hold on;
plot(x, smoothed, 'LineWidth', 2); ylabel('Signal'); grid on;
subplot(3, 1, 2); plot(x, first_deriv, 'r', 'LineWidth', 1.5); ylabel('1st derivative'); grid on;
subplot(3, 1, 3); plot(x, second_deriv, 'g', 'LineWidth', 1.5); ylabel('2nd derivative');
xlabel('x'); grid on;

library(signal)  # sgolay, sgolayfilt

set.seed(42)
x <- seq(0, 100, length.out = 500)
dx <- x[2] - x[1]
true_signal <- 0.8 * exp(-((x - 25)^2) / 20) +
               0.6 * exp(-((x - 50)^2) / 15) +
               0.9 * exp(-((x - 75)^2) / 25)
noisy_signal <- true_signal + rnorm(length(x), sd = 0.08)

# sgolayfilt(x, p = polynomial order, n = window length)
settings <- list(c(5, 2), c(11, 2), c(11, 4), c(21, 2))
par(mfrow = c(2, 2))
for (s in settings) {
  plot(x, noisy_signal, type = "l", col = "gray",
       main = sprintf("Savitzky-Golay, w = %d, p = %d", s[1], s[2]),
       xlab = "x", ylab = "Signal")
  lines(x, sgolayfilt(noisy_signal, p = s[2], n = s[1]), col = "blue", lwd = 2)
  lines(x, true_signal, lty = 2, lwd = 1.5)
}

# Derivatives: m = derivative order, ts = sample spacing
first_deriv  <- sgolayfilt(noisy_signal, p = 2, n = 11, m = 1, ts = dx)
second_deriv <- sgolayfilt(noisy_signal, p = 2, n = 11, m = 2, ts = dx)

par(mfrow = c(3, 1))
plot(x, sgolayfilt(noisy_signal, p = 2, n = 11), type = "l", col = "blue", ylab = "Signal")
plot(x, first_deriv, type = "l", col = "red", ylab = "1st derivative")
plot(x, second_deriv, type = "l", col = "darkgreen", ylab = "2nd derivative", xlab = "x")

It is the default choice when:

• Your spectrum has sharp bands whose heights and widths matter.
• You need derivatives (it is the standard method).
• You want a well-known, reproducible method that reviewers recognise.
• You work with IR, Raman, NIR or UV-Vis data.

Consider something else when:

• The signal is not smooth at all (spikes, steps): remove spikes first, for example with a median filter.
• You need to handle missing points or want a single global smoothness parameter: try Whittaker smoothing.
• You only need a quick look: a moving average is enough.

Start with p = 2. It works for most spectra. Remember that 3 gives the same smoothing.

Move to p = 4 if sharp bands are still losing height at the window you need for noise reduction, and lengthen the window accordingly.

Start with a window somewhat smaller than the FWHM of your narrowest important band (in points), then adjust:

• More noise: increase the window.
• Bands losing height: decrease the window, or raise the order.
• Always: keep $w$ comfortably larger than $p+1$.

1. Visual check: is the noise reduced without bands being flattened?
2. Residuals: original minus smoothed should look like pure noise.
3. Band metrics: compare heights and widths before and after.
4. Model performance: if the data feed a calibration model, does the cross-validated error improve?

• Raising the order without lengthening the window: less smoothing, more noise.
• Choosing p = 3 to "capture asymmetry": for smoothing it is identical to p = 2.
• A window much wider than the bands: even Savitzky-Golay will flatten them.
• Trusting the edges: the first and last half-window are noisier.
• Swapping arguments: MATLAB's sgolayfilt takes the order before the window length; SciPy takes the window first.
• Different settings for calibration and test data: always preprocess consistently.

• Savitzky-Golay smoothing fits a polynomial to each window by least squares and keeps its centre value.
• Because the fit is linear in the data, it reduces to fixed convolution weights, some of them negative.
• It keeps band heights and widths much better than averaging, at the cost of removing less noise for the same window.
• Order 0 is the moving average; orders 2 and 3 smooth identically.
• The same fit gives smoothed derivatives, which is why it is the standard tool for spectral derivatives.

• Learn about Whittaker smoothing, which treats the whole spectrum as one optimisation problem.
• Explore spectral derivatives to remove baselines and resolve overlapping bands.
• Compare all methods in the noise reduction overview.