Gaussian Smoothing

Few mathematical objects carry a person's name as firmly as the Gaussian function carries that of Carl Friedrich Gauss (1777 to 1855). Yet the bell-shaped curve was not his invention. In 1733 the French-born mathematician Abraham de Moivre, working in London, derived it as an approximation to the binomial distribution, and later included the result in the second edition of his Doctrine of Chances [1]. It is a good example of what the statistician Stephen Stigler called the law of eponymy: scientific discoveries are rarely named after the person who made them first [7].

What Gauss gave the curve was its central role in the theory of measurement errors. In 1801 the newly discovered asteroid Ceres was lost in the glare of the Sun after only a few weeks of observations. Gauss computed its orbit from those scattered measurements, and astronomers found Ceres again at the end of 1801, close to where he had predicted. In Theoria Motus Corporum Coelestium (1809) he gave his own derivation of least squares (Legendre had published the method in 1805 [8]) and justified it by assuming that measurement errors follow exactly this bell-shaped law [2]. That is why we call it the Gaussian, or normal, distribution. (The story of Ceres and least squares is told in more detail in Least Squares: optimization.)

The function itself, $e^{-x^2/(2{\sigma }^2)}$, turned out to have remarkable properties. It is infinitely smooth, its Fourier transform is again a Gaussian, and in two dimensions it splits into a product of two one-dimensional Gaussians. These properties became important when digital signal processing needed smoothing kernels: sets of weights used to average out noise. Unlike a flat moving average window, a Gaussian kernel gives the most weight to the centre and fades out smoothly.

Weighted kernels were formalised in statistics by Murray Rosenblatt (1956) and Emanuel Parzen (1962) for estimating probability densities, where the Gaussian kernel became a standard choice [3],[4]. From there the idea spread to image processing, where "Gaussian blur" is a basic operation, and to spectroscopy, where it offers a simple one-parameter alternative to the moving average.

The moving average gives every point in its window the same weight. That raises an obvious question:

Why should a point five steps away count as much as the immediate neighbour?

If the signal changes gradually, the closest neighbours are the most similar to the point we are estimating, so they deserve more influence. Points further away are still useful for cancelling noise, but they are more likely to sit on a different part of a band. This is the motivation for weighted averaging.

Gaussian smoothing uses the bell curve to set the weights: the centre point receives the highest weight, and the weight falls off smoothly with distance.

The weight given to a neighbour at distance $u$ (measured in data points) is

where $\sigma$ (sigma) controls the width of the bell. The smoothed value at point $i$ is then

In plain words: multiply each neighbour by its weight, add everything up, and divide by the sum of the weights. Dividing by the sum normalises the kernel, so that a constant signal stays exactly the same after smoothing. You may have seen the Gaussian written with a factor $1/(\sigma \sqrt{2\pi })$ in front; that factor makes the continuous curve integrate to one, but after the normalisation above it cancels and does not matter.

In practice the weights are cut off where they become negligible. About 99.7 % of the area of a Gaussian lies within $\pm 3\sigma$ of the centre, so most implementations truncate the kernel at a radius of three or four $\sigma$. For $\sigma =2$ points, truncated at $\pm 3\sigma$, the 13 normalised weights are:

offset:  -6     -5     -4     -3     -2     -1      0     +1     +2     +3     +4     +5     +6
weight: 0.002  0.009  0.027  0.065  0.121  0.176  0.200  0.176  0.121  0.065  0.027  0.009  0.002

Compare with the 7-point moving average that removes the same amount of white noise (see below), where every weight is $1/7\approx 0.143$: the Gaussian puts more weight on the centre and spreads a little weight much further out.

The parameter $\sigma$ is the standard deviation of the bell, expressed in data points:

• Small σ (about 0.5 to 1): a narrow bell, only immediate neighbours matter, light smoothing. At σ = 0.5 the centre point keeps 79 % of the weight, so the filter barely smooths at all.
• Medium σ (about 1.5 to 3): neighbours within a few points contribute, moderate smoothing.
• Large σ (4 and above): many neighbours contribute, heavy smoothing, and bands start to blur.

Use the panel below to see the kernel and its effect on different "frequencies". The lower plot needs a word of explanation: any spectrum can be described as a sum of waves. Broad bands are made of slow waves (on the left of the plot), while point-to-point noise is made of fast waves (on the right). The curve shows what fraction of each wave survives the smoothing. An ideal smoother would keep the left side at 1 and push the right side to 0.

Two things stand out:

1. The Gaussian's response falls smoothly and never goes negative. Fast waves are progressively removed. The Fourier transform of a Gaussian is again a Gaussian, which is exactly what you see.
2. The moving average's response has ripples, called sidelobes, that dip below zero. Some fast components leak through the filter, and those in the negative lobes even come out with their sign flipped. This is the origin of the slightly ragged look that moving averages sometimes leave behind [5].

For independent random noise with standard deviation ${\sigma }_{\epsilon }$, any normalised set of weights $w_j$ leaves a noise standard deviation of

For the moving average, ${\sum }_j{w}_j^2=1/w$, which gives the familiar $\sqrt{w}$ law. For a Gaussian kernel that is not too narrow, ${\sum }_j{w}_j^2\approx 1/(2\sqrt{\pi }\sigma )$. The two filters therefore remove the same amount of white noise when

So a Gaussian with $\sigma =2$ removes about as much noise as a 7-point moving average, not an 11- or 13-point one, even though its weights stretch over 13 points. This is a useful rule when you translate settings from one method to the other.

• Points beyond $3\sigma$ contribute almost nothing.
• The kernel actually used therefore spans about $6\sigma +1$ points. Libraries differ: SciPy's gaussian_filter1d truncates at $4\sigma$ by default ($8\sigma +1$ points), and MATLAB's imgaussfilt at about $2\sigma$ (about $4\sigma +1$ points).
• Example: $\sigma =2$ uses about 13 points; $\sigma =3$ about 19 points.

When reading software documentation, check whether a function asks for $\sigma$ or for a window length, and how the two are related. They are not the same thing.

Is the Gaussian simply better? It is better behaved, but it is not magic. The panel below compares the two filters at equal noise reduction, using the rule above.

What you should see:

• The noise levels on the flat baseline are practically identical, as the theory predicts.
• The band heights are close. At equal noise reduction the two kernels have a similar spread, so they broaden bands by a similar amount. For moderate and large σ (about 2 and above) the Gaussian keeps a little more of the narrowest bands, because its weight is concentrated near the centre instead of being spread evenly to the edges of a box. At the panel's starting value, σ = 3, the narrow band at point 120 keeps about 77 % of its height with the Gaussian and 71 % with the moving average; for small σ the difference disappears.
• The shapes differ. Switch on the spike: the moving average turns it into a flat-topped rectangle, the Gaussian into a smooth bell. Each filter leaves its own kernel shape behind wherever the signal changes abruptly.

So the honest summary is: Gaussian smoothing produces cleaner, more natural curves without ripples, and a modest advantage on narrow bands, but it is still an averaging filter: sharp bands are lowered and broadened by roughly the same mechanism as with a moving average. If band heights matter, you need a different idea, such as the local polynomial fits of Savitzky-Golay or the penalised fit of Whittaker.

Near the ends of a spectrum the kernel runs out of points on one side. Common strategies are to renormalise (use only the available points and divide by the sum of their weights, as the panels above do), to reflect the spectrum at the ends (the default in scipy.ndimage.gaussian_filter1d), or to repeat the last value. Each gives slightly different results in the first and last $3\sigma$ points.

As with the moving average, the right σ depends less on the noise than on the narrowest band that matters, measured in data points:

These numbers are only starting points. Two anchors behind them:

• Relate σ to your narrowest band. A Gaussian band with standard deviation $s$ (in points) becomes, after smoothing, a band with standard deviation $\sqrt{s^2+{\sigma }^2}$. Keeping $\sigma$ below about a third of $s$, that is below about FWHM/7 (because the FWHM of a Gaussian is about $2.355s$), limits the broadening to roughly 5 %.
• Check the residuals. The difference between original and smoothed data should look like pure noise.

1. Visual check: smooth, but no band merged or flattened?
2. Residuals: original minus smoothed should look like noise.
3. Feature preservation: are band positions, heights and widths still where they should be?

import numpy as np
from scipy.ndimage import gaussian_filter1d
import matplotlib.pyplot as plt

# Noisy synthetic spectrum
rng = np.random.default_rng(42)
x = np.linspace(0, 100, 500)
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))

# gaussian_filter1d(data, sigma): sigma is given in data points.
# By default the kernel is truncated at 4 sigma and the edges are reflected.
fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True, sharey=True)
for ax, sigma in zip(axes.flat, [1.0, 2.0, 3.0, 5.0]):
    smoothed = gaussian_filter1d(noisy_signal, sigma=sigma)
    ax.plot(x, noisy_signal, color='gray', alpha=0.5, label='Noisy')
    ax.plot(x, smoothed, linewidth=2, label=f'σ = {sigma}')
    ax.plot(x, true_signal, 'k--', alpha=0.7, label='True')
    ax.set_title(f'Gaussian smoothing, σ = {sigma} points')
    ax.legend()
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()


def gaussian_kernel(sigma, radius=None):
    """Normalised Gaussian weights from -radius to +radius."""
    if radius is None:
        radius = int(np.ceil(3 * sigma))
    u = np.arange(-radius, radius + 1)
    k = np.exp(-u**2 / (2 * sigma**2))
    return u, k / k.sum()


# Kernels for several sigmas
for sigma in [0.5, 1.0, 2.0, 3.0]:
    u, k = gaussian_kernel(sigma, radius=10)
    plt.plot(u, k, marker='o', label=f'σ = {sigma}')
plt.xlabel('Distance from centre (points)')
plt.ylabel('Weight')
plt.legend()
plt.show()

# Moving average with the same white-noise reduction: w ≈ 2·sqrt(pi)·sigma
sigma = 2.0
w = int(2 * round((2 * np.sqrt(np.pi) * sigma - 1) / 2) + 1)   # nearest odd number, here 7
moving_avg = np.convolve(noisy_signal, np.ones(w) / w, mode='same')
gauss = gaussian_filter1d(noisy_signal, sigma=sigma)

plt.plot(x, noisy_signal, color='gray', alpha=0.4, label='Noisy')
plt.plot(x, moving_avg, linewidth=2, label=f'Moving average (w = {w})')
plt.plot(x, gauss, '--', linewidth=2, label=f'Gaussian (σ = {sigma})')
plt.plot(x, true_signal, 'k:', label='True signal')
plt.legend()
plt.show()

% Gaussian smoothing with an explicit, normalised kernel (base MATLAB)
rng(42);
x = linspace(0, 100, 500);
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));

sigmas = [1, 2, 3, 5];
figure('Position', [100 100 1200 800]);
for k = 1:numel(sigmas)
    subplot(2, 2, k);
    plot(x, noisy_signal, 'Color', [0.6 0.6 0.6]); hold on;
    plot(x, gaussianSmooth(noisy_signal, sigmas(k)), 'LineWidth', 2);
    plot(x, true_signal, '--k', 'LineWidth', 1.5);
    title(sprintf('Gaussian smoothing, \\sigma = %g points', sigmas(k)));
    legend('Noisy', 'Smoothed', 'True', 'Location', 'best');
    grid on;
end

% Moving average with the same white-noise reduction (w ≈ 3.5 sigma)
sigma = 2;
w = 2 * round((2 * sqrt(pi) * sigma - 1) / 2) + 1;       % here 7
moving_avg = movmean(noisy_signal, w);                   % shrinks window at the edges
figure;
plot(x, noisy_signal, 'Color', [0.7 0.7 0.7]); hold on;
plot(x, moving_avg, 'LineWidth', 2);
plot(x, gaussianSmooth(noisy_signal, sigma), '--', 'LineWidth', 2);
plot(x, true_signal, ':k', 'LineWidth', 1.5);
legend('Noisy', sprintf('Moving average (w = %d)', w), 'Gaussian (\sigma = 2)', 'True');
grid on;

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

# Normalised Gaussian kernel truncated at 3 sigma
gaussian_kernel <- function(sigma) {
  r <- ceiling(3 * sigma)
  u <- -r:r
  k <- exp(-u^2 / (2 * sigma^2))
  k / sum(k)
}

# Gaussian smoothing with renormalisation at the edges
gaussian_smooth <- function(y, sigma) {
  k <- gaussian_kernel(sigma)
  num <- stats::filter(c(rep(0, length(k)), y, rep(0, length(k))), k, sides = 2)
  den <- stats::filter(c(rep(0, length(k)), rep(1, length(y)), rep(0, length(k))), k, sides = 2)
  idx <- (length(k) + 1):(length(k) + length(y))
  as.vector(num[idx] / den[idx])
}

set.seed(42)
x <- seq(0, 100, length.out = 500)
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)

par(mfrow = c(2, 2))
for (s in c(1, 2, 3, 5)) {
  plot(x, noisy_signal, type = "l", col = "gray",
       main = paste0("Gaussian smoothing, sigma = ", s), xlab = "x", ylab = "Signal")
  lines(x, gaussian_smooth(noisy_signal, s), col = "blue", lwd = 2)
  lines(x, true_signal, lty = 2, lwd = 1.5)
}

# Moving average with the same white-noise reduction (w ~ 3.5 sigma)
sigma <- 2
w <- 2 * round((2 * sqrt(pi) * sigma - 1) / 2) + 1   # here 7
moving_avg <- stats::filter(noisy_signal, rep(1 / w, w), sides = 2)
# stats::filter leaves NA in the first and last (w - 1) / 2 points;
# lines() below simply leaves those points out

par(mfrow = c(1, 1))
plot(x, noisy_signal, type = "l", col = "gray80", xlab = "x", ylab = "Signal",
     main = "Gaussian vs moving average at equal noise reduction")
lines(x, moving_avg, col = "orange", lwd = 2)
lines(x, gaussian_smooth(noisy_signal, sigma), col = "blue", lwd = 2, lty = 2)
lines(x, true_signal, lwd = 1.5, lty = 3)

Gaussian smoothing extends naturally to images and spectral images. In two dimensions the weights are

The second form shows that the 2D kernel is a product of two 1D kernels. You can therefore smooth all rows with a 1D Gaussian and then all columns, which is much cheaper than a full 2D convolution. The result looks the same in every direction (the kernel is isotropic), which matters for images. This is why Gaussian smoothing is common in image processing [6] and in hyperspectral imaging, where each wavelength channel can be smoothed spatially as an image, so that neighbouring pixels share information without mixing different wavelengths.

• Confusing σ with the window size. The kernel spans about $6\sigma +1$ points.
• Using σ too large. Bands merge and flatten.
• Forgetting to normalise the kernel. The output changes scale.
• Ignoring edge conventions. Different libraries treat the ends differently.
• Expecting a big gain in band preservation over a moving average at the same noise reduction. The main gain is in smoothness, not in resolution.

A few properties explain why the Gaussian keeps appearing in smoothing:

1. It is localised in both domains at once. Among all functions, the Gaussian has the smallest possible product of its width in position and its width in frequency. This is the same mathematics behind the uncertainty principle in quantum mechanics.
2. Its Fourier transform is again a Gaussian, so its frequency response decays smoothly without ripples.
3. It is separable in several dimensions.
4. Repeated averaging leads to it. Smoothing several times with any small positive kernel (for example a moving average) produces weights that approach a Gaussian shape, a consequence of the central limit theorem.

• Gaussian smoothing is a weighted moving average whose weights follow a bell curve of width $\sigma$.
• It reduces white noise like a moving average of about $3.5\sigma$ points, and it does so without ripples in its frequency response.
• It gives smoother, more natural-looking curves, but at the same noise reduction it preserves sharp bands only slightly better than a moving average.
• Choose $\sigma$ relative to the width of your narrowest band, and check the residuals.

• Learn how Savitzky-Golay smoothing preserves bands by fitting local polynomials.
• See how Whittaker smoothing turns smoothing into an optimisation problem.
• Explore derivatives for resolution enhancement and baseline removal.
• Read about baseline correction for removing slow drifts.