Noise Reduction

Every measurement is a mixture of two things: the information we want and a layer of random variation we did not ask for. Noise reduction is the art of removing as much of the second as possible while touching the first as little as possible. It sounds simple, and the basic tools really are simple. The difficult part is knowing which unwanted variation you are dealing with, because only some of it can be smoothed away.

The language we use for noise came into analytical chemistry from electrical engineering. During the 1940s, radar and telecommunications forced engineers to think rigorously about separating meaningful signals from background interference. Claude Shannon's 1948 paper A Mathematical Theory of Communication set the theoretical limits for sending information through a noisy channel [1], and Norbert Wiener's wartime work on optimal filtering, published in 1949 as Extrapolation, Interpolation, and Smoothing of Stationary Time Series, gave the mathematics for pulling a signal out of noise [2]. Spectroscopists soon adopted the same vocabulary to describe the performance of their instruments.

Before computers, smoothing was a physical act. You slowed down the scan, opened the slit, or put an analog RC filter (a resistor and a capacitor that damp fast fluctuations) on the detector output. Each trick improved the signal-to-noise ratio, but it cost resolution or time. Once instruments began to produce digital output, a list of numbers at evenly spaced wavelengths, smoothing could be done after the measurement, in software, without changing how the data were acquired. Savitzky and Golay's 1964 paper in Analytical Chemistry was one of the first to exploit this for spectra, and it has remained among the most cited papers in the history of Analytical Chemistry [3].

In spectroscopy and chemometrics, noise is the unwanted random fluctuation superimposed on the true signal. Measure the same sample ten times under identical conditions and you will get ten slightly different spectra. The variation between those repeats is noise.

The quality of a measurement is summarised by the signal-to-noise ratio (SNR):

where $\bar{x}$ is the mean signal at a given point and $s$ is the standard deviation of repeated measurements at that point. A higher SNR means a cleaner measurement. Engineers often express it in decibels:

so an SNR of 100 corresponds to 40 dB.

Noise in a spectrometer is not one phenomenon but a family of them [8]:

• Shot noise comes from the quantum nature of light. Photons (or electrons) arrive one by one, at random moments, so the number counted in a fixed time fluctuates. These counts follow Poisson statistics, and the fluctuation grows with the square root of the signal (so the relative noise falls as the signal grows, as we will see below).
• Thermal (Johnson-Nyquist) noise comes from the random motion of electrons in the detector and electronics. It is present even when no light reaches the detector.
• Dark current is the signal a detector produces in complete darkness, because heat alone creates some charge carriers. Its average can be subtracted, but its fluctuations remain.
• Digitisation noise appears when the analog detector signal is converted into numbers. The converter can only represent a finite number of levels, so every value is rounded.

Knowing which source dominates your measurement is the first step towards choosing the right remedy.

Noise is not just an aesthetic problem. It has concrete consequences for multivariate analysis [5],[10]:

• Calibration models can fit noise as if it were signal (overfitting), which degrades predictions on new samples.
• Classification models may find spurious differences between groups that are only noise.
• Peak detection produces false positives (noise spikes mistaken for peaks) and false negatives (real peaks hidden in noise).
• Derivative spectra amplify noise strongly, because a derivative emphasises fast point-to-point changes, and noise is exactly that.
• Variable selection methods may pick noisy, uninformative variables.

The goal of noise reduction is to suppress random fluctuations while preserving the chemical information carried by the spectral features: band positions, heights, widths and shapes.

Not every unwanted wiggle is noise in the strict sense, and the distinction decides which tool will work. Play with the panel below before reading the details: pick a type of variation and switch smoothing on.

The simplest and most common model: a random value drawn from a Gaussian distribution with zero mean is added independently to each data point.

Here $x_i$ is the true signal at point $i$, ${\epsilon }_i$ is the noise, and $\sigma$ is the same everywhere in the spectrum. This is called homoscedastic noise (from the Greek for "same scatter") because its size does not depend on the signal level. A second, separate property is independence. Noise is called white when neighbouring errors are uncorrelated, so its power is spread equally over all frequencies, by analogy with white light, which contains all colours equally. Homoscedastic describes the size of the noise; white describes its independence from point to point. The model above is both.

Thermal noise is the textbook example. Most smoothing methods are designed and analysed under this model, and they work well when it holds. You met the same idea of constant error variance in the least squares foundation.

Because photons arrive randomly, the number detected in a given time follows a Poisson distribution, whose variance equals its mean. For a signal of $I$ counts:

So stronger signals carry more absolute noise, but less relative noise: the ratio $\sqrt{I}/I=1/\sqrt{I}$ falls as the signal grows. In Raman, fluorescence or low-light UV-Vis measurements, shot noise is often the dominant source. Since smoothing applies the same strength everywhere, it can be too weak on top of intense bands and stronger than needed on the baseline.

More generally, when the noise level changes across the spectrum, for instance because the detector is less sensitive near the ends of its range, the noise is called heteroscedastic ("different scatter"). This is common in real instruments. Shot noise is one particular case of it.

Some unwanted variation is not random at all:

• Baseline drift is a slow, smooth change of the whole spectrum caused by temperature, lamp ageing or detector drift. Because it is smooth, a smoother keeps it. Use baseline correction or derivatives instead.
• Interference fringes are periodic ripples caused by light bouncing between the parallel faces of a thin cell or window. They are regular waves, not random noise, and are better removed by Fourier filtering or dedicated fringe-correction methods.
• Spikes (cosmic rays) hit CCD detectors in Raman spectroscopy and produce one or two extremely intense pixels. Averaging does not delete a spike: it smears it into a fake band. Spikes should be detected and removed first, for example with a median filter or a dedicated despiking method (see median filtering below).

In diffuse reflectance measurements (NIR, DRIFTS) light scattering changes the effective path length from sample to sample. The spectrum of a sample is then shifted and stretched relative to a reference spectrum of identical chemistry:

where the offset $a$ and the multiplier $b$ change from sample to sample. These effects are smooth, so smoothing cannot touch them. They are handled by MSC (Multiplicative Scatter Correction), EMSC or SNV (Standard Normal Variate, which centres each spectrum on its own mean and divides it by its own standard deviation) [5],[7],[11].

Before any software smoothing, the most powerful noise reduction is simply to measure more than once and average. If the noise in each scan is random and independent, averaging $n$ scans divides its standard deviation by $\sqrt{n}$:

This comes straight from the statistics of the mean: the variance of an average of $n$ independent values is ${\sigma }^2/n$ [6]. Four scans halve the noise, sixteen divide it by four, a hundred divide it by ten. The returns diminish: every further halving needs four times as many scans.

Averaging has one great virtue: it introduces no distortion of the band shapes, because each scan is averaged with other measurements of the same point, not with its neighbours. This holds provided the sample and the instrument do not change between scans: lamp drift, small wavelength shifts or a sample that reacts during the measurement do blur the average. That is why FTIR instruments co-add 16, 32 or more scans by default. Its cost is time, which is not always available in kinetics or process monitoring. When it is not, we turn to smoothing.

Smoothing exploits a simple asymmetry between signal and noise. Spectroscopic signals usually change gradually from one wavelength to the next, because they come from continuous physical phenomena such as molecular vibrations or electronic transitions. Random noise changes independently at every point. If we combine each point with its neighbours, the random errors partly cancel while the slowly varying signal survives.

The catch is that real bands are not perfectly flat over a window of neighbours. The wider the neighbourhood, the more the noise cancels, but the more the band shape is bent. Every smoothing method is a different answer to this trade-off. The four methods used most in chemometrics, from simplest to most sophisticated, are summarised below. Each has its own article with derivations, interactive panels and code.

Replace each point with the plain average of itself and its neighbours:

Strengths: extremely simple and fast, one parameter (the window size $w$).

Weaknesses: all neighbours count equally, even distant ones, so bands are flattened and broadened.

Best for: quick exploration, very broad features, and learning the idea of smoothing.

Full details: Moving Average Smoothing

A weighted average in which closer neighbours count more, following the bell-shaped Gaussian curve. The parameter $\sigma$ sets the width of the bell:

Strengths: smooth results without the frequency-response ripples (sidelobes) of the moving average, explained in the Gaussian article; one parameter, extends naturally to images.

Weaknesses: still lowers and broadens sharp bands.

Best for: natural-looking smoothing, spectral images, situations where one simple parameter is welcome.

Full details: Gaussian Smoothing

Instead of averaging, fit a low-order polynomial to each window by least squares and keep its value at the centre [3]. The whole operation reduces to fixed weights:

Strengths: preserves band heights and widths much better than averaging, and computes smoothed derivatives in the same step.

Weaknesses: two parameters (window size and polynomial order), and it removes less noise than a moving average of the same length.

Best for: IR, Raman and UV-Vis spectra where band shape matters, and derivatives.

Full details: Savitzky-Golay Smoothing

Frame smoothing as an optimisation: find the curve $\mathbf{z}$ that stays close to the data while being as smooth as possible [4]. The parameter $\lambda$ sets the balance:

where ${\Delta }^d{z}_i$ is a difference of order $d$ (for $d=2$, a measure of curvature).

Strengths: one intuitive parameter, works on the whole spectrum at once, handles missing points, and is the basis of popular baseline correction methods.

Weaknesses: needs a (sparse) linear system, and $\lambda$ spans many orders of magnitude.

Best for: general-purpose smoothing and as a building block for baseline correction.

Full details: Whittaker Smoothing

The panel below applies each method to the same noisy spectrum. Pay attention to the sharp band in the middle and to the residuals underneath. For each method, look for the setting where the residuals still look like pure noise but the error against the true signal is smallest.

You will notice that no method wins on every count. The moving average and the Gaussian remove baseline noise very efficiently but quickly flatten the sharp band. Savitzky-Golay and Whittaker can be pushed further before the residuals start to show the shape of the band.

Smoothing is the most common approach in chemometrics, but there are others, used instead of it or in combination.

Any spectrum can be written as a sum of waves of different frequencies (a Fourier transform). Broad bands and baselines are made of slow, low-frequency waves; white noise is spread over all frequencies, including the fast ones. Transforming the spectrum with the Fast Fourier Transform (FFT), damping the high frequencies and transforming back reduces noise. The practical difficulty is the cut-off: too low and real features disappear, too high and noise remains. Interference fringes, which appear as a sharp peak at one frequency (in FTIR, as a small burst at one position of the interferogram), can be removed this way, something no local smoother can do.

Wavelets decompose a signal into small wave-like pieces that are localised both in position and in scale, which suits signals with sharp local features on a smooth background. Noise mostly produces small wavelet coefficients, so it can be suppressed by setting small coefficients to zero (thresholding) and reconstructing [9]. Results can be excellent, but you must choose a wavelet family, a decomposition depth and a threshold rule, which makes it harder to tune than standard smoothing.

A median filter replaces each point with the median of its window instead of the mean. The median ignores extreme values, so a single spike surrounded by normal points disappears completely, whereas a moving average only spreads it out. Median filtering is normally used to remove spikes before smoothing. It is not a good general smoother, for two reasons. It flattens curved band tops into small plateaus. And for ordinary Gaussian noise it is less efficient than the mean: the variance of a median is about $\pi /2\approx 1.57$ times that of a mean over the same window, so it removes less noise.

For Raman spikes there are also dedicated methods that only touch the suspicious points. A popular one flags a point as a spike when its jump from the previous point is far larger than the typical jump (measured with a robust "modified z-score"), and replaces only those points by the mean of their neighbours [12].

After any noise reduction, check that it did what you intended: less noise, no distortion.

Measure the noise before and after, in a region where the true signal is known to be flat (a stretch of baseline). Two precautions make the number honest: choose a region well inside the spectrum, because every smoother behaves differently at the edges, and subtract a straight line fitted to that region first, so that a sloping baseline is not counted as noise.

For white noise, a moving average of $w$ points reduces the noise standard deviation by a factor $\sqrt{w}$, so a 9-point window should improve the SNR about threefold. Much less improvement suggests that the noise is not independent from point to point, which is common when the instrument has already smoothed the data internally.

Subtract the smoothed spectrum from the original:

and plot the residuals:

• Good smoothing: residuals look like random noise, centred on zero, with no visible pattern.
• Over-smoothing: residuals show structure, such as bumps with the shape of a band. Signal is being removed.
• Under-smoothing: the smoothed curve is still visibly noisy.

Compare key features before and after smoothing:

• Band positions should not move (with centred windows they do not, for symmetric, isolated bands; an asymmetric band or two overlapping bands can show a small shift of the apparent maximum).
• Band heights will drop somewhat with moving average and Gaussian; much less with Savitzky-Golay and Whittaker.
• Band widths (full width at half maximum, FWHM) should stay approximately constant.
• Band areas are preserved even when heights change (exactly, away from the edges, for any smoother whose weights sum to one), because smoothing redistributes intensity rather than removing it.
• Resolution: closely spaced bands should remain separate. Over-smoothing merges them.

The example below builds a synthetic spectrum with known noise, smooths it with three methods and measures what each one did: the noise left on a stretch of baseline, the resulting SNR (taking the 520 nm band height, 0.7, as the signal), and how much of that band's height was lost. The noise is measured as explained above, in a band-free region inside the spectrum (700 to 780 nm) after removing a straight-line baseline. The code for each individual method is in its own article.

import numpy as np
from scipy.ndimage import uniform_filter1d, gaussian_filter1d
from scipy.signal import savgol_filter
import matplotlib.pyplot as plt

# Synthetic spectrum: two bands on a sloped baseline, plus white noise
np.random.seed(42)
wavelength = np.linspace(400, 800, 500)            # 0.8 nm between points
true_signal = (0.7 * np.exp(-((wavelength - 520)**2) / 300) +
               0.5 * np.exp(-((wavelength - 650)**2) / 200) +
               0.001 * (wavelength - 400))
noisy_spectrum = true_signal + np.random.normal(0, 0.05, wavelength.size)

# Three smoothers with comparable settings. For the moving average,
# mode='nearest' repeats the end values instead of padding with zeros.
smoothers = {
    'Moving average (w=11)': lambda y: uniform_filter1d(y, size=11, mode='nearest'),
    'Gaussian (sigma=2)': lambda y: gaussian_filter1d(y, sigma=2.0),
    'Savitzky-Golay (w=11, p=2)': lambda y: savgol_filter(y, window_length=11, polyorder=2),
}


def baseline_noise(signal, wavelength, lo=700, hi=780):
    """Noise s.d. in a band-free region inside the spectrum, after
    removing a straight-line baseline so that the slope is not counted."""
    m = (wavelength >= lo) & (wavelength <= hi)
    coef = np.polyfit(wavelength[m], signal[m], 1)
    return np.std(signal[m] - np.polyval(coef, wavelength[m]), ddof=2)


peak = np.argmin(np.abs(wavelength - 520))         # top of the 520 nm band
peak_height = 0.7
noise_raw = baseline_noise(noisy_spectrum, wavelength)
print(f"{'Raw data':27s} noise = {noise_raw:.4f}   SNR = {peak_height / noise_raw:5.1f}")

smoothed = {}
for name, smooth in smoothers.items():
    smoothed[name] = smooth(noisy_spectrum)
    noise = baseline_noise(smoothed[name], wavelength)
    # Smoothing the noise-free signal shows the distortion on its own
    height_lost = true_signal[peak] - smooth(true_signal)[peak]
    print(f"{name:27s} noise = {noise:.4f}   SNR = {peak_height / noise:5.1f}"
          f"   improvement x{noise_raw / noise:.1f}   band height lost = {height_lost:.3f}")

# Plot: smoothed spectra (top) and residuals (bottom)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
ax1.plot(wavelength, noisy_spectrum, color='gray', alpha=0.4, linewidth=0.8, label='Noisy')
for name, s in smoothed.items():
    ax1.plot(wavelength, s, linewidth=1.8, label=name)
ax1.plot(wavelength, true_signal, 'k--', linewidth=1.2, label='True signal')
ax1.axvspan(700, 780, color='gray', alpha=0.1, label='Noise region')
ax1.set_ylabel('Intensity (a.u.)')
ax1.legend(fontsize=9)

for name, s in smoothed.items():
    ax2.plot(wavelength, noisy_spectrum - s, alpha=0.7, label=name)
ax2.axhline(0, color='black', linewidth=0.5)
ax2.set_xlabel('Wavelength (nm)')
ax2.set_ylabel('Residual')
ax2.set_title('Residuals (should look like random noise)')
ax2.legend(fontsize=9)
plt.tight_layout()
plt.show()

Raw data                    noise = 0.0534   SNR =  13.1
Moving average (w=11)       noise = 0.0147   SNR =  47.6   improvement x3.6   band height lost = 0.015
Gaussian (sigma=2)          noise = 0.0193   SNR =  36.2   improvement x2.8   band height lost = 0.006
Savitzky-Golay (w=11, p=2)  noise = 0.0239   SNR =  29.3   improvement x2.2   band height lost = 0.000

% Synthetic spectrum: two bands on a sloped baseline, plus white noise
rng(42);
wavelength = linspace(400, 800, 500);              % 0.8 nm between points
true_signal = 0.7 * exp(-((wavelength - 520).^2) / 300) + ...
              0.5 * exp(-((wavelength - 650).^2) / 200) + ...
              0.001 * (wavelength - 400);
noisy_spectrum = true_signal + 0.05 * randn(size(wavelength));

% Three smoothers with comparable settings (base MATLAB, plus sgolayfilt
% from the Signal Processing Toolbox). movmean shrinks the window at the ends.
names = {'Moving average (w=11)', 'Gaussian (sigma=2)', 'Savitzky-Golay (w=11, p=2)'};
smoothers = {@(y) movmean(y, 11), @(y) gaussianSmooth(y, 2), @(y) sgolayfilt(y, 2, 11)};

peak = find(abs(wavelength - 520) == min(abs(wavelength - 520)), 1);
peak_height = 0.7;
noise_raw = baselineNoise(noisy_spectrum, wavelength);
fprintf('%-27s noise = %.4f   SNR = %5.1f\n', 'Raw data', noise_raw, peak_height / noise_raw);

smoothed = cell(size(smoothers));
for k = 1:numel(smoothers)
    smoothed{k} = smoothers{k}(noisy_spectrum);
    noise = baselineNoise(smoothed{k}, wavelength);
    clean = smoothers{k}(true_signal);              % distortion on its own
    fprintf('%-27s noise = %.4f   SNR = %5.1f   improvement x%.1f   band height lost = %.3f\n', ...
            names{k}, noise, peak_height / noise, noise_raw / noise, true_signal(peak) - clean(peak));
end

figure('Position', [100 100 1000 700]);
subplot(2, 1, 1);
plot(wavelength, noisy_spectrum, 'Color', [0.7 0.7 0.7]); hold on;
for k = 1:numel(smoothed)
    plot(wavelength, smoothed{k}, 'LineWidth', 1.8);
end
plot(wavelength, true_signal, '--k');
ylabel('Intensity (a.u.)');
legend(['Noisy', names, 'True signal']);
subplot(2, 1, 2);
hold on;
for k = 1:numel(smoothed)
    plot(wavelength, noisy_spectrum - smoothed{k});
end
yline(0, 'k');
xlabel('Wavelength (nm)'); ylabel('Residual');
title('Residuals (should look like random noise)');
legend(names);

function s = baselineNoise(y, wavelength)
    % Noise s.d. in a band-free region inside the spectrum (700-780 nm),
    % after removing a straight-line baseline so the slope is not counted
    m = wavelength >= 700 & wavelength <= 780;
    p = polyfit(wavelength(m), y(m), 1);
    r = y(m) - polyval(p, wavelength(m));
    s = sqrt(sum(r.^2) / (numel(r) - 2));
end

function ys = gaussianSmooth(y, sigma)
    % Normalised Gaussian kernel truncated at 3 sigma, renormalised at the edges
    r = ceil(3 * sigma);
    k = exp(-(-r:r).^2 / (2 * sigma^2));
    ys = conv(y, k, 'same') ./ conv(ones(size(y)), k, 'same');
end

library(signal)  # for sgolayfilt

# Synthetic spectrum: two bands on a sloped baseline, plus white noise
set.seed(42)
wavelength <- seq(400, 800, length.out = 500)      # 0.8 nm between points
true_signal <- 0.7 * exp(-((wavelength - 520)^2) / 300) +
               0.5 * exp(-((wavelength - 650)^2) / 200) +
               0.001 * (wavelength - 400)
noisy_spectrum <- true_signal + rnorm(length(wavelength), sd = 0.05)

# Moving average whose window shrinks at the ends (like MATLAB's movmean)
moving_average <- function(y, w) {
  m <- (w - 1) %/% 2
  n <- length(y)
  sapply(seq_len(n), function(i) mean(y[max(1, i - m):min(n, i + m)]))
}

# Gaussian kernel truncated at 3 sigma, renormalised at the edges
gaussian_smooth <- function(y, sigma) {
  r <- ceiling(3 * sigma)
  k <- exp(-(-r:r)^2 / (2 * sigma^2))
  n <- length(y)
  sapply(seq_len(n), function(i) {
    idx <- (i - r):(i + r)
    ok <- idx >= 1 & idx <= n
    sum(k[ok] * y[idx[ok]]) / sum(k[ok])
  })
}

smoothers <- list(
  "Moving average (w=11)"      = function(y) moving_average(y, 11),
  "Gaussian (sigma=2)"         = function(y) gaussian_smooth(y, 2),
  "Savitzky-Golay (w=11, p=2)" = function(y) as.numeric(sgolayfilt(y, p = 2, n = 11))
)

# Noise s.d. in a band-free region inside the spectrum, after removing
# a straight-line baseline so that the slope is not counted as noise
baseline_noise <- function(y, wavelength, lo = 700, hi = 780) {
  m <- wavelength >= lo & wavelength <= hi
  fit <- lm(y[m] ~ wavelength[m])
  sqrt(sum(residuals(fit)^2) / (sum(m) - 2))
}

peak <- which.min(abs(wavelength - 520))           # top of the 520 nm band
peak_height <- 0.7
noise_raw <- baseline_noise(noisy_spectrum, wavelength)
cat(sprintf("%-27s noise = %.4f   SNR = %5.1f\n", "Raw data", noise_raw, peak_height / noise_raw))

smoothed <- lapply(smoothers, function(f) f(noisy_spectrum))
for (name in names(smoothers)) {
  noise <- baseline_noise(smoothed[[name]], wavelength)
  clean <- smoothers[[name]](true_signal)          # distortion on its own
  cat(sprintf("%-27s noise = %.4f   SNR = %5.1f   improvement x%.1f   band height lost = %.3f\n",
              name, noise, peak_height / noise, noise_raw / noise,
              true_signal[peak] - clean[peak]))
}

# Plot: smoothed spectra (top) and residuals (bottom)
cols <- c("#2563eb", "#dc2626", "#059669")
par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
plot(wavelength, noisy_spectrum, type = "l", col = "gray70",
     xlab = "Wavelength (nm)", ylab = "Intensity (a.u.)")
for (k in seq_along(smoothed)) lines(wavelength, smoothed[[k]], col = cols[k], lwd = 2)
lines(wavelength, true_signal, lty = 2)
legend("topright", c("Noisy", names(smoothed), "True signal"),
       col = c("gray70", cols, "black"), lty = c(1, 1, 1, 1, 2), cex = 0.8)
residuals_all <- sapply(smoothed, function(s) noisy_spectrum - s)
matplot(wavelength, residuals_all, type = "l", lty = 1, col = cols,
        xlab = "Wavelength (nm)", ylab = "Residual",
        main = "Residuals (should look like random noise)")
abline(h = 0)

• Noise is the random part of the unwanted variation. Drift, scatter, fringes and spikes are not random and need their own tools.
• Averaging repeated scans reduces noise by $\sqrt{n}$ without any distortion. Use it first when time allows.
• Smoothing trades noise for resolution. The methods differ in how gracefully they make that trade.
• Always judge smoothing by its residuals and by what it does to band heights and widths, not by how pretty the curve looks.