visualdynamics.core.filters

filters

Filtering, integration and differentiation of records.

These make up the shock workflow's chain: an acceleration record is low-pass filtered, integrated to velocity, and integrated again to displacement. Each step is a derived object with provenance, so a changed filter corner marks everything after it in the chain out of date, the same way a changed averaging does.

The whole record is processed, never the shock windows one at a time, for three reasons. A zero-phase filter's edge transients land at the record's ends, far from the events, where filtering each window would put them on the data that matters. Integrating each window would give each one its own arbitrary velocity baseline. And a derived record that keeps the source's abscissa keeps its marks: the shock windows and averaging frames carry over unchanged, so a report numbers the same events on the displacement as on the acceleration.

Integration and differentiation are done in the time domain. The frequency-domain route (dividing by iω) treats the record as periodic, so the FFT turns the mismatch between the first and last samples into content. The trapezoid rule and the central difference treat the record as the finite length it is. Their cost is a known amplitude roll-off near Nyquist, sin(ωT)/(ωT) for the central difference, which is why the chain filters first: with a corner at a tenth of the sample rate, the round trip from acceleration to displacement and back to acceleration comes back within about 2 % RMS (about 4 % with a corner at a fifth).

Differentiation is np.gradient (second-order central differences), not the exact inverse of the trapezoid rule. That inverse, x[i] = 2*(y[i]-y[i-1])/dt - x[i-1], has a pole at Nyquist in its recursion, and on real noise its error reached 280 % where the central difference stayed under 4 %.

References

The filter forms are textbook; these are where they and the drift a double integration invents are set out.

1. Oppenheim, A. V., & Schafer, R. W. (2010). Discrete-Time Signal Processing, 3rd ed. Prentice Hall. The bilinear transform, the filter forms and zero-phase forward-backward filtering.
2. ISO 18431-1:2005, Mechanical vibration and shock - Signal processing - Part 1: General introduction, and ISO 18431-2 on time-domain integration.
3. Boore, D. M., & Bommer, J. J. (2005). "Processing of strong-motion accelerograms: needs, options and consequences." Soil Dynamics and Earthquake Engineering. The drift a double integration invents, and what filtering does and does not cure.

Classes:

Name	Description
Filtering	The pass band a record is filtered to: its edges and order.

Functions:

Name	Description
threaded_sosfiltfilt	scipy.signal.sosfiltfilt over a record's channels, one channel
design	The filter at sample_rate, as second-order sections from
carrying_marks	Copy the source's averaging and shock windows onto a derived
filtered	The record through its filter, every channel, zero phase.
integrate	One integration, acceleration to velocity or velocity to
differentiate	One differentiation, displacement to velocity or velocity to
response	(frequencies, magnitude in dB) of the filter as the record sees

Classes

Filtering dataclass

Filtering(*, low: float | None = None, high: float | None = None, order: int = 4)

The pass band a record is filtered to: its edges and order.

low and high are the edges of the pass band, in Hz. Set high alone for a low-pass (everything below it passes), low alone for a high-pass, and both for a band-pass. Which filter it is follows from which edges are set (kind), so an invalid pairing, such as a band-pass with one corner, cannot be written. The fields are keyword-only, so a bare Filtering(200.0) is an error. At least one edge is required, edges must be above zero, low must be below high, and order must be at least 1.

The settings are frozen, like Averaging, because they are stored with the history and their fields are its staleness fingerprint; a setting that could change in place would make that fingerprint wrong.

The filter is a Butterworth applied forward and backward (sosfiltfilt). It has zero phase, so a shock's peak stays at its measured time; a causal filter would shift and smear the front by its group delay. The double pass doubles the effective order; order is the design order of the underlying filter. A band-pass of design order N has N poles per edge (scipy's convention), so each skirt rolls off like the low- or high-pass of the same order.

Methods:

Name	Description
describe	The filter in words, such as 'low-pass at 500 Hz'. The status

Attributes:

Name	Type	Description
kind	str	'low-pass', 'high-pass' or 'band-pass', from which edges are

Attributes

kind property

kind: str

'low-pass', 'high-pass' or 'band-pass', from which edges are set; never stored beside them.

Methods:

describe

describe() -> str

The filter in words, such as 'low-pass at 500 Hz'. The status line, the staleness message and the report text all use this wording.

Source code in src/visualdynamics/core/filters.py

def describe(self) -> str:
    """The filter in words, such as 'low-pass at 500 Hz'. The status
    line, the staleness message and the report text all use this
    wording."""
    if self.low is None:
        return f'low-pass at {self.high:g} Hz'
    if self.high is None:
        return f'high-pass at {self.low:g} Hz'
    return f'band-pass from {self.low:g} to {self.high:g} Hz'

Functions:

threaded_sosfiltfilt

threaded_sosfiltfilt(sos: Any, rows: Any, axis: int = -1) -> ndarray

scipy.signal.sosfiltfilt over a record's channels, one channel per thread.

The filter releases the interpreter lock, so a thread pool over the rows is about six times faster on a 32-channel record. The results are the serial ones to the bit, since each row is the same call. A single row, an array under a million values, or filtering along another axis than the last goes to one serial call.

Source code in src/visualdynamics/core/filters.py

def threaded_sosfiltfilt(sos: Any, rows: Any, axis: int = -1) -> np.ndarray:
    """`scipy.signal.sosfiltfilt` over a record's channels, one channel
    per thread.

    The filter releases the interpreter lock, so a thread pool over the
    rows is about six times faster on a 32-channel record. The results
    are the serial ones to the bit, since each row is the same call. A
    single row, an array under a million values, or filtering along
    another axis than the last goes to one serial call.
    """
    # Threaded on 2026-10-01 (Brandon: use the cores).
    from scipy.signal import sosfiltfilt

    rows = np.asarray(rows)
    if rows.ndim < 2 or axis not in (-1, rows.ndim - 1) or rows.shape[0] < 2 \
            or rows.size < 1_000_000:
        return sosfiltfilt(sos, rows, axis=axis)
    return np.stack(list(_filter_pool().map(
        lambda row: sosfiltfilt(sos, row), rows)))

design

design(filtering: Filtering, sample_rate: float) -> Any

The filter at sample_rate, as second-order sections from scipy.signal.butter.

This is the one place a Filtering becomes a scipy design: filtering a record, the response plot and the GUI previews all take their sections from here, so they agree on what the settings mean. An edge at or above the Nyquist frequency raises ValueError here, so every caller refuses it the same way.

Source code in src/visualdynamics/core/filters.py

def design(filtering: Filtering, sample_rate: float) -> Any:
    """The filter at `sample_rate`, as second-order sections from
    `scipy.signal.butter`.

    This is the one place a `Filtering` becomes a scipy design:
    filtering a record, the response plot and the GUI previews all take
    their sections from here, so they agree on what the settings mean.
    An edge at or above the Nyquist frequency raises ValueError here,
    so every caller refuses it the same way.
    """
    from scipy.signal import butter

    rate = float(sample_rate)
    for edge in (filtering.low, filtering.high):
        if edge is not None and not edge < rate / 2.0:
            raise ValueError(
                f'the corner ({edge:g} Hz) must sit below the '
                f'Nyquist frequency ({rate / 2.0:g} Hz)')
    if filtering.low is None:
        return butter(filtering.order, filtering.high, 'lowpass',
                      fs=rate, output='sos')
    if filtering.high is None:
        return butter(filtering.order, filtering.low, 'highpass',
                      fs=rate, output='sos')
    return butter(filtering.order, [filtering.low, filtering.high],
                  'bandpass', fs=rate, output='sos')

carrying_marks

carrying_marks(source: Any, result: Any) -> Any

Copy the source's averaging and shock windows onto a derived record, and return the record.

The derived record has the same abscissa, so the marks still mean what they meant: the averaging frames and shock windows land on the same samples, and a report numbers the same events on the displacement as on the acceleration, or on the modal coordinates as on the channels, since a transform through a shape set combines samples and moves none. The filtering is not carried: the result of a filter is already filtered, and offering to filter it again with the same settings would be a button that does nothing.

Source code in src/visualdynamics/core/filters.py

def carrying_marks(source: Any, result: Any) -> Any:
    """Copy the source's averaging and shock windows onto a derived
    record, and return the record.

    The derived record has the same abscissa, so the marks still mean
    what they meant: the averaging frames and shock windows land on the
    same samples, and a report numbers the same events on the
    displacement as on the acceleration, or on the modal coordinates as
    on the channels, since a transform through a shape set combines
    samples and moves none. The filtering is not carried: the result of
    a filter is already filtered, and offering to filter it again with
    the same settings would be a button that does nothing.
    """
    # Modal coordinates carry the marks too (Brandon, 2026-09-04).
    result.averaging = source.averaging
    # a tuple, matching `core.shocks.find` and the loader: the windows
    # are a fixed set for a record, and copying rather than sharing
    # keeps an edit to one record's marks off another's
    result.shocks = tuple(source.shocks) if source.shocks else source.shocks
    return result

filtered

filtered(data: TimeHistory, filtering: Filtering) -> TimeHistory

The record through its filter, every channel, zero phase.

Every channel is filtered whatever it measures: a filter acts on the waveform, not the quantity, and the drive force belongs at the same bandwidth as the responses it drove. The averaging and shock windows are carried over (carrying_marks). A record with uneven sampling, or an edge at or above Nyquist, is refused.

Source code in src/visualdynamics/core/filters.py

def filtered(data: TimeHistory, filtering: Filtering) -> TimeHistory:
    """The record through its filter, every channel, zero phase.

    Every channel is filtered whatever it measures: a filter acts on
    the waveform, not the quantity, and the drive force belongs at the
    same bandwidth as the responses it drove. The averaging and shock
    windows are carried over (`carrying_marks`). A record with uneven
    sampling, or an edge at or above Nyquist, is refused.
    """

    from .data import TimeHistory

    rate = data.sample_rate       # raises on uneven sampling, correctly
    sos = design(filtering, rate)
    result = TimeHistory(
        data.abscissa, threaded_sosfiltfilt(sos, np.real(data.ordinate)),
        **data.record_fields())
    return carrying_marks(data, result)

integrate

integrate(data: TimeHistory, drift_corner: float | None = DRIFT_CORNER) -> TimeHistory

One integration, acceleration to velocity or velocity to displacement, channel by channel.

The cumulative trapezoid rule, with drift control: each channel's mean is removed before integrating, and a second-order Butterworth high-pass at drift_corner in Hz (DRIFT_CORNER, 2 Hz) is applied forward and backward afterward. The mean is removed first because a bias integrates to a ramp, and the high-pass's edge transients grow with how far the ramp runs. Below the corner nothing in the data tells drift from true motion. Pass None to integrate raw, which says the record's own low end is trusted, drift and all.

Records of quantities the chain cannot take (forces, voltages, anything undefined) are left out of the result. A record with no acceleration or velocity channels is refused, as is a drift corner that is not between zero and Nyquist.

Source code in src/visualdynamics/core/filters.py

def integrate(data: TimeHistory,
              drift_corner: float | None = DRIFT_CORNER) -> TimeHistory:
    """One integration, acceleration to velocity or velocity to
    displacement, channel by channel.

    The cumulative trapezoid rule, with drift control: each channel's
    mean is removed before integrating, and a second-order Butterworth
    high-pass at `drift_corner` in Hz (`DRIFT_CORNER`, 2 Hz) is applied
    forward and backward afterward. The mean is removed first because a
    bias integrates to a ramp, and the high-pass's edge transients grow
    with how far the ramp runs. Below the corner nothing in the data
    tells drift from true motion. Pass ``None`` to integrate raw, which
    says the record's own low end is trusted, drift and all.

    Records of quantities the chain cannot take (forces, voltages,
    anything undefined) are left out of the result. A record with no
    acceleration or velocity channels is refused, as is a drift corner
    that is not between zero and Nyquist.
    """
    from scipy.integrate import cumulative_trapezoid
    from scipy.signal import butter

    rate = data.sample_rate
    kept = _mapped(data, INTEGRATED)
    if not kept:
        raise ValueError(
            'nothing here integrates: the record holds no acceleration '
            'or velocity channels')
    if drift_corner is not None and not 0.0 < drift_corner < rate / 2.0:
        raise ValueError(
            f'the drift corner ({drift_corner:g} Hz) must sit between '
            f'zero and the Nyquist frequency ({rate / 2.0:g} Hz)')
    rows = np.real(data.ordinate[kept])
    if drift_corner is not None:
        rows = rows - rows.mean(axis=-1, keepdims=True)
    rows = cumulative_trapezoid(rows, dx=1.0 / rate, initial=0.0, axis=-1)
    if drift_corner is not None:
        sos = butter(2, drift_corner, 'high', fs=rate, output='sos')
        rows = threaded_sosfiltfilt(sos, rows)
    return _derived(data, rows, kept, INTEGRATED)

differentiate

differentiate(data: TimeHistory) -> TimeHistory

One differentiation, displacement to velocity or velocity to acceleration, channel by channel.

Second-order central differences (np.gradient); the module docstring says why the exact inverse of the trapezoid rule is not used. There is no drift control, since differentiation removes offsets instead of growing them. Records of other quantities are left out, and a record with no displacement or velocity channels is refused.

Source code in src/visualdynamics/core/filters.py

def differentiate(data: TimeHistory) -> TimeHistory:
    """One differentiation, displacement to velocity or velocity to
    acceleration, channel by channel.

    Second-order central differences (`np.gradient`); the module
    docstring says why the exact inverse of the trapezoid rule is not
    used. There is no drift control, since differentiation removes
    offsets instead of growing them. Records of other quantities are
    left out, and a record with no displacement or velocity channels
    is refused.
    """
    rate = data.sample_rate
    kept = _mapped(data, DIFFERENTIATED)
    if not kept:
        raise ValueError(
            'nothing here differentiates: the record holds no '
            'displacement or velocity channels')
    rows = np.gradient(np.real(data.ordinate[kept]), 1.0 / rate, axis=-1)
    return _derived(data, rows, kept, DIFFERENTIATED)

response

response(filtering: Filtering, sample_rate: float, lines: int = 400, *, frequencies: Any = None) -> tuple[ndarray, ndarray]

(frequencies, magnitude in dB) of the filter as the record sees it.

This is the squared magnitude |H|², not the single-pass response. filtered runs the filter forward and backward for zero phase, so the data is multiplied by |H| twice. A Butterworth is -3 dB at its nominal corner by definition, so the zero-phase pair is -6.02 dB there, which is not what the settings suggest; for a fourth-order low-pass (the default order) the -3 dB point is at about 0.90 of the corner. Plotting |H| would draw a filter that is not the one applied.

The lines frequencies are log-spaced from RESPONSE_DECADES (3) decades below Nyquist up to Nyquist, extended down to half the lowest edge when an edge is lower than that. Edges at or above Nyquist are clamped just under it instead of refused, so a curve can be drawn while the settings are being edited. Values are floored at -240 dB.

frequencies, when given, are the frequencies the response is read at instead, in Hz, so a reading can ask for the gain at the lines it names (core.resample reads it at the span and at the lowest frequency that folds into it).

Source code in src/visualdynamics/core/filters.py

def response(filtering: Filtering, sample_rate: float,
             lines: int = 400, *, frequencies: Any = None
             ) -> tuple[np.ndarray, np.ndarray]:
    """(frequencies, magnitude in dB) of the filter as the record sees
    it.

    This is the squared magnitude |H|², not the single-pass response.
    `filtered` runs the filter forward and backward for zero phase, so
    the data is multiplied by |H| twice. A Butterworth is -3 dB at its
    nominal corner by definition, so the zero-phase pair is -6.02 dB
    there, which is not what the settings suggest; for a fourth-order
    low-pass (the default order) the -3 dB point is at about 0.90 of
    the corner. Plotting |H| would draw a filter that is not the one
    applied.

    The `lines` frequencies are log-spaced from `RESPONSE_DECADES` (3)
    decades below Nyquist up to Nyquist, extended down to half the
    lowest edge when an edge is lower than that. Edges at or above
    Nyquist are clamped just under it instead of refused, so a curve
    can be drawn while the settings are being edited. Values are
    floored at -240 dB.

    `frequencies`, when given, are the frequencies the response is
    read at instead, in Hz, so a reading can ask for the gain at the
    lines it names (`core.resample` reads it at the span and at the
    lowest frequency that folds into it).
    """
    # Squared magnitude, not the filter's own (Brandon, 2026-08-25).
    from scipy.signal import sosfreqz

    nyquist = float(sample_rate) / 2.0
    edges = [e for e in (filtering.low, filtering.high)
             if e is not None]
    if frequencies is None:
        start = min(nyquist / 10.0 ** RESPONSE_DECADES, min(edges) * 0.5)
        frequencies = np.logspace(np.log10(start), np.log10(nyquist),
                                  int(lines))
    else:
        frequencies = np.atleast_1d(np.asarray(frequencies, dtype=float))
    # edges clamped rather than refused: mid-edit the boxes pass
    # through states the filter itself would refuse, and the panel
    # wants a curve to draw against, not an exception
    safe = replace(
        filtering,
        low=(None if filtering.low is None
             else min(filtering.low, nyquist * 0.998)),
        high=(None if filtering.high is None
              else min(filtering.high, nyquist * 0.999)))
    _w, h = sosfreqz(design(safe, sample_rate), worN=frequencies,
                     fs=sample_rate)
    one_pass = np.abs(h)
    # the floor keeps log10 finite where the response underflows; at
    # −240 dB the curve is off any axis a person would draw anyway
    squared = np.maximum(one_pass ** 2, 1e-12)
    return frequencies, 20.0 * np.log10(squared)