Robust filtering for outliers and spikes

"A moving average asked one liar for his opinion and averaged it in. The median asked the whole room and believed the person standing in the middle."

A Level-Headed AI Agent

Assume you have read Chapter 4 on estimation, so terms like estimator, bias, and robustness are familiar, and Chapter 2 on how physical sensors actually fail (saturation, dropouts, impulse pickup). Those failure modes are the enemy here.

Why the mean breaks and the median does not

The core idea comes from robust statistics: the breakdown point of an estimator is the fraction of arbitrarily bad samples it can tolerate before its output can be pushed to \(\pm\infty\). The sample mean has a breakdown point of \(0\): a single infinite value makes the mean infinite. The sample median has a breakdown point of \(0.5\): you must corrupt half the window before the middle value can be forced arbitrarily far. That single fact is why a median filter shrugs off spikes that a moving average smears.

Concretely, a length-\(w\) sliding median filter replaces each sample \(x_t\) with the median of the window \(\{x_{t-k}, \dots, x_{t+k}\}\), where \(w = 2k+1\). An isolated spike is, by construction, an extreme rank in its window, so it never lands in the middle and is discarded. What survives untouched is a step edge, because on either side of the step a majority of the window agrees. This is the property linear low-pass filters cannot offer: a median filter removes impulses and preserves sharp transitions, whereas a low-pass smooths both. The price is nonlinearity, so there is no transfer function and superposition does not hold; you reason about medians in the sample domain, not the frequency domain of Section 6.4.

The Hampel filter: detect, then decide

A plain median filter has one blunt setting, the window, and it rewrites every sample, so it slightly rounds clean peaks you wanted to keep. The Hampel filter is smarter: it only replaces a sample when that sample is a statistically surprising outlier, and it leaves everything else exactly as measured. It is a decision-directed filter built on two robust statistics per window: the median \(m_t\) as the location estimate, and the median absolute deviation (MAD) as the scale estimate,

The constant \(1.4826\) rescales the MAD so that \(\hat{\sigma}_t\) estimates the standard deviation of clean Gaussian data (the MAD of a standard normal is \(1/1.4826 \approx 0.6745\)). A sample is flagged and replaced by the median when it lies more than \(n_\sigma\) robust standard deviations from the local median:

Both statistics have breakdown point \(0.5\), so the detector itself cannot be fooled by the very spikes it hunts. A single outlier barely perturbs \(m_t\) or \(\text{MAD}_t\), so it stands out cleanly and gets clipped, while the rest of the trace passes through verbatim. The threshold \(n_\sigma\) (commonly \(3\)) tunes aggressiveness, and the window \(w\) tunes how local the notion of "normal" is. This detect-then-replace logic is the small-window sibling of the change and anomaly detection you will study in Chapter 12; the difference is that here the goal is to repair the stream in place, not to raise an alarm.

import numpy as np

def hampel(x, window=7, n_sigma=3.0):
    """Robust spike removal. Replaces flagged outliers with the local median."""
    x = np.asarray(x, dtype=float)
    k = window // 2
    y = x.copy()
    n_replaced = 0
    for t in range(k, len(x) - k):
        seg = x[t - k : t + k + 1]
        m = np.median(seg)
        mad = np.median(np.abs(seg - m))
        sigma = 1.4826 * mad
        if sigma > 0 and abs(x[t] - m) > n_sigma * sigma:
            y[t] = m
            n_replaced += 1
    return y, n_replaced

rng = np.random.default_rng(0)
clean = np.sin(np.linspace(0, 8 * np.pi, 400))
noisy = clean + 0.05 * rng.standard_normal(400)
noisy[[40, 41, 190, 300]] += [6.0, -5.0, 7.0, -8.0]   # inject spikes
repaired, hits = hampel(noisy, window=7, n_sigma=3.0)
print(f"replaced {hits} samples; "
      f"max error clean vs repaired = {np.abs(clean - repaired).max():.3f}")

The code above shows the whole mechanism in about fifteen lines, which is worth doing once to understand it. In production you will not hand-roll the sliding median.

A family of rank and morphological filters

The median is one point on a spectrum of rank-order filters that sort the window and pick an order statistic. The weighted median repeats trusted samples before sorting, letting you bias toward the current sample for lower delay. Min and max filters (the two extreme ranks) are the atoms of morphological filtering: an opening (min then max) erases positive spikes and thin bright bursts, a closing (max then min) erases negative dropouts, and the two compose to strip both polarities while preserving broad structure. Morphological filters are the standard tool for baseline-wander and spike removal in ECG and other biosignals (Chapter 29), because they respect the sharp QRS peaks a linear filter would blunt. When you need the estimate to be robust and smooth, the combination filter (median first to kill spikes, then a short linear smoother from Section 6.1) is the workhorse: reject the outliers with a nonlinear stage, then attenuate the residual broadband noise with a linear one.

Tuning, latency, and what to watch for

Two knobs dominate. The window \(w\) sets robustness and delay: a wider window tolerates longer bursts (up to \(k\) consecutive bad samples for a length-\(2k+1\) median) but a causal implementation costs \(k\) samples of latency, since the current output needs \(k\) future samples. The threshold \(n_\sigma\) sets the false-positive rate: too small and you clip real transients, too large and you miss real spikes. Three failure modes deserve caution. First, if genuine outliers exceed half the window (a long burst), even the median breaks; size \(w\) to the worst expected burst. Second, a degenerate window of identical values gives \(\text{MAD}=0\), so guard the division (the code above tests sigma > 0). Third, and this is the leakage trap the book returns to often: any spike detector you tune must be fit on training data only. Choosing \(n_\sigma\) by peeking at the test set, or letting a whole-dataset MAD leak future statistics into past decisions, inflates your reported denoising quality; use the leakage-safe splits from Chapter 5. Finally, remember that these filters discard information: a "spike" might be the event you care about (an impact, a seizure spike, a gunshot). Robust filtering assumes outliers are corruption; when they are signal, detect and route them instead of deleting them.

What's Next

In Section 6.6, we turn the robustness we just gained into a real-time budget. A median or Hampel filter that needs \(k\) future samples buys its outlier rejection with latency, so the next section is about designing filters that run causally, sample-by-sample, inside a low-latency streaming loop, and about the delay that any such choice imposes before Section 6.7 makes the delay-versus-noise tradeoff explicit.