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 low-pass filter smooths everything it touches, spikes and edges alike, while a median filter smooths only what disagrees with its neighbors, so a genuine edge survives untouched. The price is nonlinearity: there is no transfer function, and 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, for instance, repeats trusted samples before sorting to 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, a closing (max then min) erases negative dropouts, and the two compose to strip both polarities while preserving broad structure. The mechanism is a size argument, not a statistical one: the min pass (erosion) suppresses any feature narrower than the structuring window, spike or genuine peak alike, and the max pass (dilation) restores the width of whatever survived. Size the window narrower than your shortest genuine feature and wider than your longest expected spike, and the opening keeps the feature while deleting the spike; get the sizing wrong and, unlike the Hampel filter, morphological filtering has no statistical safety net. It is the standard tool for baseline-wander and spike removal in ECG and other biosignals (Chapter 29), respecting sharp QRS peaks a linear filter would blunt. When you need both robustness and smoothness, chain a median stage (kill the spikes) into a short linear smoother from Section 6.1 (attenuate the residual broadband noise); this combination filter is the workhorse of practical robust preprocessing.

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, even the median breaks, so size \(w\) to the worst expected burst. Second, a degenerate window of identical values gives \(\text{MAD}=0\); guard the division (the code above tests sigma > 0). Third, and this is a leakage trap the book returns to often: fit any spike detector on training data only, since choosing \(n_\sigma\) by peeking at the test set, or letting a whole-dataset MAD leak future statistics into past decisions, inflates reported denoising quality; use the leakage-safe splits from Chapter 5. Finally, these filters discard information: a "spike" might be the event you care about (an impact, a seizure spike, a gunshot), so when outliers are signal rather than corruption, 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.