Smoothing noisy heading signals for turn detection

A raw bearing series from urban GPS is far too noisy to threshold directly, and the obvious fix — a rolling mean — is wrong twice over: it operates on a circular quantity as though it were linear, and it rounds off exactly the corners the detector exists to find. Smooth in vector space, prefer a median filter for impulsive noise, and size the window from the shortest manoeuvre you intend to detect rather than from how clean the output looks.

Why this happens

Bearing noise is the position noise divided by the step length. At a 10 m step and a 3 m positional error the bearing uncertainty is about 17°, which means a stationary-looking urban trace can swing thirty degrees between consecutive fixes with no vehicle motion at all. Accumulating that unfiltered clears any turn threshold within a few seconds.

The distribution matters as much as the magnitude. Multipath produces occasional large excursions rather than uniform jitter, so the noise is impulsive — and a mean filter smears an impulse across the whole window while a median filter removes it entirely. That single property makes the median the better default for headings, which is the opposite of the usual advice for smoothing positions. The bearings themselves come from computing bearing and heading change with NumPy.

The median keeps the turn's edge because a corner is not an outlier — it is a sustained change, and a median filter only removes values that are in a minority within the window.

Core pipeline

Convert bearings to unit vectors so every subsequent operation is linear rather than circular.
Apply a median filter first to remove impulsive multipath excursions.
Apply a short mean or Savitzky-Golay pass if residual jitter remains, sized below the shortest manoeuvre.
Convert back and re-derive the wrapped differences, never smoothing the differences themselves.

Production-ready Python implementation

PYTHON

import numpy as np
import pandas as pd
from scipy.ndimage import median_filter
from scipy.signal import savgol_filter


def smooth_heading(
    bearing_deg: pd.Series,
    median_window: int = 5,
    savgol_window: int = 0,
    polyorder: int = 2,
) -> pd.Series:
    """
    Smooth a bearing series in vector space, preserving corners.

    Parameters
    ----------
    bearing_deg : pd.Series
        Bearings in degrees. NaNs are preserved rather than interpolated —
        a missing heading is information, and filling it invents a direction.
    median_window : int
        Odd window for the impulse-removal pass. 5 at 1 Hz removes single-
        and double-fix multipath excursions without touching a real turn.
    savgol_window : int
        Optional second pass; 0 disables it. Keep it BELOW the number of
        samples in your shortest manoeuvre or the corner is rounded off.

    Returns
    -------
    pd.Series
        Smoothed bearings in [0, 360), aligned to the input index.

    Raises
    ------
    ValueError
        On even windows, a savgol window that exceeds the series, or a
        polyorder that is not smaller than the window.
    """
    if median_window % 2 == 0:
        raise ValueError("median_window must be odd.")
    if savgol_window and savgol_window % 2 == 0:
        raise ValueError("savgol_window must be odd.")
    if savgol_window and savgol_window <= polyorder:
        raise ValueError("savgol_window must exceed polyorder.")

    a = np.radians(bearing_deg.to_numpy(dtype=float))
    valid = np.isfinite(a)
    if valid.sum() < max(median_window, savgol_window or 0):
        raise ValueError("Too few valid bearings to filter with these windows.")

    # ── vector space ──────────────────────────────────────────────────
    # Smoothing sin and cos separately is the only correct way to filter a
    # circular quantity. Filtering degrees directly averages 359 and 1 to 180.
    s, c = np.sin(a), np.cos(a)

    # median_filter needs finite input; work on the valid subset and write back.
    s_f, c_f = s.copy(), c.copy()
    s_f[valid] = median_filter(s[valid], size=median_window, mode="nearest")
    c_f[valid] = median_filter(c[valid], size=median_window, mode="nearest")

    if savgol_window:
        if valid.sum() >= savgol_window:
            s_f[valid] = savgol_filter(s_f[valid], savgol_window, polyorder, mode="interp")
            c_f[valid] = savgol_filter(c_f[valid], savgol_window, polyorder, mode="interp")

    out = np.degrees(np.arctan2(s_f, c_f)) % 360.0
    out[~valid] = np.nan          # do not manufacture headings where there were none
    return pd.Series(out, index=bearing_deg.index, name="bearing_smooth_deg")


def wrapped_diff(bearing_deg: pd.Series) -> pd.Series:
    """Consecutive heading change in (-180, 180]. Re-derive AFTER smoothing."""
    raw = bearing_deg.diff().to_numpy()
    return pd.Series((raw + 180.0) % 360.0 - 180.0,
                     index=bearing_deg.index, name="dheading_deg")

Validation block

PYTHON

def validate_smoothing(raw: pd.Series, smooth: pd.Series,
                       corner_idx: int, expected_turn_deg: float = 90.0) -> None:
    """Check that noise fell and the corner survived — both, not one."""
    dr, ds = wrapped_diff(raw).abs(), wrapped_diff(smooth).abs()

    # 1. Noise reduction: the median absolute change on straight sections
    #    should fall substantially.
    quiet = ds.index.difference(ds.index[max(corner_idx - 6, 0): corner_idx + 6])
    red = 1 - ds.loc[quiet].median() / max(dr.loc[quiet].median(), 1e-9)
    assert red > 0.5, f"only {red:.0%} noise reduction — widen the median window"

    # 2. Corner preservation: the total turn across the manoeuvre must
    #    survive. This is the assertion that catches over-smoothing, and it
    #    is the one usually missing.
    window = slice(max(corner_idx - 6, 0), corner_idx + 6)
    turned = ds.iloc[window].sum()
    assert turned > 0.8 * expected_turn_deg, (
        f"corner lost: {turned:.0f}° of an expected {expected_turn_deg:.0f}° "
        "— the smoothing window is longer than the manoeuvre"
    )

    # 3. No manufactured headings where the input had none.
    assert smooth[raw.isna()].isna().all(), "smoothing filled gaps it should not have"
    print(f"OK — noise −{red:.0%}, corner retained {turned:.0f}°")

The two curves cross where the window approaches the manoeuvre duration. Tuning on the noise curve alone always lands to the right of that point, which is why the corner-retention assertion belongs in the test suite.

Filtering the sine and cosine components separately is what makes every operation on this page correct across the north discontinuity, and it hands you the consistency measure at no extra cost.

Common mistakes and gotchas

Rolling-mean on degrees. The single most common bug. Any window straddling north produces a smoothed heading pointing the opposite way.
Smoothing the differences instead of the bearings. Filtering dheading blurs the turn across neighbouring fixes and makes the cumulative angle depend on the filter rather than the road. Smooth the bearing, then differentiate.
Interpolating over missing headings. A gap in the bearing series means the vehicle was stationary or the fix was rejected. Filling it invents a direction that then feeds the accumulator.
One window for all speeds. A four-second turn is four samples at 1 Hz and one at 0.25 Hz. Express the window in seconds and convert, or the same configuration over-smooths half your fleet.
Savitzky-Golay as the first pass. A polynomial fit is pulled by a multipath spike; a median filter is not. Median first, polynomial second if at all.
Tuning until the plot looks clean. A perfectly clean heading trace is one with no turns in it. Tune against corner retention on a labelled manoeuvre.

FAQ

Median or mean filter?

Median, for headings. GPS heading noise is impulsive rather than Gaussian — a few large excursions rather than uniform jitter — and the median removes an impulse entirely while the mean spreads it over the window. Use a short mean or Savitzky-Golay pass afterwards only if a visible residual ripple remains.

How do I pick the window in seconds?

Take the shortest manoeuvre you must detect, halve it, and use that. A junction turn takes three to five seconds, so a window under two seconds is safe; a lane change takes about two, so detecting those requires a one-second window and correspondingly noisier output. There is no window that detects both cleanly, which is an argument for running two detectors rather than compromising one.

Does this help a Kalman filter?

If you are already running a position-domain Kalman filter, derive the heading from the smoothed positions rather than filtering the heading separately — filtering twice suppresses turns twice. Heading-domain smoothing is for pipelines whose positions are used raw, which is common when the geometry has to stay faithful for map matching.

Directionality & Turn Analysis — the parent reference and calibration tables.
Computing Bearing and Heading Change with NumPy — where the input series comes from.
Detecting U-turns and Directional Shifts in Fleet Data — the detector that consumes the smoothed series.
Kalman Filter vs Savitzky-Golay for Real-Time vs Batch GPS Smoothing — the position-domain equivalents.
Handling GPS Drift in Raw Trajectory Logs — the upstream cleaning that reduces how much of this is needed.

Back to Directionality & Turn Analysis