Hampel Filter for Real-Time Spike Removal in Sensor Streams

A melt-pressure transducer on an extrusion line reports a valid, slowly-varying signal for hours, then a single scan jumps 40 bar above the surrounding trace for one sample before falling straight back to baseline — electrical noise coupled onto the 4–20 mA loop, not a real process event. This page is the streaming recipe under outlier detection methods for exactly that failure class: the Hampel filter. Unlike a batch Z-score filter, which flags a sample and hands the decision to a downstream imputation stage, the Hampel identifier replaces the spike in place, in a single causal pass, using a rolling median and Median Absolute Deviation (MAD) that the spike itself never gets to contaminate. That in-place correction matters whenever an OEE pipeline consumes the value directly — a cycle-time trigger, a control-loop setpoint, a live HMI trend — rather than a flag column waiting for a separate cleaning stage to resolve later.

Hampel filter: a rejected spike replaced by the local median in a causal streaming pass A raw melt-pressure signal is plotted as a solid line with small process noise and one large upward spike. A shaded band around the local rolling median, spanning plus or minus n times sigma times MAD, brackets the normal noise but not the spike. The spike is marked rejected. A dashed output line tracks the raw signal everywhere except at the spike, where it stays flat at the local median value, marked replaced, showing the filter's in-place correction rather than a flag. Melt pressure (bar) time → median + nσ·MAD median − nσ·MAD rejected · |x−median| > nσ·MAD replaced with median raw stream Hampel output

The Hampel identifier evaluates a candidate sample against a robust local baseline built entirely from its neighbors, so one bad reading cannot inflate the very statistic used to judge it — the same defect that makes plain mean/standard-deviation thresholding fragile against the impulsive spikes covered in the modified Z-score variant of the vibration page. Where that page’s variants flag a sample and pass a NaN or a confidence score downstream, the Hampel filter’s job is narrower and more immediate: decide, on every incoming sample, whether to pass it through unchanged or substitute the local median, with no separate imputation pass required.

The Hampel identifier: rolling median, MAD, and n-sigma replacement

For a window of 2k + 1 samples centered on point xix_i, the identifier computes the median x~i\tilde{x}_i and the Median Absolute Deviation:

MADi=median(∣xi−k,\dots,xi+k−x~i∣)\text{MAD}_i = \text{median}\left(\left|x_{i-k},\dots,x_{i+k} - \tilde{x}_i\right|\right)

Scaled by the constant 1.4826, MAD is a consistent estimator of σ\sigma under a Gaussian assumption, so the rejection rule reads as an ordinary n-sigma test built on breakdown-resistant statistics:

∣xi−x~i∣>n⋅1.4826⋅MADi\impliesxi←x~i\left|x_i - \tilde{x}_i\right| > n \cdot 1.4826 \cdot \text{MAD}_i \implies x_i \leftarrow \tilde{x}_i

A batch implementation over a pandas.Series makes the mechanics explicit before the streaming version below carries the same logic forward sample by sample:

import numpy as np
import pandas as pd


def hampel_filter_batch(
    series: pd.Series,
    window: int = 7,
    n_sigma: float = 3.0,
) -> tuple[pd.Series, pd.Series]:
    """
    Offline Hampel filter over a full series. `window` is the half-width k;
    the effective window is 2k+1 samples, centered on each point. Returns
    (cleaned_series, replaced_mask).
    """
    if window < 1:
        raise ValueError("window (half-width) must be >= 1")

    median = series.rolling(window=2 * window + 1, center=True, min_periods=2 * window + 1).median()
    abs_dev = (series - median).abs()
    mad = abs_dev.rolling(window=2 * window + 1, center=True, min_periods=2 * window + 1).median()
    robust_sigma = (mad * 1.4826).replace(0.0, np.nan)  # guard: flat window -> no rejection

    deviation = (series - median).abs()
    replaced = (deviation > n_sigma * robust_sigma).fillna(False)

    cleaned = series.copy()
    cleaned[replaced] = median[replaced]
    return cleaned, replaced

center=True is the detail that separates this from the causal filter a live pipeline actually needs: a centered window looks k samples into the future, which is fine for reprocessing historian data but impossible for a value that must be published the instant it is produced.

Streaming implementation with carried window state

A real-time gateway cannot wait for future samples that have not arrived yet. The standard fix is a causal, lagged Hampel filter: keep a fixed-size buffer of the most recent 2k + 1 readings and evaluate the sample sitting k positions behind the newest arrival, which already has k samples of “future” context sitting in the buffer. This introduces a deterministic latency of k sample periods — the cost of getting a centered statistic without waiting on the wall clock — and it is the same trade the causal variant of clock drift correction makes when it delays a correction to get a stable estimate rather than reacting to every raw tick.

from __future__ import annotations

from collections import deque
from dataclasses import dataclass, field
from typing import Deque, Optional

import numpy as np


@dataclass
class HampelState:
    """Causal window buffer carried between samples, or serialized across
    async batch boundaries — the same pattern used to stitch rolling
    statistics across chunks in the vibration Z-score pipeline."""
    half_window: int
    buffer: Deque[float] = field(default_factory=deque)

    def push(self, value: float) -> None:
        self.buffer.append(value)
        max_len = 2 * self.half_window + 1
        while len(self.buffer) > max_len:
            self.buffer.popleft()


def hampel_step(state: HampelState, n_sigma: float = 3.0) -> tuple[Optional[float], bool]:
    """
    Evaluate the sample lagged `half_window` positions behind the newest
    arrival. Returns (output_value, was_replaced); output is None while the
    buffer is still filling on startup.
    """
    window_len = 2 * state.half_window + 1
    if len(state.buffer) < window_len:
        return None, False  # startup latency: not enough context yet

    values = np.fromiter(state.buffer, dtype=np.float64, count=len(state.buffer))
    center_value = values[state.half_window]

    median = float(np.median(values))
    mad = float(np.median(np.abs(values - median)))
    robust_sigma = 1.4826 * mad

    if robust_sigma < 1e-9:
        return center_value, False  # genuinely flat window; nothing to reject

    if abs(center_value - median) > n_sigma * robust_sigma:
        return median, True
    return center_value, False


class HampelStream:
    """Thin wrapper for per-sensor use inside an MQTT consumer loop."""

    def __init__(self, half_window: int = 5, n_sigma: float = 3.0) -> None:
        self._state = HampelState(half_window=half_window)
        self._n_sigma = n_sigma

    def feed(self, value: float) -> tuple[Optional[float], bool]:
        self._state.push(value)
        return hampel_step(self._state, self._n_sigma)

Two production details carry weight here. First, np.fromiter over the deque avoids repeatedly reallocating a list on every sample at high scan rates. Second, HampelState is deliberately small and serializable (a bounded deque of floats) so it can be pickled into Redis or a per-sensor key-value store between micro-batches, the same carry-forward discipline that keeps rolling statistics continuous across chunk boundaries in async ingestion — see the batch-state pattern in async batch processing. Without it, every new batch starts its buffer empty and re-pays the k-sample startup latency at every chunk edge, silently widening the delay before any correction can fire.

Tuning window and threshold per signal type

half_window (k) and n_sigma are not universal constants — they trade detection latency against sensitivity, and the right values differ by an order of magnitude between a fast, jittery flow signal and a slow thermal one. A cooling-water flow meter sampled at 10 Hz with normal turbulence noise wants a short window (k = 4–6) so single-scan electrical spikes are caught within half a second, and a looser n_sigma = 3.5–4.0 so ordinary turbulence is not mistaken for a fault. A thermocouple on a barrel zone sampled at 1 Hz changes slowly by physical necessity; a short window there rejects legitimate but rapid heater-cycling transients, so k = 8–12 with a tighter n_sigma = 2.5–3.0 catches thermocouple-wire noise spikes without smearing real set-point changes.

The failure mode in both directions is asymmetric. Too small a k starves the median/MAD estimate of samples, so it tracks the signal too closely and the filter never rejects anything — effectively a no-op. Too large a k drags in samples from a different operating regime (a load change, a changeover), inflating MAD until genuine spikes fall inside the band and pass through uncorrected. As a starting point, size k from the physical settling time of the process relative to the sample interval, then validate against a labeled window of historian data the way the parent page’s confusion-matrix replay recommends, rather than guessing a single number and shipping it.

Gotchas & anti-patterns

  • Treating n_sigma as portable across signals. A threshold tuned for a noisy flow meter under-rejects on a quiet pressure loop and over-rejects on a jittery current sensor. Tune per sensor profile, not per pipeline.
  • Ignoring the causal lag in latency-sensitive control loops. The output at time t reflects the sample from t − k scan intervals ago. Feeding that lagged value into a tight PID loop, rather than into logging or a dashboard, introduces phase error the loop was never tuned for.
  • Resetting the buffer on every deploy or batch boundary. A fresh, empty HampelState re-pays the startup latency and passes the first k samples of every restart unfiltered. Persist and reload state across process restarts, not just across batches.
  • Comparing floats for the flat-window guard with == 0. IEEE 754 doubles rarely land on exactly zero after arithmetic; use a small epsilon (robust_sigma < 1e-9 above) rather than exact equality, matching the tolerance discipline in precision and rounding limits.
  • Applying Hampel to a signal that is legitimately flat. A truly frozen or stuck-at sensor produces MAD ≈ 0 and every sample looks identical to the median — the Hampel filter has nothing to reject and silently passes a dead transducer through. That failure needs a dedicated stuck-at and frozen sensor test, not a spike filter.

Quick reference

Parameter Effect of increasing it Effect of decreasing it Typical starting point
half_window (k) More stable median/MAD; more latency (k samples) Faster response; noisier baseline, may under-reject 4–6 for fast/noisy signals, 8–12 for slow thermal
n_sigma Fewer false rejections; may miss small real spikes Catches smaller spikes; more false rejections 3.0 general default
Sample-rate mismatch N/A Fixed k in samples means variable time-window if rate drifts Re-derive k after any rate change
Buffer persistence Removes restart/batch-edge startup latency Re-pays k-sample latency at every reset Always persist across batches

Related in this section