Filtering API

skelarm.filtering

From-scratch smoothing filters for noisy reference trajectories.

A hand-taught or coarsely sampled reference is often jagged; differentiating it for a velocity/acceleration reference then amplifies the noise. These low-pass filters smooth a uniformly sampled series before resampling/differentiation. They are implemented without SciPy:

• lowpass_first_order — a one-pole RC/exponential filter.
• butterworth_lowpass — an order-n Butterworth filter built from its analog prototype poles via the bilinear transform.

Both are applied zero-phase (forward then backward) so the smoothed signal is not time-shifted, and both are seeded with steady-state initial conditions so a constant (DC) signal passes through unchanged. A series may be 1-D (N,) or multi-column (N, k) (filtered column-wise).

See docs/reference/09_trajectory_filtering.md for the theory.

FILTERS = ('none', 'lowpass', 'butterworth', 'moving_average', 'savgol') module-attribute

The supported filter kinds for :func:smooth.

butterworth_lowpass(values, dt, *, cutoff_hz, order=2)

Zero-phase Butterworth low-pass filter of the given order.

Designs the digital filter from the analog Butterworth poles (bilinear transform with frequency pre-warping) and applies it forward-and-backward (filtfilt-style) with reflected edge padding and steady-state initial conditions.

Parameters:

Name	Type	Description	Default
values	ArrayLike	Series of shape (N,) or (N, k).	required
dt	float	Sample period (s).	required
cutoff_hz	float	Cutoff frequency (Hz); must be below the Nyquist 1 / (2 dt).	required
order	int	Filter order (>= 1); keep small (<= 4).	2

Returns:

Type	Description
NDArray[float64]	The smoothed series, same shape as values.

Raises:

Type	Description
ValueError	If order < 1 or cutoff_hz is not below the Nyquist frequency.

Source code in src/skelarm/filtering.py

def butterworth_lowpass(
    values: ArrayLike,
    dt: float,
    *,
    cutoff_hz: float,
    order: int = 2,
) -> NDArray[np.float64]:
    """Zero-phase Butterworth low-pass filter of the given order.

    Designs the digital filter from the analog Butterworth poles (bilinear transform
    with frequency pre-warping) and applies it forward-and-backward (``filtfilt``-style)
    with reflected edge padding and steady-state initial conditions.

    Parameters
    ----------
    values : ArrayLike
        Series of shape ``(N,)`` or ``(N, k)``.
    dt : float
        Sample period (s).
    cutoff_hz : float
        Cutoff frequency (Hz); must be below the Nyquist ``1 / (2 dt)``.
    order : int, optional
        Filter order (``>= 1``); keep small (``<= 4``).

    Returns
    -------
    NDArray[np.float64]
        The smoothed series, same shape as ``values``.

    Raises
    ------
    ValueError
        If ``order < 1`` or ``cutoff_hz`` is not below the Nyquist frequency.
    """
    if order < 1:
        msg = f"order must be >= 1, got {order}"
        raise ValueError(msg)
    nyquist = 0.5 / dt
    if not 0.0 < cutoff_hz < nyquist:
        msg = f"cutoff_hz must be in (0, {nyquist}) (the Nyquist frequency), got {cutoff_hz}"
        raise ValueError(msg)
    b, a = _butter_lowpass_coeffs(cutoff_hz, dt, order)
    return _apply_per_column(lambda x: _filtfilt(b, a, x), values)

lowpass_first_order(values, dt, *, cutoff_hz)

Zero-phase one-pole (RC) low-pass filter.

Applies the exponential recursion y[i] = y[i-1] + a (x[i] - y[i-1]) forward and backward, with a = dt / (RC + dt) and RC = 1 / (2 pi cutoff_hz). Seeded at the boundary sample so a constant input is preserved.

Parameters:

Name	Type	Description	Default
values	ArrayLike	Series of shape (N,) or (N, k).	required
dt	float	Sample period (s).	required
cutoff_hz	float	Cutoff frequency (Hz).	required

Returns:

Type	Description
NDArray[float64]	The smoothed series, same shape as values.

Source code in src/skelarm/filtering.py

def lowpass_first_order(values: ArrayLike, dt: float, *, cutoff_hz: float) -> NDArray[np.float64]:
    """Zero-phase one-pole (RC) low-pass filter.

    Applies the exponential recursion ``y[i] = y[i-1] + a (x[i] - y[i-1])`` forward and
    backward, with ``a = dt / (RC + dt)`` and ``RC = 1 / (2 pi cutoff_hz)``. Seeded at
    the boundary sample so a constant input is preserved.

    Parameters
    ----------
    values : ArrayLike
        Series of shape ``(N,)`` or ``(N, k)``.
    dt : float
        Sample period (s).
    cutoff_hz : float
        Cutoff frequency (Hz).

    Returns
    -------
    NDArray[np.float64]
        The smoothed series, same shape as ``values``.
    """
    rc = 1.0 / (2.0 * np.pi * cutoff_hz)
    alpha = dt / (rc + dt)

    def pass1(x: NDArray[np.float64]) -> NDArray[np.float64]:
        y = np.empty_like(x)
        y[0] = x[0]
        for i in range(1, x.shape[0]):
            y[i] = y[i - 1] + alpha * (x[i] - y[i - 1])
        return y

    def both(x: NDArray[np.float64]) -> NDArray[np.float64]:
        return pass1(pass1(x)[::-1])[::-1]

    return _apply_per_column(both, values)

moving_average(values, *, window)

Centered (zero-phase) moving-average filter over an odd window of samples.

Each output is the mean of the window samples centered on it; a symmetric window introduces no phase shift and reproduces constants and linear trends in the interior. The signal is reflected at the boundaries so the window stays full near the edges.

Parameters:

Name	Type	Description	Default
values	ArrayLike	Series of shape (N,) or (N, k).	required
window	int	Window length in samples; must be odd and positive.	required

Returns:

Type	Description
NDArray[float64]	The smoothed series, same shape as values.

Raises:

Type	Description
ValueError	If window is not a positive odd integer.

Source code in src/skelarm/filtering.py

def moving_average(values: ArrayLike, *, window: int) -> NDArray[np.float64]:
    """Centered (zero-phase) moving-average filter over an odd ``window`` of samples.

    Each output is the mean of the ``window`` samples centered on it; a symmetric window
    introduces no phase shift and reproduces constants and linear trends in the interior.
    The signal is reflected at the boundaries so the window stays full near the edges.

    Parameters
    ----------
    values : ArrayLike
        Series of shape ``(N,)`` or ``(N, k)``.
    window : int
        Window length in samples; must be odd and positive.

    Returns
    -------
    NDArray[np.float64]
        The smoothed series, same shape as ``values``.

    Raises
    ------
    ValueError
        If ``window`` is not a positive odd integer.
    """
    half = _validate_window(window, polyorder=None, n=np.asarray(values).shape[0])
    kernel = np.full(2 * half + 1, 1.0 / (2 * half + 1))
    return _apply_per_column(lambda x: _correlate_centered(x, kernel, half), values)

savitzky_golay(values, *, window, polyorder=2)

Zero-phase Savitzky-Golay filter: a least-squares polynomial fit over a sliding window.

Each output is the value at the center of a degree-polyorder polynomial fit to the window samples around it. The symmetric kernel is phase-free and reproduces polynomials up to polyorder exactly in the interior, so it smooths noise while preserving peaks better than a plain average. The signal is reflected at the boundaries.

Parameters:

Name	Type	Description	Default
values	ArrayLike	Series of shape (N,) or (N, k).	required
window	int	Window length in samples; must be odd and greater than polyorder.	required
polyorder	int	Polynomial order fit in each window.	2

Returns:

Type	Description
NDArray[float64]	The smoothed series, same shape as values.

Raises:

Type	Description
ValueError	If window is not odd or not greater than polyorder.

Source code in src/skelarm/filtering.py

def savitzky_golay(values: ArrayLike, *, window: int, polyorder: int = 2) -> NDArray[np.float64]:
    """Zero-phase Savitzky-Golay filter: a least-squares polynomial fit over a sliding window.

    Each output is the value at the center of a degree-``polyorder`` polynomial fit to the
    ``window`` samples around it. The symmetric kernel is phase-free and reproduces
    polynomials up to ``polyorder`` exactly in the interior, so it smooths noise while
    preserving peaks better than a plain average. The signal is reflected at the boundaries.

    Parameters
    ----------
    values : ArrayLike
        Series of shape ``(N,)`` or ``(N, k)``.
    window : int
        Window length in samples; must be odd and greater than ``polyorder``.
    polyorder : int, optional
        Polynomial order fit in each window.

    Returns
    -------
    NDArray[np.float64]
        The smoothed series, same shape as ``values``.

    Raises
    ------
    ValueError
        If ``window`` is not odd or not greater than ``polyorder``.
    """
    half = _validate_window(window, polyorder=polyorder, n=np.asarray(values).shape[0])
    offsets = np.arange(-half, half + 1, dtype=np.float64)
    design = np.vander(offsets, polyorder + 1, increasing=True)  # columns z^0 .. z^polyorder
    kernel = np.linalg.pinv(design)[0]  # value at the window center; a symmetric kernel
    return _apply_per_column(lambda x: _correlate_centered(x, kernel, half), values)

smooth(values, dt, *, kind='none', cutoff_hz=None, order=2, window=None, polyorder=2)

Smooth a uniformly sampled series by the named filter kind.

Parameters:

Name	Type	Description	Default
values	ArrayLike	The series, shape (N,) or (N, k).	required
dt	float	Sample period in seconds (the series must be uniformly sampled).	required
kind	str	One of :data:FILTERS. "none" returns the series unchanged.	'none'
cutoff_hz	float | None	Cutoff frequency in Hz; required for "lowpass" and "butterworth".	None
order	int	Butterworth order (ignored otherwise); keep small (<= 4).	2
window	int | None	Window length in samples (odd); required for "moving_average" and "savgol".	None
polyorder	int	Savitzky-Golay polynomial order ("savgol" only); must be < window.	2

Returns:

Type	Description
NDArray[float64]	The smoothed series, same shape as values.

Raises:

Type	Description
ValueError	If kind is unknown or a parameter required by the chosen kind is missing.

Source code in src/skelarm/filtering.py

def smooth(
    values: ArrayLike,
    dt: float,
    *,
    kind: str = "none",
    cutoff_hz: float | None = None,
    order: int = 2,
    window: int | None = None,
    polyorder: int = 2,
) -> NDArray[np.float64]:
    """Smooth a uniformly sampled series by the named filter ``kind``.

    Parameters
    ----------
    values : ArrayLike
        The series, shape ``(N,)`` or ``(N, k)``.
    dt : float
        Sample period in seconds (the series must be uniformly sampled).
    kind : str, optional
        One of :data:`FILTERS`. ``"none"`` returns the series unchanged.
    cutoff_hz : float | None, optional
        Cutoff frequency in Hz; required for ``"lowpass"`` and ``"butterworth"``.
    order : int, optional
        Butterworth order (ignored otherwise); keep small (``<= 4``).
    window : int | None, optional
        Window length in samples (odd); required for ``"moving_average"`` and ``"savgol"``.
    polyorder : int, optional
        Savitzky-Golay polynomial order (``"savgol"`` only); must be ``< window``.

    Returns
    -------
    NDArray[np.float64]
        The smoothed series, same shape as ``values``.

    Raises
    ------
    ValueError
        If ``kind`` is unknown or a parameter required by the chosen kind is missing.
    """
    if kind not in FILTERS:
        msg = f"unknown filter kind {kind!r}; choose from {FILTERS}"
        raise ValueError(msg)
    arr = np.asarray(values, dtype=np.float64)
    if kind == "none":
        return arr.copy()
    if kind in {"lowpass", "butterworth"}:
        if cutoff_hz is None:
            msg = f"filter kind {kind!r} requires a cutoff_hz"
            raise ValueError(msg)
        if kind == "lowpass":
            return lowpass_first_order(arr, dt, cutoff_hz=cutoff_hz)
        return butterworth_lowpass(arr, dt, cutoff_hz=cutoff_hz, order=order)
    if window is None:
        msg = f"filter kind {kind!r} requires a window (in samples)"
        raise ValueError(msg)
    if kind == "moving_average":
        return moving_average(arr, window=window)
    return savitzky_golay(arr, window=window, polyorder=polyorder)