Moving Average Smoothing

The idea of smoothing data by averaging neighbouring values is older than its name. In 1901 the English statistician R. H. Hooker published a study relating marriage rates to trade prosperity in which he computed what he called "instantaneous averages": each year's figure was replaced by the mean of itself and the years around it, so that the long-term trend could be seen through the year-to-year fluctuations [7]. In 1909 G. U. Yule described Hooker's averages as "moving-averages" in the Journal of the Royal Statistical Society [8], and the term spread after W. I. King used it in his textbook The Elements of Statistical Method (1912) [9].

The practice itself is older still. In 1870 the actuary W. S. B. Woolhouse published a method for "graduating" mortality tables that replaced each death rate by a weighted average of the rates at neighbouring ages, to iron out irregularities caused by small numbers [11]. In 1884 the physicist J. H. Poynting averaged successive years of wheat prices and of cotton and silk imports to see the trend beneath the year-to-year fluctuations [10]. Both used the procedure long before it had a name.

What makes the moving average so durable is its simplicity. When digital signal processing appeared in the mid-20th century and engineers needed filters that the earliest computers could run, the moving average (called a boxcar filter in that world, because its weights form a flat rectangular box) was the obvious first choice. It needs only additions and one division. That is why it is the first smoothing method most scientists learn, and the natural starting point before more refined methods such as Savitzky-Golay or Whittaker smoothing.

Real spectroscopic data are noisy. When you record a fluorescence emission spectrum, an IR absorption scan or a Raman spectrum, the result is not a smooth curve. It wiggles up and down, even when the same sample is measured twice under identical conditions.

This is not a flaw of your instrument. It is physics. A spectrometer counts photons (or electrons), their arrival is random, and the detector and electronics add their own thermal noise, so even the best instruments produce noisy data. We discussed where measurement uncertainty comes from in Least Squares: your first step into chemometrics, and the different kinds of noise in the noise reduction overview.

Your brain sees through the noise and recognises a band here and a baseline there. A computer only sees numbers: 0.421, 0.398, 0.447, 0.389... Random fluctuations can create false peaks or hide real ones, they add uncertainty to band heights and areas, and they can be fitted by calibration models as if they were information. To reduce these problems we use smoothing: methods that remove random fluctuations while keeping the underlying chemical information.

The moving average rests on one observation:

Noise is independent from point to point. Signal, very often, is continuous.

Many spectroscopic signals come from physical phenomena (molecular vibrations, electronic transitions, fluorescence emission) that tend to change gradually with wavelength. A molecule's absorbance does not jump wildly from one nanometre to the next. If the absorbance at 500 nm is 0.40, the absorbance at 501 nm will be close to 0.40, not 0.95 or 0.02. Neighbouring points carry almost the same information.

Noise does not behave like that. If it pushes the value at 500 nm up, it is just as likely to push the value at 501 nm down.

Suppose the true signal at 499, 500 and 501 nm is 0.39, 0.40 and 0.41, but noise makes you measure 0.35, 0.46 and 0.38. The average of the three measurements is 0.397, very close to the true value at 500 nm. The errors (some positive, some negative) partly cancelled; the signal, nearly the same at the three points, survived.

Not every signal is continuous, though. Mass spectra or atomic emission lines can be so sharp that neighbouring points are not similar, and averaging them destroys the feature. Keep the nature of your signal in mind before smoothing anything.

Moving average smoothing (also called running mean or boxcar filter) applies this principle by replacing each point with the average of itself and its neighbours:

where:

• ${\hat{y}}_i$ is the smoothed value at point $i$,
• $w$ is the window size (the number of points averaged),
• $m$ is the number of neighbours taken on each side.

That is the whole method.

Before going further, watch the window slide along a spectrum:

Try this:

1. Click Slide the window to watch the window travel across the spectrum, building the smoothed curve point by point.
2. Use the second slider to step through the points by hand.
3. Change the window size and observe the smoothing effect.
4. Notice how larger windows give smoother curves but lower and widen the bands.
5. Move to the first or last points and see the window shrink, because there are no neighbours beyond the edge.

The ringed red points are the neighbours being averaged, and the orange ring is the smoothed value they produce.

For a spectrum with intensities $y_1,y_2,\dots ,y_n$:

Suppose you measured these intensities:

Original: [10.2, 15.1, 12.3, 19.8, 18.2, 15.9, 14.1, 18.7, 17.2, 15.3]

With a window of 3, each interior point is averaged with its two immediate neighbours (counting positions from 0):

• Position 1: $(10.2+15.1+12.3)/3=12.5$
• Position 2: $(15.1+12.3+19.8)/3=15.7$
• Position 3: $(12.3+19.8+18.2)/3=16.8$
• Position 4: $(19.8+18.2+15.9)/3=18.0$

Continuing to the end gives:

Smoothed (positions 1 to 8): [12.5, 15.7, 16.8, 18.0, 16.1, 16.2, 16.7, 17.1]

The jump from 12.3 to 19.8 has become a gentler climb. The first and last positions have no neighbour on one side; we come back to them in the section on edge effects.

We can say exactly how much, for random noise. Assume each point carries independent noise with standard deviation $\sigma$. The smoothed value is the sum of $w$ noisy values divided by $w$. Variances of independent errors add, so

and the noise standard deviation after smoothing is

A 9-point window divides the noise by 3, a 25-point window by 5 [4],[6]. This is the same square-root law that governs averaging repeated scans, with one crucial difference: repeated scans average measurements of the same point, while the moving average mixes different points. As long as the signal is nearly constant across the window that does not matter. When it is not, the signal is distorted too.

There is a stronger result hiding here. Among all averaging filters of the same length whose weights sum to one, equal weights give the smallest output variance, ${\sigma }^2/w$: the moving average is the best possible noise remover for a given window [12]. (Minimising ${\sum }_j{w}_j^2$ with ${\sum }_j{w}_j=1$ gives $w_j=1/w$ for every point.) Its weakness is not noise, but what it does to the signal.

That distortion is the price of the noise reduction. Averaging a sharp band with the baseline around it pulls its top down and spreads it sideways. The panel below makes the trade-off measurable. The metrics compare the noise left after smoothing with the $\sigma /\sqrt{w}$ prediction, and show how much of each band's height survives. Because the panel knows the true signal and the noise it added, it can measure both without the noise getting in the way.

Notice what happens as you widen the window:

• The noise follows the $1/\sqrt{w}$ law closely.
• The broad band barely changes until the window becomes large.
• The narrow band starts losing height as soon as the window approaches its width at half height (the FWHM, full width at half maximum). When the window is about as wide as the band, roughly a fifth of its height is gone.

So the right window depends less on the noise than on your narrowest feature that matters.

The same chemistry can be sampled with 100 points or with 10,000, so window sizes only make sense relative to band width:

• A spectrum with 100 points and bands about 10 to 20 points wide: start with $w=3$ or 5; very noisy data may tolerate 7.
• A high-resolution spectrum with bands about 100 points wide: start around $w=11$ to 15, and go up to 31 or more if the noise demands it.
• Very dense data (tens of thousands of points): windows of 101 or 201 points can be perfectly reasonable. The ratio to the band width is what counts.

With an odd window (3, 5, 7, ...) there are equally many neighbours on each side of the central point:

• $w=5$: 2 left, the centre, 2 right
• $w=7$: 3 left, the centre, 3 right

With an even window (4, 8, ...) the window cannot be centred on a point, so the result belongs halfway between two points. If it is stored at one of them, the whole smoothed spectrum shifts by half a point. Odd windows keep everything centred and band positions unchanged.

At the start and end of the spectrum there are not enough neighbours on one side to fill the window. Software handles this in different ways, and it is worth knowing which one yours uses:

• Leave the edges out: only the interior points with complete windows are smoothed, and the first and last $m$ points are either kept unsmoothed or marked as missing. R's stats::filter returns NA there, for example.
• Zero padding: the spectrum is extended with zeros. This is what np.convolve(..., mode='same') in NumPy and conv(..., 'same') in MATLAB do, and it pulls the edge values towards zero, which is usually wrong for spectra.
• Edge padding: the first and last values are repeated to create artificial neighbours.
• Reflection: the spectrum is mirrored at the edges.
• Shrinking window: near the edges only the neighbours that exist are used. MATLAB's movmean and smoothdata do this by default with a truncated, one-sided window (at the very first point of a 5-point window they average points 1 to 3). The panels on this page use a smaller symmetric window instead, so the very first and last points are left unsmoothed.

The interesting features are often in the middle of the spectrum, so edges are rarely critical. But if you integrate up to the end of the range, or the spectrum is short, check what your software does.

Everything so far assumes the whole data set is available. In offline spectroscopy you record the full spectrum first and process it afterwards, so a window can reach both backwards and forwards.

In real-time monitoring that is impossible: the future has not been measured yet. Typical examples are:

• PAT (Process Analytical Technology): real-time NIR or Raman monitoring of pharmaceutical manufacturing.
• Online quality control of production streams.
• Process monitoring of reactors, fermentations or polymerisations.
• Kinetic measurements, where a property is followed over time.

Here we smooth a time series (for example the concentration predicted from each new spectrum), and we need a causal, or one-sided, moving average that uses only the current point and the past:

The window now extends only backwards in time. In the panel below the window appears only to the left of the current point. Increase the window size, slide the window through the changes in the process, and watch the smoothed curve lag further behind the true signal: this delay is the fundamental cost of real-time smoothing. Here even window sizes are fine, because there is no centre to preserve.

A causal moving average delays the signal. For a steady trend, the smoothed value lags the true one by

because the average of the last $w$ time points sits at the middle of that stretch of time.

Example: with $w=9$ and one measurement every 30 seconds, the delay is 4 measurements, so 2 minutes. When the process changes, the smoothed signal shows it about 2 minutes later.

Practical implications:

• Larger windows reduce more noise but delay more.
• In feedback control, delays can make a control loop unstable.
• In alarm systems, delayed detection can be critical.

Choosing the window for real time:

• Fast processes (reactions, mixing): small windows (for example $w=3$ to 5).
• Slow processes (fermentation, curing): larger windows (for example 7 to 11) are acceptable.
• Critical alarms: favour fast detection over smoothness.
• Trending and reporting: delay matters less, so stronger smoothing is fine.

def causal_moving_average(y, window_size):
    """
    One-sided moving average for real-time applications.
    Only uses the current and past values of a time series.
    """
    smoothed = []
    for i in range(len(y)):
        start = max(0, i - window_size + 1)   # window extends backwards only
        # (the first w-1 outputs average fewer points, because there is no earlier data)
        window_values = y[start:i + 1]
        smoothed.append(sum(window_values) / len(window_values))
    return smoothed

# PAT example: one predicted concentration per NIR spectrum, every 30 s
predicted_concentration = [2.31, 2.35, 2.29, 2.40, 2.38, 2.52, 2.47]
trend = causal_moving_average(predicted_concentration, window_size=5)

In a running process you would not recompute the whole list each time: keep the last $w$ values in a buffer and update the average as each new value arrives.

You can, and it is not as silly as it sounds: applying a boxcar twice is the same as applying a single triangular weighting, and after several passes the combined weights look more and more like a bell curve. That is one route to Gaussian smoothing. But each pass also blurs more, so if what you want is simply more smoothing, choose a larger window or a better method rather than stacking passes.

Yes. If you build a calibration model or compare spectra, every spectrum (calibration, validation and test) must be processed with the same window. Inconsistent smoothing creates differences that have nothing to do with chemistry [2],[3],[5].

Check the residuals:

1. Smooth your data.
2. Subtract: $r_i=y_i-{\hat{y}}_i$.
3. Plot the residuals.

Good smoothing: the residuals look like random noise, with no pattern and centred on zero.

Over-smoothing: the residuals show structure, such as the shape of half a band. You are removing signal; use a smaller window.

Averaging a sharp band with its neighbouring baseline pulls the top down. The band becomes shorter and broader, although its area stays the same (exactly, away from the edges), because the intensity is redistributed rather than removed. This is why the moving average is not ideal for quantitative work based on the heights of sharp bands. Savitzky-Golay smoothing preserves band heights much better [1].

Here is how to implement a moving average in Python, MATLAB and R. The example spectrum has one point every 2 nm, so its main band (at 480 nm) is about 16 points wide at half height: a window of 15 is close to the FWHM of that band, so expect it to lose a clear part of its height (about 13 % here; a 7-point window loses about 3 %).

import numpy as np
import matplotlib.pyplot as plt

def moving_average(y, window_size):
    """
    Apply moving average smoothing to a 1D signal.

    Parameters:
    -----------
    y : array-like
        Input signal to be smoothed (e.g., spectrum intensities)
    window_size : int
        Number of points in the averaging window (use odd numbers!)

    Returns:
    --------
    y_smooth : array
        Smoothed signal (same length as input)
    """
    # Create the averaging kernel (array of 1/w values)
    window = np.ones(window_size) / window_size

    # Apply convolution. Note: mode='same' pads with zeros, so the first and
    # last (w-1)/2 values are pulled towards zero (see 'Edge effects')
    y_smooth = np.convolve(y, window, mode='same')

    return y_smooth

# Example: noisy fluorescence emission spectrum
np.random.seed(42)

# Wavelength range typical for fluorescence (nm)
wavelength = np.linspace(400, 600, 100)

# Simulate a realistic fluorescence spectrum with vibronic structure
true_signal = (0.85 * np.exp(-((wavelength - 480)**2) / 400) +
               0.35 * np.exp(-((wavelength - 510)**2) / 300) +
               0.15 * np.exp(-((wavelength - 545)**2) / 250))

# Add realistic shot noise (proportional to sqrt of signal)
noise = np.random.normal(0, 0.08 * np.sqrt(true_signal))
noisy_spectrum = true_signal + noise

# Try different window sizes
smooth_w5 = moving_average(noisy_spectrum, window_size=5)
smooth_w9 = moving_average(noisy_spectrum, window_size=9)
smooth_w15 = moving_average(noisy_spectrum, window_size=15)

# Plot and compare
plt.figure(figsize=(12, 6))
plt.plot(wavelength, noisy_spectrum, 'o-', alpha=0.4,
         label='Noisy data', color='gray', markersize=3)
plt.plot(wavelength, smooth_w5, linewidth=2.5,
         label='Window = 5')
plt.plot(wavelength, smooth_w9, linewidth=2.5,
         label='Window = 9')
plt.plot(wavelength, smooth_w15, linewidth=2.5,
         label='Window = 15')
plt.xlabel('Wavelength (nm)', fontsize=12)
plt.ylabel('Fluorescence Intensity (a.u.)', fontsize=12)
plt.title('Moving Average Smoothing: Window Size Comparison', fontsize=14)
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Check residuals
residuals = noisy_spectrum - smooth_w9
plt.figure(figsize=(12, 4))
plt.plot(wavelength, residuals, 'o-', markersize=3)
plt.axhline(y=0, color='red', linestyle='--', alpha=0.5)
plt.xlabel('Wavelength (nm)')
plt.ylabel('Residual (Original - Smoothed)')
plt.title('Residuals: Should look like random noise with no patterns')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

% Example: fluorescence emission spectrum
% (moving_average is defined at the end: MATLAB requires local functions
% to come after the script code)
rng(42);

% Typical fluorescence wavelength range (nm)
wavelength = linspace(400, 600, 100);

% Simulate realistic fluorescence with vibronic bands
true_signal = 0.85 * exp(-((wavelength - 480).^2) / 400) + ...
              0.35 * exp(-((wavelength - 510).^2) / 300) + ...
              0.15 * exp(-((wavelength - 545).^2) / 250);

% Add realistic shot noise
noise = 0.08 * sqrt(true_signal) .* randn(size(wavelength));
noisy_spectrum = true_signal + noise;

% Try different window sizes
smooth_w5 = moving_average(noisy_spectrum, 5);
smooth_w9 = moving_average(noisy_spectrum, 9);
smooth_w15 = moving_average(noisy_spectrum, 15);

% Visualize the comparison
figure('Position', [100 100 1000 500]);
plot(wavelength, noisy_spectrum, 'o-', 'Color', [0.7 0.7 0.7], ...
     'MarkerSize', 3, 'DisplayName', 'Noisy data');
hold on;
plot(wavelength, smooth_w5, 'LineWidth', 2.5, 'DisplayName', 'Window = 5');
plot(wavelength, smooth_w9, 'LineWidth', 2.5, 'DisplayName', 'Window = 9');
plot(wavelength, smooth_w15, 'LineWidth', 2.5, 'DisplayName', 'Window = 15');
xlabel('Wavelength (nm)', 'FontSize', 12);
ylabel('Fluorescence Intensity (a.u.)', 'FontSize', 12);
title('Moving Average: Window Size Comparison', 'FontSize', 14);
legend('Location', 'northeast', 'FontSize', 10);
grid on;

% Check residuals
figure('Position', [100 100 1000 300]);
residuals = noisy_spectrum - smooth_w9;
plot(wavelength, residuals, 'o-', 'MarkerSize', 3);
yline(0, '--r', 'LineWidth', 1.5);
xlabel('Wavelength (nm)');
ylabel('Residual');
title('Residuals: Should be random noise, no structure');
grid on;

function y_smooth = moving_average(y, window_size)
    % Apply moving average smoothing to a 1D signal
    %
    % Parameters:
    %   y - Input signal to smooth (column or row vector)
    %   window_size - Number of points in averaging window (use odd!)
    %
    % Returns:
    %   y_smooth - Smoothed signal (same size as input)

    % Create the averaging kernel
    window = ones(1, window_size) / window_size;

    % Apply convolution. 'same' pads with zeros at the edges;
    % movmean(y, window_size) shrinks the window there instead
    y_smooth = conv(y, window, 'same');
end

moving_average <- function(y, window_size) {
  # Apply moving average smoothing to a 1D signal
  #
  # Parameters:
  #   y - Input signal to smooth (numeric vector)
  #   window_size - Number of points in window (use odd numbers!)
  #
  # Returns:
  #   y_smooth - Smoothed signal

  # Create averaging kernel
  window <- rep(1 / window_size, window_size)

  # Apply convolution with centered window
  y_smooth <- stats::filter(y, window, sides = 2)

  # Handle edge NAs by using original edge values
  # (R's filter leaves NAs at edges by default)
  y_smooth[is.na(y_smooth)] <- y[is.na(y_smooth)]

  return(as.vector(y_smooth))
}

# Example: fluorescence emission spectrum
set.seed(42)

# Wavelength range for fluorescence (nm)
wavelength <- seq(400, 600, length.out = 100)

# Simulate realistic fluorescence with vibronic structure
true_signal <- 0.85 * exp(-((wavelength - 480)^2) / 400) +
               0.35 * exp(-((wavelength - 510)^2) / 300) +
               0.15 * exp(-((wavelength - 545)^2) / 250)

# Add realistic shot noise
noise <- 0.08 * sqrt(true_signal) * rnorm(length(wavelength))
noisy_spectrum <- true_signal + noise

# Compare different window sizes
smooth_w5 <- moving_average(noisy_spectrum, 5)
smooth_w9 <- moving_average(noisy_spectrum, 9)
smooth_w15 <- moving_average(noisy_spectrum, 15)

# Visualize
par(mfrow = c(1, 1), mar = c(4, 4, 3, 1))
plot(wavelength, noisy_spectrum, type = 'o', col = 'gray60',
     pch = 16, cex = 0.6,
     xlab = 'Wavelength (nm)',
     ylab = 'Fluorescence Intensity (a.u.)',
     main = 'Moving Average: Window Size Comparison')
lines(wavelength, smooth_w5, col = '#2563eb', lwd = 2.5)
lines(wavelength, smooth_w9, col = '#059669', lwd = 2.5)
lines(wavelength, smooth_w15, col = '#dc2626', lwd = 2.5)
legend('topright',
       legend = c('Noisy', 'w=5', 'w=9', 'w=15'),
       col = c('gray60', '#2563eb', '#059669', '#dc2626'),
       lwd = c(1, 2.5, 2.5, 2.5),
       pch = c(16, NA, NA, NA),
       pt.cex = 0.6,
       cex = 0.9)
grid()

# Check residuals
par(mar = c(4, 4, 3, 1))
residuals <- noisy_spectrum - smooth_w9
plot(wavelength, residuals, type = 'o', pch = 16, cex = 0.6,
     xlab = 'Wavelength (nm)', ylab = 'Residual',
     main = 'Residuals: No patterns = good smoothing')
abline(h = 0, col = 'red', lty = 2, lwd = 2)
grid()

Whatever you use it for: plot the data before and after smoothing, check the residuals for structure, use the same odd window for every sample in a study, and check how your routine treats the edges.

• The moving average replaces each point by the mean of a window of $w$ points.
• For independent random noise it divides the noise standard deviation by $\sqrt{w}$.
• The same averaging lowers and broadens any band that is not much wider than the window, so the narrowest important band sets the largest acceptable window.
• In real time, only a one-sided window is possible, and it delays the signal by about $(w-1)/2$ samples.

The moving average teaches the core ideas of every smoother: noise reduction versus feature preservation, window size, and edge effects. For most serious spectroscopic work, consider:

• Gaussian smoothing: a weighted average in which close neighbours matter more, giving smoother results with the same simplicity.
• Savitzky-Golay smoothing: fits local polynomials instead of averaging, preserving band heights much better, and computes derivatives.
• Whittaker smoothing: an optimisation approach that balances fidelity to the data against smoothness with a single parameter.