API reference

Quantity of motion

micromotion.qom

Quantity of motion.

The signature measure of this research programme: the average speed of a body part, in millimetres per second, restricted to the micromotion band.

There is one definition and three named variants. The point of naming them is that they answer different questions and previously differed only by an undocumented line of filter code. Any figure or paper should be able to say which variant produced its numbers.

Definition

Band-limit each axis to 0.3-10 Hz, bring it to velocity, band-limit again, take the Euclidean norm across axes, and report the mean in mm/s.

Acceleration is brought to velocity by integration, position by differentiation. The second band-limiting is not cosmetic. Integrating a signal with any residual offset produces a ramp that dominates the result, and differentiating amplifies the high-frequency noise that the first pass was meant to exclude.

Variants

raw The band as defined. Contains respiration, the ballistocardiac impulse and postural sway together. This is what the deposited files report unless they say otherwise. compensated Respiration and cardiac activity removed: the lower edge is raised to 0.5 Hz and the per-recording cardiac peak is notched out. What is left is postural micromotion. tilt_corrected For a single accelerometer only. A body-worn accelerometer cannot distinguish leaning from translating: tilting into gravity produces an acceleration with no displacement. Where a gyroscope is available the tilt component is estimated and removed. Measured directly on the fNIRS session, tilt inflates raw QoM by 1.56x.

G module-attribute

G = 9.80665

Standard gravity, m/s^2.

Present because accelerometers export in g and the band-limited result must be in SI before integration. Getting this wrong is not hypothetical: every phone quantity of motion in this project was 9.80665x too large until 2026-07-28.

BANDS module-attribute

BANDS = {'micromotion': filters.BAND, 'optical_legacy': filters.OPTICAL_LEGACY_BAND}

The two conventions in use, by name.

micromotion is 0.3-10 Hz and is the only one that can be applied to every sensor, so it is the one a cross-collection comparison must use. optical_legacy is the 10 Hz low-pass behind the published championship figures; it retains sub-0.3 Hz postural drift and reads about 15 per cent higher.

QomResult dataclass

Quantity of motion, with the series it was reduced from.

Source code in src/micromotion/qom.py

@dataclass
class QomResult:
    """Quantity of motion, with the series it was reduced from."""

    mean_mm_s: float
    median_mm_s: float
    speed: np.ndarray = field(repr=False)
    fs: float = 0.0
    variant: str = "raw"
    cardiac_hz: float = float("nan")
    n_samples: int = 0
    edge_samples: int = 0

    def binned(self, bin_s: float = 5.0):
        """Average speed in fixed-width bins, as a DataFrame.

        The final bin is usually partial and is flagged rather than dropped, because a
        short bin is not comparable with a full one and silently including it inflated the
        deposited five-second series three- to fourteenfold.
        """
        import pandas as pd

        idx = (np.arange(self.n_samples) / self.fs // bin_s).astype(int)
        counts = np.bincount(idx)
        means = np.bincount(idx, self.speed) / counts
        full = int(round(bin_s * self.fs))
        edge = np.where(counts < full, "partial", "ok").astype(object)
        n_edge = max(1, int(np.ceil(self.edge_samples / (bin_s * self.fs))))
        for i in list(range(n_edge)) + list(range(len(counts) - n_edge, len(counts))):
            if 0 <= i < len(counts) and edge[i] == "ok":
                edge[i] = "filter_transient"
        return pd.DataFrame(
            {
                "time_s": np.arange(len(counts)) * bin_s,
                "qom_mm_s": means,
                "n_samples": counts,
                "edge": edge,
            }
        )

binned

binned(bin_s: float = 5.0)

Average speed in fixed-width bins, as a DataFrame.

The final bin is usually partial and is flagged rather than dropped, because a short bin is not comparable with a full one and silently including it inflated the deposited five-second series three- to fourteenfold.

Source code in src/micromotion/qom.py

def binned(self, bin_s: float = 5.0):
    """Average speed in fixed-width bins, as a DataFrame.

    The final bin is usually partial and is flagged rather than dropped, because a
    short bin is not comparable with a full one and silently including it inflated the
    deposited five-second series three- to fourteenfold.
    """
    import pandas as pd

    idx = (np.arange(self.n_samples) / self.fs // bin_s).astype(int)
    counts = np.bincount(idx)
    means = np.bincount(idx, self.speed) / counts
    full = int(round(bin_s * self.fs))
    edge = np.where(counts < full, "partial", "ok").astype(object)
    n_edge = max(1, int(np.ceil(self.edge_samples / (bin_s * self.fs))))
    for i in list(range(n_edge)) + list(range(len(counts) - n_edge, len(counts))):
        if 0 <= i < len(counts) and edge[i] == "ok":
            edge[i] = "filter_transient"
    return pd.DataFrame(
        {
            "time_s": np.arange(len(counts)) * bin_s,
            "qom_mm_s": means,
            "n_samples": counts,
            "edge": edge,
        }
    )

speed_from_acceleration

speed_from_acceleration(acc, fs: float, unit: str = 'm/s^2', lo: float = filters.BAND[0], hi: float = filters.BAND[1], notch_hz: float | None = None, integrate: str = 'rectangle') -> np.ndarray

Band-limited speed, in mm/s, from acceleration.

acc is (n_samples, n_axes). unit is "m/s^2" or "g".

integrate selects the quadrature rule, and the choice is not cosmetic. Both are in use in this corpus: the StillStanding365 and fNIRS pipelines integrate with the trapezoid rule, the Stillness2025, Taqasim and 2024 championship pipelines with a rectangle sum. They differ by about 0.26 per cent on real phone data, which is small but is a systematic bias rather than noise -- the rectangle rule lags the signal by half a sample.

Neither rule is universally right here, so the default is the one that keeps this package's own published numbers self-consistent: "rectangle" is what the harmonised cross-collection table and every figure derived from it were computed with, and it reproduces the deposited Taqasim value (93.140 against 93.091 mm/s) where the trapezoid rule gives 93.405. Pass integrate="trapezoid" to reproduce StillStanding365 and fNIRS, whose deposited pipelines use it.

Which rule the project should standardise on is an open question, deliberately not settled by this default.

Source code in src/micromotion/qom.py

def speed_from_acceleration(
    acc,
    fs: float,
    unit: str = "m/s^2",
    lo: float = filters.BAND[0],
    hi: float = filters.BAND[1],
    notch_hz: float | None = None,
    integrate: str = "rectangle",
) -> np.ndarray:
    """Band-limited speed, in mm/s, from acceleration.

    ``acc`` is (n_samples, n_axes). ``unit`` is ``"m/s^2"`` or ``"g"``.

    ``integrate`` selects the quadrature rule, and the choice is not cosmetic. Both are in
    use in this corpus: the StillStanding365 and fNIRS pipelines integrate with the trapezoid
    rule, the Stillness2025, Taqasim and 2024 championship pipelines with a rectangle sum.
    They differ by about 0.26 per cent on real phone data, which is small but is a systematic
    bias rather than noise -- the rectangle rule lags the signal by half a sample.

    Neither rule is universally right here, so the default is the one that keeps this
    package's own published numbers self-consistent: ``"rectangle"`` is what the harmonised
    cross-collection table and every figure derived from it were computed with, and it
    reproduces the deposited Taqasim value (93.140 against 93.091 mm/s) where the trapezoid
    rule gives 93.405. Pass ``integrate="trapezoid"`` to reproduce StillStanding365 and
    fNIRS, whose deposited pipelines use it.

    Which rule the project should standardise on is an open question, deliberately not
    settled by this default.
    """
    a = _to_2d(acc)
    if unit == "g":
        a = a * G
    elif unit != "m/s^2":
        raise ValueError(f"unknown acceleration unit {unit!r}; use 'm/s^2' or 'g'")
    a = filters.bandpass(a, fs, lo, hi)
    if notch_hz:
        a = filters.notch(a, fs, notch_hz)
    if integrate == "trapezoid":
        v = cumulative_trapezoid(a, dx=1.0 / fs, initial=0, axis=0)
    elif integrate == "rectangle":
        v = np.cumsum(a, axis=0) / fs
    else:
        raise ValueError(f"unknown rule {integrate!r}; use 'trapezoid' or 'rectangle'")
    v = filters.bandpass(v, fs, lo, hi)
    return np.linalg.norm(v, axis=1) * MM_PER_M

derivative

derivative(x, fs: float) -> np.ndarray

Fourth-order central difference along the first axis.

The two-point central difference that numpy.gradient computes has the frequency response sin(w*dt)/dt rather than w, so it increasingly under-reads towards Nyquist. That does not matter at 200 Hz, where the top of the micromotion band is a twentieth of Nyquist, but it matters at the 20 Hz common rate where the band edge is Nyquist itself: on real optical data the two-point rule loses 4.9 per cent of the quantity of motion across that resampling, and this rule loses 2.0 per cent.

A spectral derivative would be exact for a truly band-limited signal and is not used here, because these recordings do not begin and end at the same value and the implied wraparound step adds broadband energy that differentiation then amplifies. Measured on the same file, it inflated the result by 24 per cent.

Source code in src/micromotion/qom.py

def derivative(x, fs: float) -> np.ndarray:
    """Fourth-order central difference along the first axis.

    The two-point central difference that ``numpy.gradient`` computes has the frequency
    response sin(w*dt)/dt rather than w, so it increasingly under-reads towards Nyquist.
    That does not matter at 200 Hz, where the top of the micromotion band is a twentieth of
    Nyquist, but it matters at the 20 Hz common rate where the band edge is Nyquist itself:
    on real optical data the two-point rule loses 4.9 per cent of the quantity of motion
    across that resampling, and this rule loses 2.0 per cent.

    A spectral derivative would be exact for a truly band-limited signal and is not used
    here, because these recordings do not begin and end at the same value and the implied
    wraparound step adds broadband energy that differentiation then amplifies. Measured on
    the same file, it inflated the result by 24 per cent.
    """
    x = _to_2d(x)
    if len(x) < 5:
        return np.gradient(x, 1.0 / fs, axis=0)
    v = np.empty_like(x)
    v[2:-2] = (x[:-4] - 8 * x[1:-3] + 8 * x[3:-1] - x[4:]) * (fs / 12.0)
    v[:2], v[-2:] = v[2], v[-3]
    return v

speed_from_position

speed_from_position(pos, fs: float, unit: str = 'mm', lo: float = filters.BAND[0], hi: float = filters.BAND[1]) -> np.ndarray

Band-limited speed, in mm/s, from position.

pos is (n_samples, n_axes), normally the three coordinates of one optical marker. unit is "mm" or "m".

Source code in src/micromotion/qom.py

def speed_from_position(
    pos,
    fs: float,
    unit: str = "mm",
    lo: float = filters.BAND[0],
    hi: float = filters.BAND[1],
) -> np.ndarray:
    """Band-limited speed, in mm/s, from position.

    ``pos`` is (n_samples, n_axes), normally the three coordinates of one optical marker.
    ``unit`` is ``"mm"`` or ``"m"``.
    """
    p = _to_2d(pos)
    if unit == "m":
        p = p * MM_PER_M
    elif unit != "mm":
        raise ValueError(f"unknown position unit {unit!r}; use 'mm' or 'm'")
    p = filters.bandpass(p, fs, lo, hi)
    v = derivative(p, fs)
    v = filters.bandpass(v, fs, lo, hi)
    return np.linalg.norm(v, axis=1)

qom

qom(data, fs: float, kind: str = 'acceleration', unit: str | None = None, variant: str = 'raw', band: str = 'micromotion', gyro=None, integrate: str = 'rectangle') -> QomResult

Quantity of motion for one recording.

Parameters:

Name	Type	Description	Default
data		(n_samples, n_axes) acceleration or position.	required
fs	float	Measured sampling rate. Use the rate measured from the timestamps, not the nominal one; see :func:micromotion.resample.measured_rate.	required
kind	str	"acceleration" or "position".	'acceleration'
unit	str | None	Defaults to "m/s^2" for acceleration and "mm" for position.	None
variant	str	"raw", "compensated" or "tilt_corrected".	'raw'
band	str	"micromotion" or "optical_legacy". See :data:BANDS.	'micromotion'
gyro		(n_samples, 3) angular velocity in rad/s. Required by tilt_corrected.	None

Source code in src/micromotion/qom.py

def qom(
    data,
    fs: float,
    kind: str = "acceleration",
    unit: str | None = None,
    variant: str = "raw",
    band: str = "micromotion",
    gyro=None,
    integrate: str = "rectangle",
) -> QomResult:
    """Quantity of motion for one recording.

    Parameters
    ----------
    data
        (n_samples, n_axes) acceleration or position.
    fs
        Measured sampling rate. Use the rate measured from the timestamps, not the nominal
        one; see :func:`micromotion.resample.measured_rate`.
    kind
        ``"acceleration"`` or ``"position"``.
    unit
        Defaults to ``"m/s^2"`` for acceleration and ``"mm"`` for position.
    variant
        ``"raw"``, ``"compensated"`` or ``"tilt_corrected"``.
    band
        ``"micromotion"`` or ``"optical_legacy"``. See :data:`BANDS`.
    gyro
        (n_samples, 3) angular velocity in rad/s. Required by ``tilt_corrected``.
    """
    x = _to_2d(data)
    if band not in BANDS:
        raise ValueError(f"unknown band {band!r}; use one of {sorted(BANDS)}")
    lo, hi = BANDS[band]
    hz = float("nan")

    if lo is None and kind == "acceleration":
        raise ValueError(
            "the 'optical_legacy' band has no lower edge, and an accelerometer cannot be "
            "integrated without one: gravity is a DC term and any residual offset becomes "
            "a ramp. Use band='micromotion' for accelerometer data."
        )

    if kind == "acceleration":
        unit = unit or "m/s^2"
        if variant == "compensated":
            mag = np.linalg.norm(x, axis=1)
            hz = cardiac_peak(mag, fs)
            speed = speed_from_acceleration(x, fs, unit, lo=0.5, hi=hi, notch_hz=hz,
                                            integrate=integrate)
        elif variant == "tilt_corrected":
            if gyro is None:
                raise ValueError("variant 'tilt_corrected' needs a gyroscope signal")
            speed = speed_from_acceleration(
                remove_tilt(x, gyro, fs, unit), fs, "m/s^2", lo, hi
            )
        elif variant == "raw":
            speed = speed_from_acceleration(x, fs, unit, lo, hi, integrate=integrate)
        else:
            raise ValueError(f"unknown variant {variant!r}")
    elif kind == "position":
        if variant != "raw":
            raise ValueError(
                f"variant {variant!r} applies to accelerometers only; optical position "
                "measures displacement directly, so there is no tilt ambiguity and no "
                "integration drift to compensate"
            )
        speed = speed_from_position(x, fs, unit or "mm", lo, hi)
    else:
        raise ValueError(f"unknown kind {kind!r}; use 'acceleration' or 'position'")

    return QomResult(
        mean_mm_s=float(np.mean(speed)),
        median_mm_s=float(np.median(speed)),
        speed=speed,
        fs=fs,
        variant=variant,
        cardiac_hz=hz,
        n_samples=len(speed),
        edge_samples=filters.edge_transient_samples(fs, lo),
    )

remove_tilt

remove_tilt(acc, gyro, fs: float, unit: str = 'm/s^2') -> np.ndarray

Subtract the gravity component that rotation moves between axes.

The sensor's orientation is tracked by integrating the gyroscope, the gravity vector is rotated into the sensor frame at each sample, and what remains is translation. The integration drifts, so the estimated gravity direction is high-passed back towards the measured one; this is a complementary filter, not an attitude estimator, and it is adequate only because the body barely moves.

Source code in src/micromotion/qom.py

def remove_tilt(acc, gyro, fs: float, unit: str = "m/s^2") -> np.ndarray:
    """Subtract the gravity component that rotation moves between axes.

    The sensor's orientation is tracked by integrating the gyroscope, the gravity vector is
    rotated into the sensor frame at each sample, and what remains is translation. The
    integration drifts, so the estimated gravity direction is high-passed back towards the
    measured one; this is a complementary filter, not an attitude estimator, and it is
    adequate only because the body barely moves.
    """
    a = _to_2d(acc) * (G if unit == "g" else 1.0)
    w = _to_2d(gyro)
    if w.shape[0] != a.shape[0]:
        raise ValueError("gyroscope and accelerometer must have the same length")

    g_hat = np.zeros_like(a)
    g_hat[0] = a[0] / (np.linalg.norm(a[0]) + 1e-12)
    dt = 1.0 / fs
    tau = 1.0                       # s; trust the accelerometer beyond this
    alpha = tau / (tau + dt)
    for i in range(1, len(a)):
        # rotate the previous estimate by -omega*dt (the frame turns, the vector does not)
        gyro_step = g_hat[i - 1] - np.cross(w[i], g_hat[i - 1]) * dt
        meas = a[i] / (np.linalg.norm(a[i]) + 1e-12)
        g = alpha * gyro_step + (1 - alpha) * meas
        g_hat[i] = g / (np.linalg.norm(g) + 1e-12)

    g_mag = float(np.median(np.linalg.norm(a, axis=1)))
    return a - g_hat * g_mag

tilt_fraction

tilt_fraction(acc, gyro, fs: float, unit: str = 'm/s^2') -> dict

How much of a body-worn accelerometer's quantity of motion is tilt rather than travel.

An accelerometer cannot tell leaning from moving: rotating in the gravity field produces an acceleration with no displacement. Where a gyroscope is present the rotation is known, so the gravity component it accounts for can be removed and the two compared.

Measured on the fNIRS session this returns 1.56, meaning the raw figure is a little over half again the translational one. Do not read a single session as a population value; the point is that the inflation is measurable rather than assumed, and the assumption in circulation was 1.3 to 1.5.

Source code in src/micromotion/qom.py

def tilt_fraction(acc, gyro, fs: float, unit: str = "m/s^2") -> dict:
    """How much of a body-worn accelerometer's quantity of motion is tilt rather than travel.

    An accelerometer cannot tell leaning from moving: rotating in the gravity field produces
    an acceleration with no displacement. Where a gyroscope is present the rotation is known,
    so the gravity component it accounts for can be removed and the two compared.

    Measured on the fNIRS session this returns 1.56, meaning the raw figure is a little over
    half again the translational one. Do not read a single session as a population value; the
    point is that the inflation is measurable rather than assumed, and the assumption in
    circulation was 1.3 to 1.5.
    """
    raw = qom(acc, fs, kind="acceleration", unit=unit).mean_mm_s
    corrected = qom(remove_tilt(acc, gyro, fs, unit), fs,
                    kind="acceleration", unit="m/s^2").mean_mm_s
    return {"raw_mm_s": raw, "translation_mm_s": corrected,
            "inflation": raw / corrected if corrected else float("nan"),
            "tilt_fraction": 1 - corrected / raw if raw else float("nan")}

envelope

envelope(x, fs, smooth=1.0, normalize=True)

Smooth, optionally z-scored envelope of a signal: Savitzky-Golay smoothing (order 2, window smooth seconds) followed by standardisation. Used to compare motion/audio envelopes across sources on a common, amplitude-free scale.

Source: Westney-comparisons study (Jensenius).

Args: x (np.ndarray): Input 1-D signal. fs (float): Sampling rate of the signal (Hz). smooth (float, optional): Smoothing window in seconds. None or 0 disables smoothing. Defaults to 1.0. normalize (bool, optional): If True, z-score the result. Defaults to True.

Returns: np.ndarray: The smoothed (and optionally z-scored) envelope, same length as the input.

Source code in src/micromotion/qom.py

def envelope(x, fs, smooth=1.0, normalize=True):
    """
    Smooth, optionally z-scored envelope of a signal: Savitzky-Golay
    smoothing (order 2, window `smooth` seconds) followed by
    standardisation. Used to compare motion/audio envelopes across sources
    on a common, amplitude-free scale.

    Source: Westney-comparisons study (Jensenius).

    Args:
        x (np.ndarray): Input 1-D signal.
        fs (float): Sampling rate of the signal (Hz).
        smooth (float, optional): Smoothing window in seconds. None or 0 disables
            smoothing. Defaults to 1.0.
        normalize (bool, optional): If True, z-score the result. Defaults to True.

    Returns:
        np.ndarray: The smoothed (and optionally z-scored) envelope, same length
            as the input.
    """
    from scipy.signal import savgol_filter
    x = np.asarray(x, float)
    if smooth:
        w = max(3, int(smooth * fs) | 1)
        if len(x) > w:
            x = savgol_filter(x, w, 2)
    if normalize:
        x = (x - x.mean()) / (x.std() + 1e-9)
    return x

bin_series

bin_series(x, fs, bin_s=1.0)

Mean of consecutive, non-overlapping bins of a signal (e.g. a per-second quantity-of-motion envelope from a per-frame speed series). Trailing samples that do not fill a whole bin are dropped.

Source: stillstanding study (Jensenius); also used in the Westney-comparisons study as a per-second envelope.

Args: x (np.ndarray): Input 1-D signal. fs (float): Sampling rate of the signal (Hz). bin_s (float, optional): Bin length in seconds. Defaults to 1.0.

Returns: np.ndarray: One mean value per bin (empty if the signal is shorter than two bins).

Source code in src/micromotion/qom.py

def bin_series(x, fs, bin_s=1.0):
    """
    Mean of consecutive, non-overlapping bins of a signal (e.g. a per-second
    quantity-of-motion envelope from a per-frame speed series). Trailing
    samples that do not fill a whole bin are dropped.

    Source: stillstanding study (Jensenius); also used in the
    Westney-comparisons study as a per-second envelope.

    Args:
        x (np.ndarray): Input 1-D signal.
        fs (float): Sampling rate of the signal (Hz).
        bin_s (float, optional): Bin length in seconds. Defaults to 1.0.

    Returns:
        np.ndarray: One mean value per bin (empty if the signal is shorter than
            two bins).
    """
    x = np.asarray(x, float)
    step = max(1, int(round(fs * bin_s)))
    n = len(x) // step
    if n < 2:
        return np.array([])
    return x[:n * step].reshape(n, step).mean(axis=1)

band_limited_qom

band_limited_qom(pos, fs, lo=0.3, hi=15.0, order=4, auto_decimate=True)

Band-limited quantity of motion from a position trajectory: the position is band-pass filtered (zero phase) to [lo, hi] Hz and the QoM is the per-frame speed, i.e. the Euclidean norm of the first difference times the sampling rate (units of the input per second).

For very low bands relative to the sampling rate (band edge below about fs/40), a direct high-order band-pass is numerically fragile; in that regime the trajectory is first decimated (zero phase) so the band sits comfortably in the new Nyquist range, then filtered with a second-order section (SOS) band-pass. This is the "slow sway" regime (e.g. 0.1-0.5 Hz postural sway from 100 Hz mocap). Set auto_decimate=False to force the direct filter.

Source: stillstanding study and Westney-comparisons study (Jensenius) -- this unifies the band-limited QoM cores used on mocap markers (mm), MediaPipe landmarks (px) and slow postural sway across both studies.

Args: pos (np.ndarray): Position trajectory of shape (N,) or (N, D) (e.g. D=2 image coordinates or D=3 mocap coordinates). Non-finite samples are linearly interpolated per dimension. fs (float): Sampling rate of the trajectory (Hz). lo (float, optional): Lower band edge (Hz). Defaults to 0.3. hi (float, optional): Upper band edge (Hz), clipped to 0.9 x Nyquist. Defaults to 15.0. order (int, optional): Butterworth order of the direct band-pass. Defaults to 4. auto_decimate (bool, optional): Enable the decimate+SOS low-band regime. Defaults to True.

Returns: tuple: (speed, fs_out) where speed is the per-frame speed series (length N-1, or shorter when decimated) and fs_out is its sampling rate (equal to fs unless decimated). speed is empty (and fs_out equals the input fs) when the input has fewer than int(fs) + 5 samples, or when it still contains non-finite samples after per-dimension interpolation (i.e. a dimension had fewer than 3 finite samples to interpolate from). In the auto-decimate regime, speed is also empty (with fs_out the decimated rate) when decimation leaves fewer than ~30 samples -- too few for a stable SOS band-pass.

Raises: ValueError: If the band is invalid, i.e. does not satisfy 0 < lo < hi <= 0.45*fs (after hi is clipped to 0.9 x Nyquist).

Source code in src/micromotion/qom.py

def band_limited_qom(pos, fs, lo=0.3, hi=15.0, order=4, auto_decimate=True):
    """
    Band-limited quantity of motion from a position trajectory: the position
    is band-pass filtered (zero phase) to `[lo, hi]` Hz and the QoM is the
    per-frame speed, i.e. the Euclidean norm of the first difference times
    the sampling rate (units of the input per second).

    For very low bands relative to the sampling rate (band edge below about
    fs/40), a direct high-order band-pass is numerically fragile; in that
    regime the trajectory is first decimated (zero phase) so the band sits
    comfortably in the new Nyquist range, then filtered with a second-order
    section (SOS) band-pass. This is the "slow sway" regime (e.g. 0.1-0.5 Hz
    postural sway from 100 Hz mocap). Set `auto_decimate=False` to force the
    direct filter.

    Source: stillstanding study and Westney-comparisons study (Jensenius) --
    this unifies the band-limited QoM cores used on mocap markers (mm),
    MediaPipe landmarks (px) and slow postural sway across both studies.

    Args:
        pos (np.ndarray): Position trajectory of shape (N,) or (N, D) (e.g. D=2
            image coordinates or D=3 mocap coordinates). Non-finite samples are
            linearly interpolated per dimension.
        fs (float): Sampling rate of the trajectory (Hz).
        lo (float, optional): Lower band edge (Hz). Defaults to 0.3.
        hi (float, optional): Upper band edge (Hz), clipped to 0.9 x Nyquist.
            Defaults to 15.0.
        order (int, optional): Butterworth order of the direct band-pass.
            Defaults to 4.
        auto_decimate (bool, optional): Enable the decimate+SOS low-band regime.
            Defaults to True.

    Returns:
        tuple: `(speed, fs_out)` where `speed` is the per-frame speed series
            (length N-1, or shorter when decimated) and `fs_out` is its
            sampling rate (equal to `fs` unless decimated). `speed` is empty
            (and `fs_out` equals the input `fs`) when the input has fewer
            than `int(fs) + 5` samples, or when it still contains non-finite
            samples after per-dimension interpolation (i.e. a dimension had
            fewer than 3 finite samples to interpolate from). In the
            auto-decimate regime, `speed` is also empty (with `fs_out` the
            decimated rate) when decimation leaves fewer than ~30 samples --
            too few for a stable SOS band-pass.

    Raises:
        ValueError: If the band is invalid, i.e. does not satisfy
            `0 < lo < hi <= 0.45*fs` (after `hi` is clipped to 0.9 x Nyquist).
    """
    from scipy import signal
    pos = np.asarray(pos, float)
    if pos.ndim == 1:
        pos = pos[:, None]
    pos = np.column_stack([_interp_nans(pos[:, i]) for i in range(pos.shape[1])])

    hi_eff = min(hi, 0.9 * fs / 2)
    if not (0 < lo < hi_eff <= 0.45 * fs + 1e-9):
        raise ValueError("band must satisfy 0 < lo < hi <= 0.45*fs")

    if len(pos) < int(fs) + 5 or not np.isfinite(pos).all():
        return np.array([]), fs

    if auto_decimate and fs / hi_eff >= 40:
        q = int(min(13, fs // (20 * hi_eff)))
        pos = signal.decimate(pos, q, axis=0, zero_phase=True)
        fs_out = fs / q
        if len(pos) < 30:
            return np.array([]), fs_out
        sos = signal.butter(2, [lo / (fs_out / 2), hi_eff / (fs_out / 2)],
                            btype="band", output="sos")
        filtered = signal.sosfiltfilt(sos, pos, axis=0)
    else:
        fs_out = fs
        b, a = signal.butter(order, [lo / (fs / 2), hi_eff / (fs / 2)], btype="band")
        filtered = signal.filtfilt(b, a, pos, axis=0)
    speed = np.linalg.norm(np.diff(filtered, axis=0), axis=1) * fs_out
    return speed, fs_out

accel_to_speed

accel_to_speed(acc, fs, highpass=0.3, order=2, normalize_gravity=False)

Integrated speed from a 3-axis accelerometer: each axis is high-pass filtered (removing gravity and DC), integrated to velocity, high-pass filtered again (killing integration drift), and the speed is the Euclidean norm of the velocity (m/s for input in m/s^2).

Source: stillstanding study (Jensenius) -- the "corpus method" for integrated quantity of motion from chest-worn accelerometers.

Args: acc (np.ndarray): Acceleration of shape (N, 3) in m/s^2 (or raw counts with normalize_gravity=True). fs (float): Sampling rate (Hz). highpass (float, optional): High-pass cutoff (Hz) used both before and after integration. Defaults to 0.3. order (int, optional): Butterworth order of the high-pass filters. Defaults to 2. normalize_gravity (bool, optional): If True, rescale the raw input so that the median vector magnitude equals 1 g (9.80665 m/s^2) before filtering -- useful for uncalibrated sensors whose resting output should be gravity. Defaults to False.

Returns: np.ndarray: Speed series of length N (m/s).

Source code in src/micromotion/qom.py

def accel_to_speed(acc, fs, highpass=0.3, order=2, normalize_gravity=False):
    """
    Integrated speed from a 3-axis accelerometer: each axis is high-pass
    filtered (removing gravity and DC), integrated to velocity, high-pass
    filtered again (killing integration drift), and the speed is the
    Euclidean norm of the velocity (m/s for input in m/s^2).

    Source: stillstanding study (Jensenius) -- the "corpus method" for
    integrated quantity of motion from chest-worn accelerometers.

    Args:
        acc (np.ndarray): Acceleration of shape (N, 3) in m/s^2 (or raw counts
            with `normalize_gravity=True`).
        fs (float): Sampling rate (Hz).
        highpass (float, optional): High-pass cutoff (Hz) used both before and
            after integration. Defaults to 0.3.
        order (int, optional): Butterworth order of the high-pass filters.
            Defaults to 2.
        normalize_gravity (bool, optional): If True, rescale the raw input so
            that the median vector magnitude equals 1 g (9.80665 m/s^2) before
            filtering -- useful for uncalibrated sensors whose resting output
            should be gravity. Defaults to False.

    Returns:
        np.ndarray: Speed series of length N (m/s).
    """
    from scipy.signal import butter, filtfilt
    G = 9.80665
    acc = np.asarray(acc, float)
    if normalize_gravity:
        norm = np.linalg.norm(acc, axis=1)
        acc = acc / np.median(norm) * G
    b, a = butter(order, highpass / (fs / 2), btype="high")
    acc = filtfilt(b, a, acc, axis=0)
    vel = np.cumsum(acc, axis=0) / fs
    vel = filtfilt(b, a, vel, axis=0)
    return np.linalg.norm(vel, axis=1)

group_qom

group_qom(points, fs, lo=0.3, hi=15.0, **kwargs)

Mean band-limited quantity of motion over a group of markers/landmarks, plus the group's mean speed envelope: each trajectory is passed through band_limited_qom and the per-trajectory speeds are averaged.

Source: stillstanding study and Westney-comparisons study (Jensenius) -- per-body-part QoM (head, shoulders, arms, wrists) from mocap markers and pose landmarks.

Args: points (np.ndarray): Trajectories of shape (N, M, D): N frames, M markers/landmarks, D spatial dimensions. fs (float): Sampling rate (Hz). lo (float, optional): Lower band edge (Hz). Defaults to 0.3. hi (float, optional): Upper band edge (Hz). Defaults to 15.0. **kwargs: Passed on to band_limited_qom.

Returns: tuple: (qom, speed, fs_out) where qom is the mean speed across markers and time (NaN if no marker yields a valid series), speed is the group's mean per-frame speed series, and fs_out its sampling rate.

Source code in src/micromotion/qom.py

def group_qom(points, fs, lo=0.3, hi=15.0, **kwargs):
    """
    Mean band-limited quantity of motion over a group of markers/landmarks,
    plus the group's mean speed envelope: each trajectory is passed through
    `band_limited_qom` and the per-trajectory speeds are averaged.

    Source: stillstanding study and Westney-comparisons study (Jensenius) --
    per-body-part QoM (head, shoulders, arms, wrists) from mocap markers and
    pose landmarks.

    Args:
        points (np.ndarray): Trajectories of shape (N, M, D): N frames, M
            markers/landmarks, D spatial dimensions.
        fs (float): Sampling rate (Hz).
        lo (float, optional): Lower band edge (Hz). Defaults to 0.3.
        hi (float, optional): Upper band edge (Hz). Defaults to 15.0.
        **kwargs: Passed on to `band_limited_qom`.

    Returns:
        tuple: `(qom, speed, fs_out)` where `qom` is the mean speed across
            markers and time (NaN if no marker yields a valid series), `speed`
            is the group's mean per-frame speed series, and `fs_out` its
            sampling rate.
    """
    points = np.asarray(points, float)
    speeds, fs_out = [], fs
    for m in range(points.shape[1]):
        sp, fs_out = band_limited_qom(points[:, m, :], fs, lo=lo, hi=hi, **kwargs)
        if len(sp):
            speeds.append(sp)
    if not speeds:
        return np.nan, np.array([]), fs_out
    L = min(len(s) for s in speeds)
    mean_speed = np.mean([s[:L] for s in speeds], axis=0)
    return float(np.mean([s.mean() for s in speeds])), mean_speed, fs_out

pose_qom

pose_qom(landmarks, fs, lo=0.3, hi=5.0, **kwargs)

Band-limited quantity of motion of 2-D pose landmarks (px/s): a thin wrapper around group_qom with the band used for image-space pose trajectories (0.3-5 Hz), where higher bands are dominated by landmark jitter rather than motion.

Source: Westney-comparisons study (Jensenius).

Args: landmarks (np.ndarray): Landmark trajectories of shape (N, L, 2) in pixels (a single landmark of shape (N, 2) is also accepted). fs (float): Sampling rate (Hz, e.g. video frame rate). lo (float, optional): Lower band edge (Hz). Defaults to 0.3. hi (float, optional): Upper band edge (Hz). Defaults to 5.0. **kwargs: Passed on to band_limited_qom.

Returns: tuple: (qom, speed, fs_out) as in group_qom.

Source code in src/micromotion/qom.py

def pose_qom(landmarks, fs, lo=0.3, hi=5.0, **kwargs):
    """
    Band-limited quantity of motion of 2-D pose landmarks (px/s): a thin
    wrapper around `group_qom` with the band used for image-space pose
    trajectories (0.3-5 Hz), where higher bands are dominated by landmark
    jitter rather than motion.

    Source: Westney-comparisons study (Jensenius).

    Args:
        landmarks (np.ndarray): Landmark trajectories of shape (N, L, 2) in
            pixels (a single landmark of shape (N, 2) is also accepted).
        fs (float): Sampling rate (Hz, e.g. video frame rate).
        lo (float, optional): Lower band edge (Hz). Defaults to 0.3.
        hi (float, optional): Upper band edge (Hz). Defaults to 5.0.
        **kwargs: Passed on to `band_limited_qom`.

    Returns:
        tuple: `(qom, speed, fs_out)` as in `group_qom`.
    """
    landmarks = np.asarray(landmarks, float)
    if landmarks.ndim == 2:
        landmarks = landmarks[:, None, :]
    return group_qom(landmarks, fs, lo=lo, hi=hi, **kwargs)

body_scale

body_scale(landmarks, upper=(11, 12), lower=(23, 24))

Body-size scale (in the landmarks' own units, e.g. pixels) as the median torso length: the distance from the midpoint of the upper landmarks (shoulders) to the midpoint of the lower landmarks (hips). The torso length is preferred over shoulder width because it stays robust in a profile view, where the shoulder width collapses.

The default indices are MediaPipe Pose landmarks (11/12 shoulders, 23/24 hips).

Source: Westney-comparisons study (Jensenius).

Args: landmarks (np.ndarray): Landmark trajectories of shape (N, L, C) with C >= 2; only the first two coordinates are used. upper (tuple, optional): Indices of the two shoulder landmarks. Defaults to (11, 12). lower (tuple, optional): Indices of the two hip landmarks. Defaults to (23, 24).

Returns: float: Median torso length (NaN if no finite frames).

Source code in src/micromotion/qom.py

def body_scale(landmarks, upper=(11, 12), lower=(23, 24)):
    """
    Body-size scale (in the landmarks' own units, e.g. pixels) as the median
    torso length: the distance from the midpoint of the `upper` landmarks
    (shoulders) to the midpoint of the `lower` landmarks (hips). The torso
    length is preferred over shoulder width because it stays robust in a
    profile view, where the shoulder width collapses.

    The default indices are MediaPipe Pose landmarks (11/12 shoulders,
    23/24 hips).

    Source: Westney-comparisons study (Jensenius).

    Args:
        landmarks (np.ndarray): Landmark trajectories of shape (N, L, C) with
            C >= 2; only the first two coordinates are used.
        upper (tuple, optional): Indices of the two shoulder landmarks.
            Defaults to (11, 12).
        lower (tuple, optional): Indices of the two hip landmarks.
            Defaults to (23, 24).

    Returns:
        float: Median torso length (NaN if no finite frames).
    """
    landmarks = np.asarray(landmarks, float)
    um = (landmarks[:, upper[0], :2] + landmarks[:, upper[1], :2]) / 2
    lm = (landmarks[:, lower[0], :2] + landmarks[:, lower[1], :2]) / 2
    d = np.linalg.norm(um - lm, axis=1)
    d = d[np.isfinite(d)]
    return float(np.median(d)) if len(d) else np.nan

normalized_qom

normalized_qom(landmarks, fs, scale=None, lo=0.3, hi=5.0, upper=(11, 12), lower=(23, 24), **kwargs)

Body-scale-normalised quantity of motion (body-lengths per second): the pose QoM divided by the performer's own body scale (median torso length, see body_scale). Being dimensionless, this is invariant to camera framing/zoom and comparable across recordings.

Source: Westney-comparisons study (Jensenius) -- framing-invariant with/without-audience comparison of a pianist's motion.

Args: landmarks (np.ndarray): Landmark trajectories of shape (N, L, 2). fs (float): Sampling rate (Hz). scale (float, optional): Precomputed body scale. Defaults to None (which computes body_scale(landmarks, upper, lower)). lo (float, optional): Lower band edge (Hz). Defaults to 0.3. hi (float, optional): Upper band edge (Hz). Defaults to 5.0. upper (tuple, optional): Shoulder landmark indices for body_scale. Defaults to (11, 12). lower (tuple, optional): Hip landmark indices for body_scale. Defaults to (23, 24). **kwargs: Passed on to band_limited_qom.

Returns: tuple: (qom, speed, fs_out) as in group_qom, with both qom and speed divided by the body scale. When scale is non-finite (e.g. body_scale found no finite torso-length sample) or not strictly positive (degenerate, coincident upper/lower landmarks), division would otherwise silently propagate NaN/inf through qom and speed; instead both are explicitly returned as NaN (qom as a NaN scalar, speed as an all-NaN array of the same shape) so the invalid-scale case is unambiguous rather than merely inferred from the arithmetic.

Source code in src/micromotion/qom.py

def normalized_qom(landmarks, fs, scale=None, lo=0.3, hi=5.0,
                   upper=(11, 12), lower=(23, 24), **kwargs):
    """
    Body-scale-normalised quantity of motion (body-lengths per second):
    the pose QoM divided by the performer's own body scale (median torso
    length, see `body_scale`). Being dimensionless, this is invariant to
    camera framing/zoom and comparable across recordings.

    Source: Westney-comparisons study (Jensenius) -- framing-invariant
    with/without-audience comparison of a pianist's motion.

    Args:
        landmarks (np.ndarray): Landmark trajectories of shape (N, L, 2).
        fs (float): Sampling rate (Hz).
        scale (float, optional): Precomputed body scale. Defaults to None (which
            computes `body_scale(landmarks, upper, lower)`).
        lo (float, optional): Lower band edge (Hz). Defaults to 0.3.
        hi (float, optional): Upper band edge (Hz). Defaults to 5.0.
        upper (tuple, optional): Shoulder landmark indices for `body_scale`. Defaults to (11, 12).
        lower (tuple, optional): Hip landmark indices for `body_scale`. Defaults to (23, 24).
        **kwargs: Passed on to `band_limited_qom`.

    Returns:
        tuple: `(qom, speed, fs_out)` as in `group_qom`, with both `qom` and
            `speed` divided by the body scale. When `scale` is non-finite
            (e.g. `body_scale` found no finite torso-length sample) or not
            strictly positive (degenerate, coincident upper/lower landmarks),
            division would otherwise silently propagate NaN/inf through
            `qom` and `speed`; instead both are explicitly returned as NaN
            (`qom` as a NaN scalar, `speed` as an all-NaN array of the same
            shape) so the invalid-scale case is unambiguous rather than
            merely inferred from the arithmetic.
    """
    landmarks = np.asarray(landmarks, float)
    if scale is None:
        scale = body_scale(landmarks, upper=upper, lower=lower)
    qom, speed, fs_out = pose_qom(landmarks, fs, lo=lo, hi=hi, **kwargs)
    if not np.isfinite(scale) or scale <= 0:
        return float("nan"), np.full_like(speed, np.nan), fs_out
    return qom / scale, speed / scale, fs_out

grid_qom

grid_qom(frames, grid=(6, 4), region=(0.0, 1.0, 0.0, 1.0), threshold=8.0)

Spatial grid quantity of motion from a stack of grayscale frames: the absolute inter-frame difference is thresholded (small differences set to zero to suppress sensor noise) and averaged within each cell of a grid[0] x grid[1] grid laid over region, yielding one motion time series per cell plus a per-cell mean-motion heatmap.

Source: Westney-comparisons study (Jensenius) -- audience-region motion mapping in a concert hall.

Args: frames (np.ndarray): Grayscale frames of shape (T, H, W). grid (tuple, optional): Grid size (columns, rows). Defaults to (6, 4). region (tuple, optional): Region of interest as fractions (x0, x1, y0, y1) of the frame. Defaults to the full frame. threshold (float, optional): Absolute-difference threshold below which pixel changes are zeroed (0-255 scale). Defaults to 8.0.

Returns: tuple: (series, heat) where series has shape (T-1, rows*cols) (cells in row-major order) and heat has shape (rows, cols) with each cell's time-mean motion.

Raises: ValueError: If frames is not 3-D (T, H, W), as in _motionanalysis.motiongram_data.

Source code in src/micromotion/qom.py

def grid_qom(frames, grid=(6, 4), region=(0.0, 1.0, 0.0, 1.0), threshold=8.0):
    """
    Spatial grid quantity of motion from a stack of grayscale frames: the
    absolute inter-frame difference is thresholded (small differences set to
    zero to suppress sensor noise) and averaged within each cell of a
    `grid[0]` x `grid[1]` grid laid over `region`, yielding one motion time
    series per cell plus a per-cell mean-motion heatmap.

    Source: Westney-comparisons study (Jensenius) -- audience-region motion
    mapping in a concert hall.

    Args:
        frames (np.ndarray): Grayscale frames of shape (T, H, W).
        grid (tuple, optional): Grid size (columns, rows). Defaults to (6, 4).
        region (tuple, optional): Region of interest as fractions
            (x0, x1, y0, y1) of the frame. Defaults to the full frame.
        threshold (float, optional): Absolute-difference threshold below which
            pixel changes are zeroed (0-255 scale). Defaults to 8.0.

    Returns:
        tuple: `(series, heat)` where `series` has shape (T-1, rows*cols)
            (cells in row-major order) and `heat` has shape (rows, cols) with
            each cell's time-mean motion.

    Raises:
        ValueError: If `frames` is not 3-D (T, H, W), as in
            `_motionanalysis.motiongram_data`.
    """
    frames = np.asarray(frames, dtype=np.float32)
    if frames.ndim != 3:
        raise ValueError("grid_qom expects frames of shape (T, H, W)")
    T, H, W = frames.shape
    gx, gy = grid
    x0, x1, y0, y1 = region
    xs = np.linspace(int(x0 * W), int(x1 * W), gx + 1).astype(int)
    ys = np.linspace(int(y0 * H), int(y1 * H), gy + 1).astype(int)
    d = np.abs(np.diff(frames, axis=0))
    d[d < threshold] = 0.0
    series = np.empty((T - 1, gy * gx), dtype=np.float32)
    for r in range(gy):
        for c in range(gx):
            cell = d[:, ys[r]:ys[r + 1], xs[c]:xs[c + 1]]
            series[:, r * gx + c] = cell.mean(axis=(1, 2))
    heat = series.mean(axis=0).reshape(gy, gx)
    return series, heat

Filters and bands

micromotion.filters

Band-limiting for the micromotion band.

One definition, used everywhere. Across the Still Standing repository 46 scripts defined their own filter and they did not all agree; the differences were invisible in the output and moved quantity of motion by up to 10 per cent.

The canonical band is 0.3-10 Hz, a zero-phase Butterworth of order 4 applied as second-order sections. The lower edge sits below the respiratory rate and above the postural drift that integration turns into a ramp; the upper edge is where accelerometer noise begins to dominate the signal of a person standing still.

BAND module-attribute

BAND = (0.2, 10.0)

The micromotion band, in Hz.

The lower edge was 0.3 Hz until 2026-07-29, inherited rather than chosen. Swept across seven optical datasets, 665 recordings, the between-dataset spread is 3.2 per cent at 0.15 Hz, 2.1 at 0.20, 2.7 at 0.25, 6.2 at 0.30 and 10.1 at 0.40. 0.2 Hz is a clear optimum, and it is also where the 20 Hz origin dataset stops being an outlier in either direction.

A lower edge this close to DC has to be checked against integration drift, since that is what the edge is there to control. It survives: on accelerometer data the ratio of mean to median speed, which rises if drift is leaking in, is 2.07 at 0.2 Hz against 2.00 at 0.3 -- flat. What does change is the absolute value, by about 73 per cent, because more of the low-frequency content is retained.

See deposit/_analysis/osdb_qom/REPORT.md.

OPTICAL_LEGACY_BAND module-attribute

OPTICAL_LEGACY_BAND = (None, 10.0)

A 10 Hz low-pass with no lower edge: the convention behind the published championship quantity of motion.

It is kept because those numbers are in print, not because it is interchangeable with :data:BAND. Optical position is an absolute measurement, so sub-0.3 Hz postural drift in it is real movement and there is no reason to discard it. A body-worn accelerometer cannot offer the same choice: gravity is a DC term, and integrating any residual offset produces a ramp that swamps the result. So the lower edge is optional for position and mandatory for acceleration.

The two are not the same measure. On the 2015 championship the band-pass reads 15.5 per cent below the low-pass. Any table that puts optical and accelerometer collections side by side must therefore use :data:BAND throughout, and say so.

NARROW_BAND_RATIO module-attribute

NARROW_BAND_RATIO = 40.0

Warn when the sampling rate exceeds the upper band edge by more than this.

A band-pass whose edges sit very close to zero in normalised frequency is numerically fragile, and how fragile depends on how it is realised. Second-order sections, which this module uses, stay accurate far longer than the transfer-function form -- but not forever, and a caller designing their own filter should be warned before they are bitten.

The failure is silent and it is not hypothetical. A 0.1-0.5 Hz third-order band-pass at 250 Hz, written in the usual butter(3, [lo/ny, hi/ny]) transfer-function form, has a largest pole radius of 0.9979 and a measured passband gain of 0.84 at 0.15 Hz where it should be 0.99. Nothing raises, nothing looks wrong, and every amplitude downstream is 16 per cent low. The same design as second-order sections gives 0.9875.

The fix when this warns is to decimate first, so the band sits comfortably inside the new Nyquist range, then filter.

NYQUIST_MARGIN module-attribute

NYQUIST_MARGIN = 0.99

How close to Nyquist an upper band edge may sit, as a fraction.

A filter designed right at Nyquist has no transition band left, so the edge is pulled in. The default is conservative. The Still Standing corpus uses 0.999 throughout, which matters whenever the band edge is already near Nyquist -- a 10 Hz low-pass on 20 Hz data becomes 9.9 Hz here and 9.99 Hz there, and on real optical data that moved quantity of motion by 7e-5. Pass margin to match whatever convention the surrounding analysis uses.

bandpass

bandpass(x, fs: float, lo: float | None = BAND[0], hi: float = BAND[1], order: int = ORDER, margin: float = NYQUIST_MARGIN)

Zero-phase band-limiting along the first axis.

lo=None gives a pure low-pass, which is the :data:OPTICAL_LEGACY_BAND convention.

Source code in src/micromotion/filters.py

def bandpass(x, fs: float, lo: float | None = BAND[0], hi: float = BAND[1],
             order: int = ORDER, margin: float = NYQUIST_MARGIN):
    """Zero-phase band-limiting along the first axis.

    ``lo=None`` gives a pure low-pass, which is the :data:`OPTICAL_LEGACY_BAND` convention.
    """
    if lo is None:
        return lowpass(x, fs, hi, order, margin)
    wl, wh = _edges(fs, lo, hi, margin)
    sos = signal.butter(order, [wl, wh], btype="band", output="sos")
    return signal.sosfiltfilt(sos, np.asarray(x, float), axis=0)

lowpass

lowpass(x, fs: float, fc: float = BAND[1], order: int = ORDER, margin: float = NYQUIST_MARGIN)

Zero-phase low-pass along the first axis.

margin is how close to Nyquist fc may sit; see :data:NYQUIST_MARGIN.

Source code in src/micromotion/filters.py

def lowpass(x, fs: float, fc: float = BAND[1], order: int = ORDER,
            margin: float = NYQUIST_MARGIN):
    """Zero-phase low-pass along the first axis.

    ``margin`` is how close to Nyquist ``fc`` may sit; see :data:`NYQUIST_MARGIN`.
    """
    ny = fs / 2.0
    fc = min(fc, ny * margin)
    sos = signal.butter(order, fc / ny, btype="low", output="sos")
    return signal.sosfiltfilt(sos, np.asarray(x, float), axis=0)

highpass

highpass(x, fs: float, fc: float = BAND[0], order: int = ORDER)

Zero-phase high-pass along the first axis.

Used where the upper edge is meaningless because the rate is already near the band limit, and for gravity removal when no low-pass is wanted.

Source code in src/micromotion/filters.py

def highpass(x, fs: float, fc: float = BAND[0], order: int = ORDER):
    """Zero-phase high-pass along the first axis.

    Used where the upper edge is meaningless because the rate is already near the band
    limit, and for gravity removal when no low-pass is wanted.
    """
    ny = fs / 2.0
    if fc <= 0 or fc >= ny:
        raise ValueError(f"cutoff {fc} Hz is not below Nyquist {ny} Hz")
    sos = signal.butter(order, fc / ny, btype="high", output="sos")
    return signal.sosfiltfilt(sos, np.asarray(x, float), axis=0)

notch

notch(x, fs: float, f0: float, q: float = 6.0)

Zero-phase notch at f0 Hz.

Used to remove the cardiac peak when isolating postural micromotion. f0 is normally found with :func:micromotion.spectral.cardiac_peak.

Source code in src/micromotion/filters.py

def notch(x, fs: float, f0: float, q: float = 6.0):
    """Zero-phase notch at ``f0`` Hz.

    Used to remove the cardiac peak when isolating postural micromotion. ``f0`` is
    normally found with :func:`micromotion.spectral.cardiac_peak`.
    """
    if not np.isfinite(f0) or f0 <= 0 or f0 >= fs / 2:
        return np.asarray(x, float)
    b, a = signal.iirnotch(f0 / (fs / 2), Q=q)
    return signal.filtfilt(b, a, np.asarray(x, float), axis=0)

edge_transient_samples

edge_transient_samples(fs: float, lo: float | None = BAND[0], order: int = ORDER) -> int

Samples at each end that filtfilt contaminates.

A conservative estimate: the impulse response of the low edge, doubled for the forward-backward pass. Callers should either trim this or flag it, as the deposited five-second binning does with its edge column.

Source code in src/micromotion/filters.py

def edge_transient_samples(fs: float, lo: float | None = BAND[0], order: int = ORDER) -> int:
    """Samples at each end that ``filtfilt`` contaminates.

    A conservative estimate: the impulse response of the low edge, doubled for the
    forward-backward pass. Callers should either trim this or flag it, as the deposited
    five-second binning does with its ``edge`` column.
    """
    return int(np.ceil(2 * order * fs / (lo if lo else BAND[1])))

Sampling rates and resampling

micromotion.resample

Rate measurement and resampling.

Two rules, both learned the hard way.

Measure the rate, do not read it. Nominal rates in this corpus are wrong by up to 4.4 per cent, and one record's documented rate was out by a factor of 37.

Downsample, never upsample. Upsampling invents structure between samples, and every method that reads across scales treats the invention as real. On 2026-07-28 an analysis that upsampled 20 Hz data to 25 Hz produced multifractal widths up to 6.6 where the plausible range is around 1. Nothing failed; the numbers were simply wrong.

COMMON_RATE module-attribute

COMMON_RATE = 20.0

Hz. The cross-collection analysis rate.

Not a compromise: the micromotion band stops at 10 Hz, so 20 Hz is its Nyquist rate and everything above it is discarded by the band-pass anyway. It is also the native rate of the slowest collection, so reaching it never requires upsampling. See deposit/SAMPLING_RATES.md.

measured_rate

measured_rate(t) -> float

Sampling rate in Hz, from a timestamp vector.

Sample count over elapsed span, deliberately not the reciprocal of the median interval. Where timestamps are rounded to whole milliseconds the intervals become a mixture of adjacent integers and their median is a quantisation artefact: on Stillness2025 that route returns exactly 250 Hz for a recording that runs at 256, and on the pre-study phones it returned 636 Hz for a stream arriving at 106.

Source code in src/micromotion/resample.py

def measured_rate(t) -> float:
    """Sampling rate in Hz, from a timestamp vector.

    Sample count over elapsed span, deliberately not the reciprocal of the median interval.
    Where timestamps are rounded to whole milliseconds the intervals become a mixture of
    adjacent integers and their median is a quantisation artefact: on Stillness2025 that
    route returns exactly 250 Hz for a recording that runs at 256, and on the pre-study
    phones it returned 636 Hz for a stream arriving at 106.
    """
    t = np.asarray(t, float)
    if len(t) < 2:
        raise ValueError("need at least two timestamps")
    span = t[-1] - t[0]
    if span <= 0:
        raise ValueError("timestamps do not increase")
    return (len(t) - 1) / span

rate_quality

rate_quality(t) -> dict

How regular a timestamp vector actually is.

Returns the measured rate, the jitter, the largest gap, and counts of duplicated and backward timestamps. The balance-board files carry 4.7 per cent duplicates and 83 backward steps, and one pre-study file is 19 per cent covered because of a single 132-second gap; neither is visible from the rate alone.

Source code in src/micromotion/resample.py

def rate_quality(t) -> dict:
    """How regular a timestamp vector actually is.

    Returns the measured rate, the jitter, the largest gap, and counts of duplicated and
    backward timestamps. The balance-board files carry 4.7 per cent duplicates and 83
    backward steps, and one pre-study file is 19 per cent covered because of a single
    132-second gap; neither is visible from the rate alone.
    """
    t = np.asarray(t, float)
    d = np.diff(t)
    fs = measured_rate(t)
    nominal = 1.0 / fs
    return {
        "measured_rate_hz": fs,
        "n_samples": len(t),
        "duration_s": float(t[-1] - t[0]),
        "jitter_cv": float(np.std(d) / np.mean(d)) if np.mean(d) else float("nan"),
        "max_gap_s": float(d.max()) if len(d) else float("nan"),
        "gap_ratio": float(d.max() / nominal) if len(d) and nominal else float("nan"),
        "n_duplicate_t": int((d == 0).sum()),
        "n_backward_t": int((d < 0).sum()),
        "coverage": float(d[d <= 5 * nominal].sum() / (t[-1] - t[0])) if len(d) else 1.0,
    }

to_rate

to_rate(x, fs_in: float, fs_out: float = COMMON_RATE)

Anti-alias resample to fs_out, refusing to upsample.

Raises rather than upsampling. If a series genuinely cannot reach the target it does not belong in a comparison built at that rate, and interpolating it in would corrupt the comparison silently.

Source code in src/micromotion/resample.py

def to_rate(x, fs_in: float, fs_out: float = COMMON_RATE):
    """Anti-alias resample to ``fs_out``, refusing to upsample.

    Raises rather than upsampling. If a series genuinely cannot reach the target it does
    not belong in a comparison built at that rate, and interpolating it in would corrupt
    the comparison silently.
    """
    if fs_out > fs_in + 1e-9:
        raise ValueError(
            f"refusing to upsample {fs_in:.4g} Hz to {fs_out:.4g} Hz. "
            "Upsampling invents structure between samples; exclude this series instead."
        )
    x = np.asarray(x, float)
    if abs(fs_in - fs_out) < 1e-9:
        return x

    # Polyphase FIR, not signal.resample. The FFT resampler treats the series as periodic,
    # so on optical position data -- which carries a large offset and does not begin and end
    # at the same value -- it rings across the whole recording. That moved quantity of
    # motion by up to 8 per cent and did so non-monotonically in the rate, which is exactly
    # the kind of silent error a harmonised table cannot carry.
    up, down = Fraction(fs_out / fs_in).limit_denominator(1000).as_integer_ratio()
    if up == 0:
        raise ValueError(f"rate ratio {fs_out}/{fs_in} is too extreme to resample")
    return signal.resample_poly(x, up, down, axis=0, padtype="line")

regularize

regularize(t, x, fs_out: float | None = None, max_gap_s: float | None = None)

Put an irregularly sampled signal onto a uniform grid.

Sorts, drops duplicate and backward timestamps, then interpolates linearly. Samples that fall inside a gap longer than max_gap_s are returned as NaN rather than bridged, so that a 132-second hole cannot be mistaken for 132 seconds of stillness.

Written for the Wii balance board, whose timestamps are unsorted, 4.7 per cent duplicated, and arrive at a median 61.8 Hz that varies between files.

Source code in src/micromotion/resample.py

def regularize(t, x, fs_out: float | None = None, max_gap_s: float | None = None):
    """Put an irregularly sampled signal onto a uniform grid.

    Sorts, drops duplicate and backward timestamps, then interpolates linearly. Samples
    that fall inside a gap longer than ``max_gap_s`` are returned as NaN rather than
    bridged, so that a 132-second hole cannot be mistaken for 132 seconds of stillness.

    Written for the Wii balance board, whose timestamps are unsorted, 4.7 per cent
    duplicated, and arrive at a median 61.8 Hz that varies between files.
    """
    t = np.asarray(t, float)
    x = np.asarray(x, float)
    x2 = x[:, None] if x.ndim == 1 else x

    order = np.argsort(t, kind="stable")
    t, x2 = t[order], x2[order]
    keep = np.concatenate([[True], np.diff(t) > 0])
    t, x2 = t[keep], x2[keep]
    if len(t) < 2:
        raise ValueError("fewer than two usable timestamps after cleaning")

    fs_out = fs_out or measured_rate(t)
    grid = np.arange(t[0], t[-1], 1.0 / fs_out)
    out = np.column_stack([np.interp(grid, t, x2[:, i]) for i in range(x2.shape[1])])

    if max_gap_s is not None:
        d = np.diff(t)
        for i in np.flatnonzero(d > max_gap_s):
            out[(grid > t[i]) & (grid < t[i + 1])] = np.nan

    return grid, (out[:, 0] if x.ndim == 1 else out)

interpolate_gaps

interpolate_gaps(x, max_gap: int = 200)

Bridge short runs of NaN, leave long ones alone.

Filters cannot run across missing samples, so gaps have to be handled before anything else. Bridging a dropped frame is reconstruction; bridging a 469-second hole is invention, and the difference is only a matter of degree, which is why the threshold is explicit and the long gaps stay NaN for the caller to exclude.

Works column by column. Leading and trailing gaps are never filled, since there is nothing on one side to interpolate from.

Source code in src/micromotion/resample.py

def interpolate_gaps(x, max_gap: int = 200):
    """Bridge short runs of NaN, leave long ones alone.

    Filters cannot run across missing samples, so gaps have to be handled before anything
    else. Bridging a dropped frame is reconstruction; bridging a 469-second hole is
    invention, and the difference is only a matter of degree, which is why the threshold is
    explicit and the long gaps stay NaN for the caller to exclude.

    Works column by column. Leading and trailing gaps are never filled, since there is
    nothing on one side to interpolate from.
    """
    a = np.array(x, float)
    one_d = a.ndim == 1
    if one_d:
        a = a[:, None]
    for j in range(a.shape[1]):
        v = a[:, j]
        bad = np.isnan(v)
        if not bad.any() or bad.all():
            continue
        good = np.flatnonzero(~bad)
        idx = np.flatnonzero(bad)
        for run in np.split(idx, np.flatnonzero(np.diff(idx) != 1) + 1):
            if len(run) <= max_gap and run[0] > 0 and run[-1] < len(v) - 1:
                v[run] = np.interp(run, good, v[good])
    return a[:, 0] if one_d else a

gap_report

gap_report(x, fs: float) -> dict

Where the missing data is, and how it is distributed.

A single missing fraction hides the distinction that matters: one per cent scattered evenly is a usable recording, and one per cent in a single block in the middle is two recordings.

Source code in src/micromotion/resample.py

def gap_report(x, fs: float) -> dict:
    """Where the missing data is, and how it is distributed.

    A single missing fraction hides the distinction that matters: one per cent scattered
    evenly is a usable recording, and one per cent in a single block in the middle is two
    recordings.
    """
    a = np.asarray(x, float)
    bad = np.isnan(a)
    if a.ndim > 1:
        bad = bad.any(axis=1)
    if not bad.any():
        return {"missing_frac": 0.0, "n_gaps": 0, "longest_gap_s": 0.0}
    idx = np.flatnonzero(bad)
    runs = np.split(idx, np.flatnonzero(np.diff(idx) != 1) + 1)
    return {
        "missing_frac": float(bad.mean()),
        "n_gaps": len(runs),
        "longest_gap_s": float(max(len(r) for r in runs) / fs),
        "median_gap_s": float(np.median([len(r) for r in runs]) / fs),
    }

Reading files

micromotion.io

Readers for the layouts this corpus actually uses.

Each returns a :class:~micromotion.record.MotionRecord. Gap sentinels become NaN, units are recorded rather than assumed, and the sampling rate is measured wherever the file carries a timebase.

:func:read dispatches on content, not on the extension, because the extension lies: the balance-board dumps are headerless and space-delimited whatever they are called, and the Qualisys family puts three different header shapes behind one name.

Y_UP_COLLECTIONS module-attribute

Y_UP_COLLECTIONS = ('Standstill2019',)

Editions whose vertical axis is Y.

The 2019 championship was recorded on OptiTrack and converted, and the conversion left the vertical on Y where every other optical file in the corpus is Z-up. Nothing in the file says so; only the data dictionary does.

BOARD_MM module-attribute

BOARD_MM = (433.0, 238.0)

Wii Balance Board sensing area, width by depth, in millimetres.

Multiplying the normalised centre of pressure by these gives millimetres, which is what makes balance data comparable with the optical collections.

read_qualisys

read_qualisys(path: str, drop_gaps: bool = True) -> MotionRecord

Qualisys or Qualisys-style TSV, in all three header shapes found in the corpus.

The shapes differ in what follows the ten KEY<TAB>value metadata lines:

• a <marker> X column-name row, then data (2012, 2015, 2017, 2018, 2019, HpSp);
• nothing, data begins immediately (2022, Bishop 2020);
• a Frame/Time column-name row, then data (MocapNoiseFloor, Solberg 2016).

The shape is detected by trying to parse the eleventh line as numbers, so a file that was exported differently still reads correctly.

Source code in src/micromotion/io.py

def read_qualisys(path: str, drop_gaps: bool = True) -> MotionRecord:
    """Qualisys or Qualisys-style TSV, in all three header shapes found in the corpus.

    The shapes differ in what follows the ten ``KEY<TAB>value`` metadata lines:

    * a ``<marker> X`` column-name row, then data (2012, 2015, 2017, 2018, 2019, HpSp);
    * nothing, data begins immediately (2022, Bishop 2020);
    * a ``Frame``/``Time`` column-name row, then data (MocapNoiseFloor, Solberg 2016).

    The shape is detected by trying to parse the eleventh line as numbers, so a file that
    was exported differently still reads correctly.
    """
    lines = _decode(path)
    meta, n_header = {}, 0
    for i, ln in enumerate(lines[:12]):
        parts = ln.rstrip("\n").split("\t")
        key = parts[0].strip()
        if key in QUALISYS_KEYS:
            meta[key] = [p for p in parts[1:] if p.strip() != ""]
            n_header = i + 1
        else:
            break

    if "MARKER_NAMES" not in meta:
        raise ValueError(f"{path} has no MARKER_NAMES line; not a Qualisys export")
    markers = [m.strip() for m in meta["MARKER_NAMES"]]
    fs_header = float(meta["FREQUENCY"][0])

    def is_numeric(line: str) -> bool:
        try:
            float(line.split("\t")[0])
            return True
        except (ValueError, IndexError):
            return False

    has_frame_time = False
    if is_numeric(lines[n_header]):
        skip = n_header                       # shape 2: no column-name row
    else:
        skip = n_header + 1                   # shape 1 or 3
        first = lines[n_header].split("\t")[0].strip()
        has_frame_time = first.lower() == "frame"

    df = pd.read_csv(path, sep="\t", skiprows=skip, header=None, encoding="latin-1")
    if df.iloc[:, -1].isna().all():
        df = df.iloc[:, :-1]                  # the trailing tab's phantom column
    arr = df.to_numpy(float)

    t = None
    if has_frame_time:
        t = arr[:, 1].copy()
        arr = arr[:, 2:]

    n_expected = 3 * len(markers)
    if arr.shape[1] != n_expected:
        # Trust the columns, not the count. One HpSp file declares 22 markers and carries 7.
        n_actual = arr.shape[1] // 3
        markers = markers[:n_actual] if n_actual <= len(markers) else [
            *markers, *[f"unnamed{i}" for i in range(len(markers), n_actual)]
        ]
        arr = arr[:, : 3 * len(markers)]

    if drop_gaps:
        # A zero triplet is the Qualisys gap code, not a marker at the origin. In HpSp that
        # is 878 948 lines; read as data it would report a marker teleporting to (0,0,0).
        block = arr.reshape(len(arr), -1, 3)
        gap = np.all(block == 0.0, axis=2)
        block[gap] = np.nan
        arr = block.reshape(len(arr), -1)

    channels = [f"{m} {ax}" for m in markers for ax in "XYZ"]
    vertical = "Y" if any(c in path for c in Y_UP_COLLECTIONS) else "Z"
    fs = measured_rate(t) if t is not None and len(t) > 1 else fs_header

    return MotionRecord(
        data=arr,
        fs=fs,
        channels=channels,
        kind="position",
        unit="mm",
        vertical=vertical,
        t=t,
        source=path,
        meta={"header": meta, "nominal_fs": fs_header, "n_markers": len(markers)},
    )

read_sverm

read_sverm(path: str) -> MotionRecord

Sverm plain-header export: Time then s<n>_head_{x,y,z}.

A separate reader because it shares nothing with the Qualisys export but the units. It silently returned zero series to a corpus analysis until one was written for it.

Source code in src/micromotion/io.py

def read_sverm(path: str) -> MotionRecord:
    """Sverm plain-header export: ``Time`` then ``s<n>_head_{x,y,z}``.

    A separate reader because it shares nothing with the Qualisys export but the units. It
    silently returned zero series to a corpus analysis until one was written for it.
    """
    df = pd.read_csv(path, sep="\t")
    if "Time" not in df.columns:
        raise ValueError(f"{path} has no Time column; not a Sverm export")
    t = df["Time"].to_numpy(float)
    cols = [c for c in df.columns if c != "Time"]
    return MotionRecord(
        data=df[cols].to_numpy(float),
        fs=measured_rate(t),
        channels=[c.replace("_", " ") for c in cols],
        kind="position",
        unit="mm",
        t=t,
        source=path,
        meta={"subjects": sorted({c.split("_")[0] for c in cols})},
    )

read_ax3

read_ax3(path: str) -> MotionRecord

Axivity AX3 export with a ts,x,y,z header: Standstill2024 and Taqasim.

The rate is measured from the timestamps and differs meaningfully between files: it is a property of the physical logger, spanning 191.3-207.7 Hz across the 2024 units, and it differs between the two Taqasim sites in the same direction as the effect that record's headline claim reports.

Source code in src/micromotion/io.py

def read_ax3(path: str) -> MotionRecord:
    """Axivity AX3 export with a ``ts,x,y,z`` header: Standstill2024 and Taqasim.

    The rate is measured from the timestamps and differs meaningfully between files: it is
    a property of the physical logger, spanning 191.3-207.7 Hz across the 2024 units, and
    it differs between the two Taqasim sites in the same direction as the effect that
    record's headline claim reports.
    """
    df = pd.read_csv(path, sep="\t")
    ts = pd.to_datetime(df["ts"])
    t = (ts - ts.iloc[0]).dt.total_seconds().to_numpy()
    return MotionRecord(
        data=df[["x", "y", "z"]].to_numpy(float),
        fs=measured_rate(t),
        channels=["x", "y", "z"],
        kind="acceleration",
        unit="g",
        t=t,
        t0=ts.iloc[0],
        source=path,
        meta={"device_slot": (re.match(r"[A-Z]+(\d+)", os.path.basename(path)[:-4]) or [None, None])[1]},
    )

read_phone

read_phone(path: str, trim_clap_s: float = 0.0) -> MotionRecord

Physics Toolbox phone log, cleaned tab-separated form.

ax/ay/az are linear acceleration in m/s^2. They are not in g, whatever an older version of the data dictionary said; reading them as g inflated every quantity of motion in this project by 9.80665 until 2026-07-28. Only gF* is in g.

trim_clap_s drops that many seconds from each end. The StillStanding365 recordings open and close with a synchronisation clap which is a timing marker, not movement, and the deposited pipeline trims 12 s.

Source code in src/micromotion/io.py

def read_phone(path: str, trim_clap_s: float = 0.0) -> MotionRecord:
    """Physics Toolbox phone log, cleaned tab-separated form.

    ``ax``/``ay``/``az`` are linear acceleration in m/s^2. They are not in g, whatever an
    older version of the data dictionary said; reading them as g inflated every quantity of
    motion in this project by 9.80665 until 2026-07-28. Only ``gF*`` is in g.

    ``trim_clap_s`` drops that many seconds from each end. The StillStanding365 recordings
    open and close with a synchronisation clap which is a timing marker, not movement, and
    the deposited pipeline trims 12 s.
    """
    df = pd.read_csv(path, sep="\t")
    t = df["time"].to_numpy(float)
    m = np.ones(len(t), bool)
    if trim_clap_s:
        m = (t >= trim_clap_s) & (t <= t[-1] - trim_clap_s)
    cols = [c for c in ("ax", "ay", "az") if c in df.columns]
    data = df[cols].to_numpy(float)[m]

    # Missing samples are written as exact zeros, not NaN. Every file opens with a few.
    zero_rows = np.all(data == 0.0, axis=1)
    data[zero_rows] = np.nan

    return MotionRecord(
        data=data,
        fs=measured_rate(t),
        channels=cols,
        kind="acceleration",
        unit="m/s^2",
        t=t[m],
        source=path,
        meta={"trimmed_s": trim_clap_s, "n_zero_rows": int(zero_rows.sum()),
              "extra": {c: df[c].to_numpy(float)[m] for c in df.columns
                        if c in ("wx", "wy", "wz", "gFx", "gFy", "gFz")}},
    )

read_equivital

read_equivital(path: str) -> MotionRecord

Equivital physiology CSV: Stillness2025 accelerometer, ECG, respiration or RR.

Four of the five files per participant are delimited by comma-and-space and the fifth by a bare comma, so the separator is handled rather than assumed. The accelerometer is in raw counts, calibrated to g by the median vector magnitude, since a person standing still averages one g.

The rate is measured from the timestamps by span over count. Taking the median interval instead returns exactly 250 Hz for a recording that runs at 256, because the timestamps are rounded to whole milliseconds; that error was live in the deposited quantity of motion until 2026-07-28.

Source code in src/micromotion/io.py

def read_equivital(path: str) -> MotionRecord:
    """Equivital physiology CSV: Stillness2025 accelerometer, ECG, respiration or RR.

    Four of the five files per participant are delimited by comma-and-space and the fifth by
    a bare comma, so the separator is handled rather than assumed. The accelerometer is in
    raw counts, calibrated to g by the median vector magnitude, since a person standing
    still averages one g.

    The rate is measured from the timestamps by span over count. Taking the median interval
    instead returns exactly 250 Hz for a recording that runs at 256, because the timestamps
    are rounded to whole milliseconds; that error was live in the deposited quantity of
    motion until 2026-07-28.
    """
    df = pd.read_csv(path, skipinitialspace=True)
    tcol = df.columns[0]
    ts = pd.to_datetime(df[tcol], format="mixed", utc=True)
    t = (ts - ts.iloc[0]).dt.total_seconds().to_numpy()
    cols = list(df.columns[1:])

    # copy(): pandas returns a read-only view under copy-on-write, and the rail masking
    # below writes in place. Without this, reading a respiration file raises.
    data = df[cols].to_numpy(float).copy()
    unit, kind = "counts", "acceleration"
    if any("Accelerometer" in c for c in cols):
        mag = np.median(np.linalg.norm(data, axis=1))
        data = data / mag                      # counts -> g
        unit = "g"
    elif "Breathing" in cols:
        kind, unit = "respiration", "adc"
        # 0 and 1023 are the rails of a 10-bit converter, not waveform.
        data[(data <= 0) | (data >= 1023)] = np.nan

    return MotionRecord(
        data=data,
        fs=measured_rate(t),
        channels=cols,
        kind=kind,
        unit=unit,
        t=t,
        t0=ts.iloc[0],
        source=path,
    )

read_balance_board

read_balance_board(path: str) -> MotionRecord

Wii balance board dump: headerless, space-delimited, irregularly sampled.

Eight columns: timestamp in milliseconds, four corner load cells, total load, and centre of pressure in x and y, all normalised to 0-1.

Samples with no load carry a centre of pressure of exactly (0.5, 0.5), the midpoint of the board, which is a plausible-looking value for nobody standing on it. Those are returned as NaN. They are 13.8 per cent of the HpSp balance data.

Source code in src/micromotion/io.py

def read_balance_board(path: str) -> MotionRecord:
    """Wii balance board dump: headerless, space-delimited, irregularly sampled.

    Eight columns: timestamp in milliseconds, four corner load cells, total load, and
    centre of pressure in x and y, all normalised to 0-1.

    Samples with no load carry a centre of pressure of exactly (0.5, 0.5), the midpoint of
    the board, which is a plausible-looking value for nobody standing on it. Those are
    returned as NaN. They are 13.8 per cent of the HpSp balance data.
    """
    arr = np.loadtxt(path)
    if arr.ndim != 2 or arr.shape[1] < 8:
        raise ValueError(f"{path} is not an 8-column balance dump")
    t = arr[:, 0] / 1000.0
    data = arr[:, 6:8].copy()
    load = arr[:, 5]
    data[(load <= 0) | ((data[:, 0] == 0.5) & (data[:, 1] == 0.5))] = np.nan
    return MotionRecord(
        data=data,
        fs=measured_rate(t) if len(t) > 1 else float("nan"),
        channels=["cop_x", "cop_y"],
        kind="position",
        unit="normalised",
        t=t,
        source=path,
        meta={"load": load, "cells": arr[:, 1:5], "irregular": True},
    )

read_fnirs

read_fnirs(path: str) -> MotionRecord

Artinis Brite wide export: 96 haemoglobin channels plus two IMUs.

Only the accelerometer and gyroscope are returned as data; the haemoglobin channels are kept in meta because they are not motion. The accelerometer is in milli-g.

Column 0 is a sample index, not time, so the rate comes from the header. The absolute start can be decoded from the TimeStampHi/Lo pair if it is needed.

Source code in src/micromotion/io.py

def read_fnirs(path: str) -> MotionRecord:
    """Artinis Brite wide export: 96 haemoglobin channels plus two IMUs.

    Only the accelerometer and gyroscope are returned as ``data``; the haemoglobin channels
    are kept in ``meta`` because they are not motion. The accelerometer is in milli-g.

    Column 0 is a sample index, not time, so the rate comes from the header. The absolute
    start can be decoded from the ``TimeStampHi``/``Lo`` pair if it is needed.
    """
    df = pd.read_csv(path, sep="\t")
    acc = [c for c in df.columns if c.startswith("ACC_")]
    gyr = [c for c in df.columns if c.startswith("GYR_")]
    if not acc:
        raise ValueError(f"{path} has no ACC_ columns; not an fNIRS export")
    dev = acc[0].split("_")[-2]
    acc = [c for c in acc if dev in c][:3]
    gyr = [c for c in gyr if dev in c][:3]

    legend = os.path.join(os.path.dirname(path), "channel_legend.txt")
    fs = 75.0
    if os.path.exists(legend):
        m = re.search(r"sample rate:\s*([\d.]+)", _io.open(legend, encoding="latin-1").read())
        if m:
            fs = float(m.group(1))

    return MotionRecord(
        data=df[acc].to_numpy(float) / 1000.0,     # milli-g -> g
        fs=fs,
        channels=acc,
        kind="acceleration",
        unit="g",
        source=path,
        meta={
            "device": dev,
            "gyro_deg_s": df[gyr].to_numpy(float) if gyr else None,
            "haemoglobin": [c for c in df.columns if c.endswith(("O2Hb", "HHb"))],
        },
    )

sniff

sniff(path: str) -> str

Identify a file's layout from its first lines.

Source code in src/micromotion/io.py

def sniff(path: str) -> str:
    """Identify a file's layout from its first lines."""
    head = "".join(_decode(path)[:3])
    first = head.split("\n")[0]
    if first.startswith("NO_OF_FRAMES"):
        return "qualisys"
    if first.startswith("Time\t") and "_head_" in first:
        return "sverm"
    if first.startswith("ts\tx\ty\tz"):
        return "ax3"
    if first.startswith("time\t") and "\tax\t" in first:
        return "phone"
    if first.startswith("sample_index\t"):
        return "fnirs"
    if first.startswith("DateTime"):
        return "equivital"
    if re.fullmatch(r"[\d.\s eE+-]+", first) and len(first.split()) == 8:
        return "balance"
    raise ValueError(f"cannot identify the layout of {path}")

read

read(path: str, **kw) -> MotionRecord

Read any corpus motion file, dispatching on content.

Source code in src/micromotion/io.py

def read(path: str, **kw) -> MotionRecord:
    """Read any corpus motion file, dispatching on content."""
    return _READERS[sniff(path)](path, **kw)

micromotion.record

The common structure every reader returns.

The corpus holds sixteen distinct file layouts across fourteen collections, and each one has already cost an analysis at least once: an axis convention that differs in a single edition, a rate that three documents disagree about, a gap code that looks like a valid measurement. A reader's job is to absorb all of that and hand back the same object regardless.

MotionRecord dataclass

One recording, in a form the rest of the package can use without special cases.

Attributes:

Name	Type	Description
data	ndarray	(n_samples, n_channels) float array. Gaps are NaN, never a sentinel.
fs	float	Sampling rate in Hz, measured from the file where a timebase exists and taken from the header only where it does not.
channels	list[str]	Column names, one per column of data.
kind	str	"position" or "acceleration": what the numbers are, which decides whether quantity of motion differentiates or integrates.
unit	str	"mm", "m", "g", "m/s^2" or "counts".
vertical	str	Which axis letter is up. Not always Z: the 2019 championship is Y-up, alone in the corpus.
t	ndarray | None	Timestamps in seconds where the file carries them, otherwise None.
t0	object | None	Absolute start time where the file carries one, otherwise None.
source	str	Path the record was read from.
meta	dict	Anything else the header held, unmodified.

Source code in src/micromotion/record.py

@dataclass
class MotionRecord:
    """One recording, in a form the rest of the package can use without special cases.

    Attributes
    ----------
    data
        (n_samples, n_channels) float array. Gaps are NaN, never a sentinel.
    fs
        Sampling rate in Hz, measured from the file where a timebase exists and taken from
        the header only where it does not.
    channels
        Column names, one per column of ``data``.
    kind
        ``"position"`` or ``"acceleration"``: what the numbers are, which decides whether
        quantity of motion differentiates or integrates.
    unit
        ``"mm"``, ``"m"``, ``"g"``, ``"m/s^2"`` or ``"counts"``.
    vertical
        Which axis letter is up. Not always Z: the 2019 championship is Y-up, alone in the
        corpus.
    t
        Timestamps in seconds where the file carries them, otherwise None.
    t0
        Absolute start time where the file carries one, otherwise None.
    source
        Path the record was read from.
    meta
        Anything else the header held, unmodified.
    """

    data: np.ndarray = field(repr=False)
    fs: float
    channels: list[str]
    kind: str = "position"
    unit: str = "mm"
    vertical: str = "Z"
    t: np.ndarray | None = field(default=None, repr=False)
    t0: object | None = None
    source: str = ""
    meta: dict = field(default_factory=dict)

    def __post_init__(self):
        self.data = np.asarray(self.data, float)
        if self.data.ndim == 1:
            self.data = self.data[:, None]
        if len(self.channels) != self.data.shape[1]:
            raise ValueError(
                f"{len(self.channels)} channel names for {self.data.shape[1]} columns "
                f"in {self.source}"
            )

    @property
    def n_samples(self) -> int:
        return self.data.shape[0]

    @property
    def duration_s(self) -> float:
        return self.n_samples / self.fs

    @property
    def markers(self) -> list[str]:
        """Marker names, for records whose channels are ``<marker> X/Y/Z`` triplets."""
        seen, out = set(), []
        for c in self.channels:
            name = c.rsplit(" ", 1)[0]
            if name not in seen:
                seen.add(name)
                out.append(name)
        return out

    def marker(self, name: str) -> np.ndarray:
        """The (n_samples, 3) block for one marker, by name.

        Reading marker columns by name rather than by position is not a nicety. Six HpSp
        files break the documented 22-marker order, and the README's positional quick-start
        mis-assigned markers in every one of them until it was rewritten.
        """
        idx = [i for i, c in enumerate(self.channels) if c.rsplit(" ", 1)[0] == name]
        if not idx:
            raise KeyError(f"no marker {name!r} in {self.source}; have {self.markers}")
        return self.data[:, idx]

    def missing_fraction(self) -> float:
        """Proportion of the array that is NaN."""
        return float(np.isnan(self.data).mean())

    def qom(self, **kw):
        """Quantity of motion for this record, with kind and unit already filled in.

        Passing a whole record computes across every channel at once, which is rarely what
        you want for a multi-marker file; select a marker first.
        """
        from .qom import qom as _qom

        kw.setdefault("kind", self.kind)
        kw.setdefault("unit", self.unit)
        return _qom(self.data, self.fs, **kw)

markers property

markers: list[str]

Marker names, for records whose channels are <marker> X/Y/Z triplets.

marker

marker(name: str) -> np.ndarray

The (n_samples, 3) block for one marker, by name.

Reading marker columns by name rather than by position is not a nicety. Six HpSp files break the documented 22-marker order, and the README's positional quick-start mis-assigned markers in every one of them until it was rewritten.

Source code in src/micromotion/record.py

def marker(self, name: str) -> np.ndarray:
    """The (n_samples, 3) block for one marker, by name.

    Reading marker columns by name rather than by position is not a nicety. Six HpSp
    files break the documented 22-marker order, and the README's positional quick-start
    mis-assigned markers in every one of them until it was rewritten.
    """
    idx = [i for i, c in enumerate(self.channels) if c.rsplit(" ", 1)[0] == name]
    if not idx:
        raise KeyError(f"no marker {name!r} in {self.source}; have {self.markers}")
    return self.data[:, idx]

missing_fraction

missing_fraction() -> float

Proportion of the array that is NaN.

Source code in src/micromotion/record.py

def missing_fraction(self) -> float:
    """Proportion of the array that is NaN."""
    return float(np.isnan(self.data).mean())

qom

qom(**kw)

Quantity of motion for this record, with kind and unit already filled in.

Passing a whole record computes across every channel at once, which is rarely what you want for a multi-marker file; select a marker first.

Source code in src/micromotion/record.py

def qom(self, **kw):
    """Quantity of motion for this record, with kind and unit already filled in.

    Passing a whole record computes across every channel at once, which is rarely what
    you want for a multi-marker file; select a marker first.
    """
    from .qom import qom as _qom

    kw.setdefault("kind", self.kind)
    kw.setdefault("unit", self.unit)
    return _qom(self.data, self.fs, **kw)

Spectral

micromotion.spectral

Spectral helpers: finding the physiological peaks in a motion signal.

A body-worn accelerometer on someone standing still picks up two rhythms that are not movement in the intentional sense. Respiration sits at roughly 0.2-0.4 Hz and the ballistocardiac impulse, the recoil of the heart ejecting blood, at roughly 0.8-1.8 Hz. Both are inside or adjacent to the micromotion band, so isolating postural motion means locating them first.

CARDIAC_BAND module-attribute

CARDIAC_BAND = (0.7, 2.2)

Hz. 42-132 bpm, which covers rest through mild exertion.

RESPIRATORY_BAND module-attribute

RESPIRATORY_BAND = (0.1, 0.5)

Hz. 6-30 breaths per minute.

cardiac_peak

cardiac_peak(x, fs: float, window_s: float = 60.0) -> float

Dominant frequency in the cardiac band, in Hz.

Pass the acceleration magnitude. Multiply by 60 for beats per minute. Returns NaN if the recording is too short for the band to be resolved.

Source code in src/micromotion/spectral.py

def cardiac_peak(x, fs: float, window_s: float = 60.0) -> float:
    """Dominant frequency in the cardiac band, in Hz.

    Pass the acceleration magnitude. Multiply by 60 for beats per minute. Returns NaN if
    the recording is too short for the band to be resolved.
    """
    return _peak(x, fs, CARDIAC_BAND, window_s)

respiratory_peak

respiratory_peak(x, fs: float, window_s: float = 120.0) -> float

Dominant frequency in the respiratory band, in Hz.

Needs a long window: at 0.2 Hz, a 60-second segment holds twelve cycles, which is few enough that the estimate wanders.

Source code in src/micromotion/spectral.py

def respiratory_peak(x, fs: float, window_s: float = 120.0) -> float:
    """Dominant frequency in the respiratory band, in Hz.

    Needs a long window: at 0.2 Hz, a 60-second segment holds twelve cycles, which is few
    enough that the estimate wanders.
    """
    return _peak(x, fs, RESPIRATORY_BAND, window_s)

band_power

band_power(x, fs: float, band: tuple[float, float], window_s: float = 60.0) -> float

Integrated power between two frequencies.

Source code in src/micromotion/spectral.py

def band_power(x, fs: float, band: tuple[float, float], window_s: float = 60.0) -> float:
    """Integrated power between two frequencies."""
    x = np.asarray(x, float)
    nper = int(min(len(x), fs * window_s))
    f, p = signal.welch(signal.detrend(x), fs, nperseg=nper)
    m = (f >= band[0]) & (f <= band[1])
    return float(np.trapezoid(p[m], f[m])) if m.sum() > 1 else float("nan")

spectral_peak

spectral_peak(x, fs: float, band: tuple[float, float], window_s: float = 60.0) -> dict

Peak frequency in a band, with a signal-to-noise ratio against the band's own median.

The ratio is what makes the peak reportable. Every spectrum has a maximum somewhere inside any band you choose; whether it rises above the surrounding noise decides whether a rhythm is present at all. Below about 3, treat the peak as absent rather than weak.

Source code in src/micromotion/spectral.py

def spectral_peak(x, fs: float, band: tuple[float, float],
                  window_s: float = 60.0) -> dict:
    """Peak frequency in a band, with a signal-to-noise ratio against the band's own median.

    The ratio is what makes the peak reportable. Every spectrum has a maximum somewhere
    inside any band you choose; whether it rises above the surrounding noise decides whether
    a rhythm is present at all. Below about 3, treat the peak as absent rather than weak.
    """
    x = np.asarray(x, float)
    if len(x) < fs * 10:
        return {"freq": float("nan"), "power": float("nan"), "snr": float("nan")}
    nper = int(min(len(x), fs * window_s))
    f, p = signal.welch(signal.detrend(x), fs, nperseg=nper)
    m = (f >= band[0]) & (f <= band[1])
    if not m.any():
        return {"freq": float("nan"), "power": float("nan"), "snr": float("nan")}
    k = int(np.argmax(p[m]))
    med = float(np.median(p[m]))
    return {"freq": float(f[m][k]), "power": float(p[m][k]),
            "snr": float(p[m][k] / med) if med > 0 else float("nan")}

band_rms

band_rms(x, fs: float, band: tuple[float, float]) -> float

Root-mean-square amplitude within a band, in the input units.

Rate-independent, unlike an integrated measure, which makes it the right thing to quote when comparing instruments that sample at different rates.

Source code in src/micromotion/spectral.py

def band_rms(x, fs: float, band: tuple[float, float]) -> float:
    """Root-mean-square amplitude within a band, in the input units.

    Rate-independent, unlike an integrated measure, which makes it the right thing to quote
    when comparing instruments that sample at different rates.
    """
    from .filters import bandpass

    return float(np.sqrt(np.mean(bandpass(np.asarray(x, float), fs, *band) ** 2)))

band_power_fraction

band_power_fraction(x, fs: float, bands: dict, window_s: float = 60.0) -> dict

Proportion of total power falling in each named band.

Pass something like {"respiratory": (0.1, 0.5), "cardiac": (0.7, 2.2)}. This is how the finding that roughly 38 per cent of the chest-phone signal power sits at the heart rate was established, which is in turn why the compensated quantity-of-motion variant exists.

Source code in src/micromotion/spectral.py

def band_power_fraction(x, fs: float, bands: dict, window_s: float = 60.0) -> dict:
    """Proportion of total power falling in each named band.

    Pass something like ``{"respiratory": (0.1, 0.5), "cardiac": (0.7, 2.2)}``. This is how
    the finding that roughly 38 per cent of the chest-phone signal power sits at the heart
    rate was established, which is in turn why the compensated quantity-of-motion variant
    exists.
    """
    x = np.asarray(x, float)
    nper = int(min(len(x), fs * window_s))
    f, p = signal.welch(signal.detrend(x), fs, nperseg=nper)
    total = float(np.trapezoid(p, f))
    out = {}
    for name, (lo, hi) in bands.items():
        m = (f >= lo) & (f <= hi)
        out[name] = (float(np.trapezoid(p[m], f[m]) / total)
                     if m.sum() > 1 and total > 0 else float("nan"))
    return out

mean_frequency

mean_frequency(x, fs: float, band: tuple[float, float] = (0.1, 5.0), window_s: float = 60.0) -> float

Power-weighted mean frequency within a band, in Hz.

A different statistic from the peak: it moves when energy shifts between frequencies even if the dominant one does not change, so it is the more sensitive of the two to a gradual change in how a person is standing.

Source code in src/micromotion/spectral.py

def mean_frequency(x, fs: float, band: tuple[float, float] = (0.1, 5.0),
                   window_s: float = 60.0) -> float:
    """Power-weighted mean frequency within a band, in Hz.

    A different statistic from the peak: it moves when energy shifts between frequencies even
    if the dominant one does not change, so it is the more sensitive of the two to a gradual
    change in how a person is standing.
    """
    x = np.asarray(x, float)
    nper = int(min(len(x), fs * window_s))
    f, p = signal.welch(signal.detrend(x), fs, nperseg=nper)
    m = (f >= band[0]) & (f <= band[1])
    return float(np.sum(f[m] * p[m]) / np.sum(p[m])) if p[m].sum() > 0 else float("nan")

detect_breaths

detect_breaths(x, fs: float, band: tuple[float, float] = RESPIRATORY_BAND, prominence_sd: float = 0.4, min_period_s: float = 2.0) -> dict

Breath timing from a respiration belt or a band-limited motion signal.

Returns peak and trough times and the cycle durations between them. Peaks and troughs are both returned because inhale and exhale are not symmetric, and a rate alone hides that.

Source code in src/micromotion/spectral.py

def detect_breaths(x, fs: float, band: tuple[float, float] = RESPIRATORY_BAND,
                   prominence_sd: float = 0.4, min_period_s: float = 2.0) -> dict:
    """Breath timing from a respiration belt or a band-limited motion signal.

    Returns peak and trough times and the cycle durations between them. Peaks and troughs are
    both returned because inhale and exhale are not symmetric, and a rate alone hides that.
    """
    from .filters import bandpass

    y = bandpass(np.asarray(x, float), fs, *band)
    kw = dict(distance=max(1, int(min_period_s * fs)),
              prominence=prominence_sd * np.std(y))
    peaks, _ = signal.find_peaks(y, **kw)
    troughs, _ = signal.find_peaks(-y, **kw)
    cycles = np.diff(peaks) / fs if len(peaks) > 1 else np.array([])
    return {"peaks_s": peaks / fs, "troughs_s": troughs / fs, "cycle_s": cycles,
            "rate_per_min": float(60.0 / np.median(cycles)) if len(cycles) else float("nan"),
            "n_breaths": len(peaks)}

detect_breaths_adaptive

detect_breaths_adaptive(x, fs: float, band: tuple[float, float] = RESPIRATORY_BAND, vel_frac: float = 0.55, baseline_hz: float = 0.2) -> dict

Breath detection that rejects chest movement which is not breathing.

Peak-and-trough detection, which :func:detect_breaths does, counts any sufficiently prominent bump. On a belt worn by someone standing that includes postural sway, weight shifts and swallows, all of which look like small breaths.

This asks a different question. A breath is a sustained rise in chest expansion whose velocity exceeds a threshold set from the signal's own positive-derivative mean, and which crosses an adaptive baseline -- a heavily low-passed copy of the signal rather than zero. A rise that never crosses that baseline did not start from an exhaled state and is discarded. That rejection step is what removes the sway.

Returns inspiration and expiration onsets, cycle durations and the rate. Expiration onset is taken as the end of the rise, which assumes passive expiration.

After Finn Upham's respiration work, reimplemented and used with permission.

Do not reach for this by default, on the evidence available here. It was added on the expectation that rejecting non-breath rises would beat plain peak detection, and measured against it that expectation did not hold. On twelve HpSp respiration-belt recordings the two agree: median error against the spectral estimate 1.82 breaths per minute for both. On eight Stillness2025 chest accelerometers, which is the case the rejection step was supposed to help, it is markedly worse -- median error against the same participant's belt 10.3 breaths per minute against 3.6 for :func:detect_breaths, over-counting throughout.

That may be the parameters rather than the idea; the velocity threshold is derived from a belt's amplitude distribution and an accelerometer's is not the same shape. It is kept because the approach is sound in its original setting and because someone should be able to tune it, not because it is currently the better detector.

Source code in src/micromotion/spectral.py

def detect_breaths_adaptive(x, fs: float, band: tuple[float, float] = RESPIRATORY_BAND,
                            vel_frac: float = 0.55, baseline_hz: float = 0.2) -> dict:
    """Breath detection that rejects chest movement which is not breathing.

    Peak-and-trough detection, which :func:`detect_breaths` does, counts any sufficiently
    prominent bump. On a belt worn by someone standing that includes postural sway, weight
    shifts and swallows, all of which look like small breaths.

    This asks a different question. A breath is a sustained rise in chest expansion whose
    velocity exceeds a threshold set from the signal's own positive-derivative mean, and which
    crosses an adaptive baseline -- a heavily low-passed copy of the signal rather than zero.
    A rise that never crosses that baseline did not start from an exhaled state and is
    discarded. That rejection step is what removes the sway.

    Returns inspiration and expiration onsets, cycle durations and the rate. Expiration onset
    is taken as the end of the rise, which assumes passive expiration.

    After Finn Upham's respiration work, reimplemented and used with permission.

    Do not reach for this by default, on the evidence available here. It was added on the
    expectation that rejecting non-breath rises would beat plain peak detection, and measured
    against it that expectation did not hold. On twelve HpSp respiration-belt recordings the
    two agree: median error against the spectral estimate 1.82 breaths per minute for both.
    On eight Stillness2025 chest accelerometers, which is the case the rejection step was
    supposed to help, it is markedly worse -- median error against the same participant's belt
    10.3 breaths per minute against 3.6 for :func:`detect_breaths`, over-counting throughout.

    That may be the parameters rather than the idea; the velocity threshold is derived from a
    belt's amplitude distribution and an accelerometer's is not the same shape. It is kept
    because the approach is sound in its original setting and because someone should be able
    to tune it, not because it is currently the better detector.
    """
    from .filters import bandpass, lowpass

    x = np.asarray(x, float)
    y = bandpass(x, fs, *band)
    d = np.gradient(y, 1.0 / fs)

    pos = d[d > 0]
    if not len(pos):
        return {"inspiration_s": np.array([]), "expiration_s": np.array([]),
                "cycle_s": np.array([]), "rate_per_min": float("nan"), "n_breaths": 0}
    thresh = pos.mean() * vel_frac

    # Crossings of an adaptive baseline, not of zero: the belt drifts.
    base = lowpass(y, fs, min(baseline_hz, fs / 2 * 0.4), order=2)
    above = y > base

    rising = d > thresh
    edges = np.diff(rising.astype(int))
    starts = np.flatnonzero(edges == 1) + 1
    ends = np.flatnonzero(edges == -1) + 1
    if len(ends) and len(starts) and ends[0] < starts[0]:
        ends = ends[1:]
    n = min(len(starts), len(ends))
    starts, ends = starts[:n], ends[:n]

    insp, expi = [], []
    for s, e in zip(starts, ends):
        # A genuine breath rises through the baseline; a sway bump does not.
        if e > s and above[s:e].any() and not above[s:e].all():
            insp.append(s / fs)
            expi.append(e / fs)
    insp, expi = np.asarray(insp), np.asarray(expi)
    cycles = np.diff(insp) if len(insp) > 1 else np.array([])
    return {
        "inspiration_s": insp, "expiration_s": expi, "cycle_s": cycles,
        "rate_per_min": float(60.0 / np.median(cycles)) if len(cycles) else float("nan"),
        "n_breaths": len(insp),
    }

Dynamics

micromotion.dynamics

Scaling and nonlinearity measures for postural time series.

These describe how movement is structured in time rather than how much of it there is. They are the part of this package with reuse value beyond the standstill corpus, and they are also the part where a subtly wrong implementation produces a plausible number rather than an error. Every function here is covered by a test against a process whose answer is known in advance; see tests/test_dynamics.py.

One warning carried over from the corpus. These methods read across scales, so they are the ones that upsampling corrupts. Resample downwards only, and use :func:micromotion.resample.to_rate, which refuses to do otherwise.

iaaft

iaaft(x, iters: int = 100, rng=None) -> np.ndarray

Iterative amplitude-adjusted Fourier transform surrogate.

Preserves both the amplitude distribution and the power spectrum while destroying any nonlinear structure, so it is the null hypothesis "this series is a linear Gaussian process observed through a monotonic transform".

Source code in src/micromotion/dynamics.py

def iaaft(x, iters: int = 100, rng=None) -> np.ndarray:
    """Iterative amplitude-adjusted Fourier transform surrogate.

    Preserves both the amplitude distribution and the power spectrum while destroying any
    nonlinear structure, so it is the null hypothesis "this series is a linear Gaussian
    process observed through a monotonic transform".
    """
    rng = rng or np.random.default_rng()
    x = np.asarray(x, float)
    n = len(x)
    amp = np.abs(np.fft.rfft(x))
    srt = np.sort(x)
    y = rng.permutation(x)
    for _ in range(iters):
        Y = np.fft.rfft(y)
        y = np.fft.irfft(amp * np.exp(1j * np.angle(Y)), n=n)
        y = srt[np.argsort(np.argsort(y))]
    return y

phase_surrogate

phase_surrogate(x, rng=None) -> np.ndarray

Phase-randomised surrogate: preserves the spectrum, not the distribution.

Cheaper than :func:iaaft and a weaker null, since a non-Gaussian amplitude distribution alone can make a linear series look nonlinear against it.

Source code in src/micromotion/dynamics.py

def phase_surrogate(x, rng=None) -> np.ndarray:
    """Phase-randomised surrogate: preserves the spectrum, not the distribution.

    Cheaper than :func:`iaaft` and a weaker null, since a non-Gaussian amplitude
    distribution alone can make a linear series look nonlinear against it.
    """
    rng = rng or np.random.default_rng()
    x = np.asarray(x, float)
    X = np.fft.rfft(x)
    ph = rng.uniform(0, 2 * np.pi, len(X))
    ph[0] = 0
    if len(x) % 2 == 0:
        ph[-1] = 0
    return np.fft.irfft(np.abs(X) * np.exp(1j * ph), n=len(x))

circular_shift_surrogate

circular_shift_surrogate(x, rng=None, margin: float = 0.1) -> np.ndarray

Rotate a series in time.

Preserves everything about the series itself and destroys only its alignment with another, so it is the right null for a correlation between two recordings and the wrong one for a property of a single series.

Source code in src/micromotion/dynamics.py

def circular_shift_surrogate(x, rng=None, margin: float = 0.1) -> np.ndarray:
    """Rotate a series in time.

    Preserves everything about the series itself and destroys only its alignment with
    another, so it is the right null for a correlation between two recordings and the wrong
    one for a property of a single series.
    """
    rng = rng or np.random.default_rng()
    x = np.asarray(x)
    n = len(x)
    k = rng.integers(int(margin * n), int((1 - margin) * n))
    return np.roll(x, k)

surrogate_test

surrogate_test(x, statistic, n: int = 99, method=iaaft, rng=None) -> dict

Compare a statistic against its surrogate distribution.

Returns the observed value, the surrogate mean and standard deviation, a z score and a two-sided p value. The p value uses the standard (count + 1) / (n + 1) form, which cannot return zero.

Source code in src/micromotion/dynamics.py

def surrogate_test(x, statistic, n: int = 99, method=iaaft, rng=None) -> dict:
    """Compare a statistic against its surrogate distribution.

    Returns the observed value, the surrogate mean and standard deviation, a z score and a
    two-sided p value. The p value uses the standard ``(count + 1) / (n + 1)`` form, which
    cannot return zero.
    """
    rng = rng or np.random.default_rng()
    obs = float(statistic(np.asarray(x, float)))
    null = np.array([statistic(method(x, rng=rng)) for _ in range(n)])
    sd = null.std()
    z = (obs - null.mean()) / (sd + 1e-30)
    p = (np.sum(np.abs(null - null.mean()) >= abs(obs - null.mean())) + 1) / (n + 1)
    return {"observed": obs, "null_mean": float(null.mean()), "null_sd": float(sd),
            "z": float(z), "p": float(p), "n_surrogates": n}

trev

trev(x, tau: int = 1) -> float

Time-reversal asymmetry.

A linear Gaussian process looks the same run backwards; most nonlinear ones do not. The statistic is the skew of the lag-tau increments, normalised by their variance.

Source code in src/micromotion/dynamics.py

def trev(x, tau: int = 1) -> float:
    """Time-reversal asymmetry.

    A linear Gaussian process looks the same run backwards; most nonlinear ones do not. The
    statistic is the skew of the lag-``tau`` increments, normalised by their variance.
    """
    x = np.asarray(x, float)
    d = x[tau:] - x[:-tau]
    return float(np.mean(d**3) / (np.mean(d**2) ** 1.5 + 1e-30))

dfa

dfa(x, smin: int | None = None, smax: int | None = None, nsc: int = 18, order: int = 1) -> dict

Detrended fluctuation analysis.

Returns the scaling exponent alpha. For reference, 0.5 is white noise, 1.0 is pink noise, and 1.5 is Brownian motion. Values above 0.5 mean the series persists: a deviation tends to be followed by more of the same.

Source code in src/micromotion/dynamics.py

def dfa(x, smin: int | None = None, smax: int | None = None, nsc: int = 18,
        order: int = 1) -> dict:
    """Detrended fluctuation analysis.

    Returns the scaling exponent ``alpha``. For reference, 0.5 is white noise, 1.0 is pink
    noise, and 1.5 is Brownian motion. Values above 0.5 mean the series persists: a
    deviation tends to be followed by more of the same.
    """
    x = np.asarray(x, float)
    x = x[np.isfinite(x)]
    n = len(x)
    smin = smin or 8
    smax = smax or n // 4
    if n < 100 or smax <= smin:
        return {"alpha": float("nan"), "scales": np.array([]), "F": np.array([])}
    Y = np.cumsum(x - x.mean())
    scales = np.unique(np.round(np.logspace(np.log10(smin), np.log10(smax), nsc)).astype(int))
    F = []
    for s in scales:
        nseg = n // s
        if nseg < 4:
            F.append(np.nan)
            continue
        seg = Y[: nseg * s].reshape(nseg, s)
        t = np.arange(s)
        V = np.polynomial.polynomial.polyvander(t, order)
        coef, *_ = np.linalg.lstsq(V, seg.T, rcond=None)
        resid = seg.T - V @ coef
        F.append(np.sqrt(np.mean(resid**2)))
    F = np.array(F, float)
    m = np.isfinite(F) & (F > 0)
    if m.sum() < 5:
        return {"alpha": float("nan"), "scales": scales, "F": F}
    alpha = float(np.polyfit(np.log(scales[m]), np.log(F[m]), 1)[0])
    return {"alpha": alpha, "scales": scales, "F": F}

mfdfa

mfdfa(x, qs=None, smin: int = 16, smax: int | None = None, nsc: int = 20, order: int = 1) -> dict | None

Multifractal detrended fluctuation analysis.

Returns h(q), the singularity spectrum, the generalised Hurst exponent h2 and the spectrum width. A width near zero means one scaling exponent describes the whole series; a wide spectrum means different parts of it scale differently.

A width above about 2 for postural data is a signal to check the preprocessing rather than to celebrate: widths up to 6.6 were once produced entirely by upsampling.

Source code in src/micromotion/dynamics.py

def mfdfa(x, qs=None, smin: int = 16, smax: int | None = None, nsc: int = 20,
          order: int = 1) -> dict | None:
    """Multifractal detrended fluctuation analysis.

    Returns ``h(q)``, the singularity spectrum, the generalised Hurst exponent ``h2`` and
    the spectrum ``width``. A width near zero means one scaling exponent describes the whole
    series; a wide spectrum means different parts of it scale differently.

    A width above about 2 for postural data is a signal to check the preprocessing rather
    than to celebrate: widths up to 6.6 were once produced entirely by upsampling.
    """
    qs = np.arange(-5, 5.1, 0.5) if qs is None else np.asarray(qs, float)
    x = np.asarray(x, float)
    x = x[np.isfinite(x)]
    N = len(x)
    smax = smax or N // 8
    if smax <= smin:
        return None
    Y = np.cumsum(x - x.mean())
    scales = np.unique(np.round(np.logspace(np.log10(smin), np.log10(smax), nsc)).astype(int))
    Fq = np.full((len(qs), len(scales)), np.nan)
    for si, s in enumerate(scales):
        nseg = N // s
        if nseg < 4:
            continue
        seg = Y[: nseg * s].reshape(nseg, s)
        t = np.arange(s)
        V = np.polynomial.polynomial.polyvander(t, order)
        coef, *_ = np.linalg.lstsq(V, seg.T, rcond=None)
        F2 = np.mean((seg.T - V @ coef) ** 2, axis=0)
        F2 = F2[F2 > 0]
        if len(F2) < 4:
            continue
        for qi, q in enumerate(qs):
            if abs(q) < 1e-9:
                Fq[qi, si] = np.exp(0.5 * np.mean(np.log(F2)))
            else:
                Fq[qi, si] = np.mean(F2 ** (q / 2)) ** (1 / q)
    h = np.full(len(qs), np.nan)
    for qi in range(len(qs)):
        m = np.isfinite(Fq[qi])
        if m.sum() >= 5:
            h[qi] = np.polyfit(np.log(scales[m]), np.log(Fq[qi][m]), 1)[0]
    ok = np.isfinite(h)
    if ok.sum() < 5:
        return None
    q, hh = qs[ok], h[ok]
    tau = q * hh - 1
    alpha = np.gradient(tau, q)
    return {
        "qs": q, "h": hh, "alpha": alpha, "f": q * alpha - tau,
        "h2": float(np.interp(2, q, hh)),
        "width": float(np.nanmax(alpha) - np.nanmin(alpha)),
    }

sda

sda(x, y=None, fs: float = 25.0, maxlag: float = 10.0, nlags: int = 60) -> dict

Stabilogram diffusion analysis, after Collins and De Luca.

Postural sway behaves like two different processes at two timescales: over short intervals it drifts away from where it was, and over longer ones it is pulled back. The crossover between them is found by fitting two lines to the log-log mean-square displacement and taking their intersection.

Returns short- and long-term Hurst exponents and diffusion coefficients, and the critical time and displacement. Hs above 0.5 with Hl below it is the normal pattern: open-loop drift, then closed-loop correction.

Source code in src/micromotion/dynamics.py

def sda(x, y=None, fs: float = 25.0, maxlag: float = 10.0, nlags: int = 60) -> dict:
    """Stabilogram diffusion analysis, after Collins and De Luca.

    Postural sway behaves like two different processes at two timescales: over short
    intervals it drifts away from where it was, and over longer ones it is pulled back. The
    crossover between them is found by fitting two lines to the log-log mean-square
    displacement and taking their intersection.

    Returns short- and long-term Hurst exponents and diffusion coefficients, and the
    critical time and displacement. ``Hs`` above 0.5 with ``Hl`` below it is the normal
    pattern: open-loop drift, then closed-loop correction.
    """
    x = np.asarray(x, float)
    y = np.zeros_like(x) if y is None else np.asarray(y, float)
    n = len(x)
    L = min(int(maxlag * fs), n - 5)
    if L < 10:
        return {"Hs": float("nan"), "Hl": float("nan")}
    lags = np.unique(np.round(np.logspace(0, np.log10(L), nlags)).astype(int))
    msd = np.array([np.nanmean((x[k:] - x[:-k]) ** 2 + (y[k:] - y[:-k]) ** 2) for k in lags])
    ok = np.isfinite(msd) & (msd > 0)
    out = fit_two_region(lags[ok] / fs, msd[ok])
    out["lags_s"], out["msd"] = lags[ok] / fs, msd[ok]
    return out

fit_two_region

fit_two_region(dt, msd) -> dict

Fit two straight lines to a log-log curve and locate their crossover.

Generic: any measure with a scaling break can use it.

Source code in src/micromotion/dynamics.py

def fit_two_region(dt, msd) -> dict:
    """Fit two straight lines to a log-log curve and locate their crossover.

    Generic: any measure with a scaling break can use it.
    """
    ld, lm = np.log10(dt), np.log10(msd)
    best = None
    for i in range(4, len(dt) - 4):
        s1 = np.polyfit(ld[: i + 1], lm[: i + 1], 1)
        s2 = np.polyfit(ld[i:], lm[i:], 1)
        r = (np.sum((lm[: i + 1] - np.polyval(s1, ld[: i + 1])) ** 2)
             + np.sum((lm[i:] - np.polyval(s2, ld[i:])) ** 2))
        if best is None or r < best[0]:
            best = (r, i, s1, s2)
    if best is None:
        return {"Hs": float("nan"), "Hl": float("nan")}
    _, i, s1, s2 = best
    xc = (s2[1] - s1[1]) / (s1[0] - s2[0])
    return {"Hs": s1[0] / 2, "Hl": s2[0] / 2, "Ds": 10 ** s1[1] / 2, "Dl": 10 ** s2[1] / 2,
            "dtc": 10**xc, "msdc": 10 ** np.polyval(s1, xc), "idx": int(i)}

sampen

sampen(x, m: int = 2, r: float | None = None) -> float

Sample entropy: how unpredictable the series is.

The negative log probability that two segments which match for m samples still match for m + 1. Higher is less regular. r defaults to 0.2 standard deviations.

Source code in src/micromotion/dynamics.py

def sampen(x, m: int = 2, r: float | None = None) -> float:
    """Sample entropy: how unpredictable the series is.

    The negative log probability that two segments which match for ``m`` samples still match
    for ``m + 1``. Higher is less regular. ``r`` defaults to 0.2 standard deviations.
    """
    x = np.asarray(x, float)
    x = x[np.isfinite(x)]
    n = len(x)
    if n < m + 20:
        return float("nan")
    r = 0.2 * np.std(x) if r is None else r
    if r <= 0:
        return float("nan")

    def count(mm):
        emb = np.lib.stride_tricks.sliding_window_view(x, mm)[: n - m]
        tree = cKDTree(emb)
        return tree.count_neighbors(tree, r, p=np.inf) - len(emb)

    a, b = count(m + 1), count(m)
    return float(-np.log(a / b)) if a > 0 and b > 0 else float("nan")

ami

ami(x, maxlag: int = 100, bins: int = 32) -> np.ndarray

Average mutual information against lag, for choosing an embedding delay.

Source code in src/micromotion/dynamics.py

def ami(x, maxlag: int = 100, bins: int = 32) -> np.ndarray:
    """Average mutual information against lag, for choosing an embedding delay."""
    x = np.asarray(x, float)
    out = []
    for k in range(1, maxlag + 1):
        c, _, _ = np.histogram2d(x[:-k], x[k:], bins=bins)
        p = c / c.sum()
        px, py = p.sum(1, keepdims=True), p.sum(0, keepdims=True)
        m = p > 0
        out.append(float(np.sum(p[m] * np.log(p[m] / (px @ py)[m]))))
    return np.array(out)

first_ami_minimum

first_ami_minimum(x, maxlag: int = 100, smooth: int = 9) -> int

The conventional embedding delay: the first local minimum of mutual information.

The curve is smoothed first. Estimated from a histogram it is noisy enough that the first local minimum of the raw curve is usually spurious: on a clean sine, whose true answer is a quarter period, the unsmoothed rule returns a delay of 2.

Source code in src/micromotion/dynamics.py

def first_ami_minimum(x, maxlag: int = 100, smooth: int = 9) -> int:
    """The conventional embedding delay: the first local minimum of mutual information.

    The curve is smoothed first. Estimated from a histogram it is noisy enough that the
    first local minimum of the raw curve is usually spurious: on a clean sine, whose true
    answer is a quarter period, the unsmoothed rule returns a delay of 2.
    """
    a = ami(x, maxlag)
    if smooth > 1 and len(a) > smooth:
        k = np.ones(smooth) / smooth
        a = np.convolve(a, k, mode="same")
        a[: smooth // 2] = a[smooth // 2]
        a[-(smooth // 2):] = a[-(smooth // 2) - 1]
    for i in range(1, len(a) - 1):
        if a[i] < a[i - 1] and a[i] < a[i + 1]:
            return i + 1
    return int(np.argmin(a)) + 1

embed

embed(x, dim: int, tau: int) -> np.ndarray

Time-delay embedding. x may be one series or several columns.

Source code in src/micromotion/dynamics.py

def embed(x, dim: int, tau: int) -> np.ndarray:
    """Time-delay embedding. ``x`` may be one series or several columns."""
    x = np.asarray(x, float)
    if x.ndim == 1:
        x = x[:, None]
    n = len(x) - (dim - 1) * tau
    if n <= 0:
        raise ValueError("series is too short for this embedding")
    return np.concatenate([x[i * tau: i * tau + n] for i in range(dim)], axis=1)

rqa

rqa(x, dim: int = 3, tau: int | None = None, rr: float = 0.05, lmin: int = 2) -> dict

Recurrence quantification, at fixed recurrence rate.

The threshold is solved per plot so that exactly rr of pairs count as recurrent. That matters: with a fixed absolute threshold, determinism partly measures how tightly packed the trajectory is, so two recordings of different amplitude are not comparable.

Pass several columns to get multidimensional recurrence quantification, which is the form used for between-body coupling.

Leave tau unset unless you have a reason. A delay of 1 makes consecutive embedding vectors share all but one coordinate, so recurrences chain into diagonals that reflect the embedding rather than the dynamics: white noise embedded at tau=1 reports a determinism of 0.57, against 0.08 at a properly chosen delay.

Source code in src/micromotion/dynamics.py

def rqa(x, dim: int = 3, tau: int | None = None, rr: float = 0.05,
        lmin: int = 2) -> dict:
    """Recurrence quantification, at fixed recurrence rate.

    The threshold is solved per plot so that exactly ``rr`` of pairs count as recurrent.
    That matters: with a fixed absolute threshold, determinism partly measures how tightly
    packed the trajectory is, so two recordings of different amplitude are not comparable.

    Pass several columns to get multidimensional recurrence quantification, which is the
    form used for between-body coupling.

    Leave ``tau`` unset unless you have a reason. A delay of 1 makes consecutive embedding
    vectors share all but one coordinate, so recurrences chain into diagonals that reflect
    the embedding rather than the dynamics: white noise embedded at ``tau=1`` reports a
    determinism of 0.57, against 0.08 at a properly chosen delay.
    """
    x = np.asarray(x, float)
    if x.ndim == 1:
        tau = tau or first_ami_minimum(x, min(100, len(x) // 10))
        emb = embed(x, dim, tau)
    else:
        emb = x
        tau = tau or 1
    d = np.linalg.norm(emb[:, None, :] - emb[None, :, :], axis=2)
    eps = np.quantile(d[np.triu_indices_from(d, 1)], rr)
    R = d <= eps

    n = len(R)
    lengths = []
    for k in range(-(n - lmin), n - lmin + 1):
        if k == 0:
            continue
        diag = np.diagonal(R, k)
        run = 0
        for v in diag:
            if v:
                run += 1
            else:
                if run >= lmin:
                    lengths.append(run)
                run = 0
        if run >= lmin:
            lengths.append(run)
    lengths = np.array(lengths)
    n_rec = R.sum() - n
    det = lengths.sum() / n_rec if n_rec and len(lengths) else 0.0
    p = np.bincount(lengths)[lmin:] if len(lengths) else np.array([])
    p = p[p > 0] / p.sum() if p.sum() else np.array([])
    return {
        "RR": float(n_rec / (n * n - n)),
        "DET": float(det),
        "Lmax": int(lengths.max()) if len(lengths) else 0,
        "Lmean": float(lengths.mean()) if len(lengths) else 0.0,
        "ENTR": float(-np.sum(p * np.log(p))) if len(p) else 0.0,
        "eps": float(eps), "dim": dim, "tau": int(tau),
    }

plv

plv(a, b, fs: float, band: tuple[float, float] | None = None) -> dict

Phase-locking value between two signals.

Returns the locking strength from 0 to 1 and the preferred phase difference in radians. Band-pass first if the signals are broadband; locking is only meaningful within a band where each signal has a well-defined phase.

Source code in src/micromotion/dynamics.py

def plv(a, b, fs: float, band: tuple[float, float] | None = None) -> dict:
    """Phase-locking value between two signals.

    Returns the locking strength from 0 to 1 and the preferred phase difference in radians.
    Band-pass first if the signals are broadband; locking is only meaningful within a band
    where each signal has a well-defined phase.
    """
    from .filters import bandpass

    a = np.asarray(a, float)
    b = np.asarray(b, float)
    if band:
        a, b = bandpass(a, fs, *band), bandpass(b, fs, *band)
    pa = np.angle(_signal.hilbert(a - a.mean()))
    pb = np.angle(_signal.hilbert(b - b.mean()))
    z = np.exp(1j * (pa - pb))
    return {"plv": float(np.abs(z.mean())), "preferred_phase": float(np.angle(z.mean()))}

apen

apen(x, m: int = 2, r: float | None = None) -> float

Approximate entropy.

The older sibling of :func:sampen, and biased: it counts each template as matching itself, which pulls the estimate towards regularity, and the bias grows as the series gets shorter. Provided because the balance literature reports it and comparisons need it. For new work prefer :func:sampen, which drops the self-match.

Source code in src/micromotion/dynamics.py

def apen(x, m: int = 2, r: float | None = None) -> float:
    """Approximate entropy.

    The older sibling of :func:`sampen`, and biased: it counts each template as matching
    itself, which pulls the estimate towards regularity, and the bias grows as the series
    gets shorter. Provided because the balance literature reports it and comparisons need
    it. For new work prefer :func:`sampen`, which drops the self-match.
    """
    x = np.asarray(x, float)
    x = x[np.isfinite(x)]
    n = len(x)
    if n < m + 20:
        return float("nan")
    r = 0.2 * np.std(x) if r is None else r
    if r <= 0:
        return float("nan")

    def phi(mm):
        emb = np.lib.stride_tricks.sliding_window_view(x, mm)
        tree = cKDTree(emb)
        counts = np.array(tree.query_ball_point(emb, r, p=np.inf, return_length=True))
        return float(np.mean(np.log(counts / len(emb))))

    return float(phi(m) - phi(m + 1))

dcca

dcca(a, b, scales=None, order: int = 1) -> dict

Detrended cross-correlation between two non-stationary series.

An ordinary correlation between two signals that each wander is dominated by the wandering. This detrends both inside windows of many sizes and correlates what is left, giving a coefficient per timescale.

That per-scale answer is the point: two bodies can be uncorrelated second to second and correlated over a minute, and a single number cannot express it. Returns rho against scales, running from -1 to 1.

Source code in src/micromotion/dynamics.py

def dcca(a, b, scales=None, order: int = 1) -> dict:
    """Detrended cross-correlation between two non-stationary series.

    An ordinary correlation between two signals that each wander is dominated by the
    wandering. This detrends both inside windows of many sizes and correlates what is left,
    giving a coefficient per timescale.

    That per-scale answer is the point: two bodies can be uncorrelated second to second and
    correlated over a minute, and a single number cannot express it. Returns ``rho`` against
    ``scales``, running from -1 to 1.
    """
    a = np.asarray(a, float)
    b = np.asarray(b, float)
    n = min(len(a), len(b))
    a, b = a[:n], b[:n]
    m = np.isfinite(a) & np.isfinite(b)
    if m.sum() < 100:
        return {"scales": np.array([]), "rho": np.array([])}
    a, b = a[m], b[m]
    n = len(a)
    A, B = np.cumsum(a - a.mean()), np.cumsum(b - b.mean())
    if scales is None:
        scales = np.unique(np.round(np.logspace(np.log10(10), np.log10(n // 4), 16))
                           ).astype(int)
    rho = []
    for s in scales:
        nseg = n // s
        if nseg < 4:
            rho.append(np.nan)
            continue
        t = np.arange(s)
        V = np.polynomial.polynomial.polyvander(t, order)
        sa = A[: nseg * s].reshape(nseg, s).T
        sb = B[: nseg * s].reshape(nseg, s).T
        ra = sa - V @ np.linalg.lstsq(V, sa, rcond=None)[0]
        rb = sb - V @ np.linalg.lstsq(V, sb, rcond=None)[0]
        fab = np.mean(ra * rb)
        faa, fbb = np.mean(ra * ra), np.mean(rb * rb)
        rho.append(fab / np.sqrt(faa * fbb) if faa > 0 and fbb > 0 else np.nan)
    return {"scales": np.asarray(scales), "rho": np.asarray(rho, float)}

Postural geometry

micromotion.posture

The shape and extent of standing still.

Quantity of motion says how fast a body part moved. These say where it went: how far it strayed, what area it covered, whether it swayed along one line or wandered in all directions.

The pair is more informative than either alone. Two people with the same quantity of motion can occupy regions differing several-fold, because speed and extent are close to independent in this data — how much someone moves predicts only about a quarter of how large a region they occupy.

sway_geometry

sway_geometry(xy) -> dict

Principal axis, anisotropy and dispersion of a two-dimensional trace.

Give it the horizontal coordinates of a marker, or a centre-of-pressure track.

anisotropy runs from 0 to 1: zero when the sway is equally wide in every direction, approaching one when it collapses onto a line. There are two other definitions of this quantity in circulation in this project — one where 1 means isotropic and one reporting only the angle — so the definition is stated rather than assumed. axis_deg is the direction of greatest sway, and is axial: 10 degrees and 190 degrees describe the same posture. Test it with :func:micromotion.circular.rayleigh_axial, never with the ordinary Rayleigh test.

Source code in src/micromotion/posture.py

def sway_geometry(xy) -> dict:
    """Principal axis, anisotropy and dispersion of a two-dimensional trace.

    Give it the horizontal coordinates of a marker, or a centre-of-pressure track.

    ``anisotropy`` runs from 0 to 1: zero when the sway is equally wide in every direction,
    approaching one when it collapses onto a line. There are two other definitions of this
    quantity in circulation in this project — one where 1 means isotropic and one reporting
    only the angle — so the definition is stated rather than assumed. ``axis_deg`` is the
    direction of greatest sway, and is axial: 10 degrees and 190 degrees describe the same
    posture. Test it with :func:`micromotion.circular.rayleigh_axial`, never with the
    ordinary Rayleigh test.
    """
    p = np.asarray(xy, float)
    p = p[np.isfinite(p).all(axis=1)][:, :2]
    if len(p) < 3:
        return {"anisotropy": float("nan"), "axis_deg": float("nan")}
    c = p - p.mean(axis=0)
    _, s, vt = np.linalg.svd(c, full_matrices=False)
    lam = s**2 / max(len(c) - 1, 1)
    major, minor = float(lam[0]), float(lam[1])
    return {
        "anisotropy": float(1 - minor / major) if major > 0 else float("nan"),
        "axis_deg": float(np.degrees(np.arctan2(vt[0, 1], vt[0, 0])) % 180.0),
        "sd_major": float(np.sqrt(major)),
        "sd_minor": float(np.sqrt(minor)),
        "rms_radius": float(np.sqrt(np.mean(np.sum(c**2, axis=1)))),
        "range_major": float(np.ptp(c @ vt[0])),
        "range_minor": float(np.ptp(c @ vt[1])),
    }

ellipse_area_95

ellipse_area_95(xy) -> float

Area of the 95 per cent confidence ellipse, in the input units squared.

The standard summary of postural extent in the balance literature, which makes it the number to report when comparing against clinical work.

Source code in src/micromotion/posture.py

def ellipse_area_95(xy) -> float:
    """Area of the 95 per cent confidence ellipse, in the input units squared.

    The standard summary of postural extent in the balance literature, which makes it the
    number to report when comparing against clinical work.
    """
    p = np.asarray(xy, float)
    p = p[np.isfinite(p).all(axis=1)][:, :2]
    if len(p) < 3:
        return float("nan")
    cov = np.cov(p.T)
    return float(np.pi * _stats.chi2.ppf(0.95, 2) * np.sqrt(max(np.linalg.det(cov), 0)))

path_length

path_length(xy, fs: float | None = None) -> dict

Total distance travelled, and the rate at which it accumulated.

Unfiltered and undifferentiated, so it is not a quantity of motion and is not comparable with one: sensor noise adds to path length monotonically, which means a noisier instrument reports a longer path for an identical movement. Reported because the balance literature uses it, and because it correlates with head quantity of motion at 0.61 in this corpus, which is worth knowing but is not an equivalence.

Source code in src/micromotion/posture.py

def path_length(xy, fs: float | None = None) -> dict:
    """Total distance travelled, and the rate at which it accumulated.

    Unfiltered and undifferentiated, so it is not a quantity of motion and is not comparable
    with one: sensor noise adds to path length monotonically, which means a noisier
    instrument reports a longer path for an identical movement. Reported because the balance
    literature uses it, and because it correlates with head quantity of motion at 0.61 in
    this corpus, which is worth knowing but is not an equivalence.
    """
    p = np.asarray(xy, float)
    p = p[np.isfinite(p).all(axis=1)]
    if len(p) < 2:
        return {"path": float("nan"), "path_rate": float("nan")}
    d = float(np.sum(np.linalg.norm(np.diff(p, axis=0), axis=1)))
    return {"path": d, "path_rate": d / (len(p) / fs) if fs else float("nan")}

dispersion_radius

dispersion_radius(xy, quantile: float = 0.95) -> float

Radius containing a given proportion of the samples, about the mean position.

A robust alternative to the ellipse area when the trace has excursions: a single lean that lasts two seconds changes the ellipse considerably and this hardly at all.

Source code in src/micromotion/posture.py

def dispersion_radius(xy, quantile: float = 0.95) -> float:
    """Radius containing a given proportion of the samples, about the mean position.

    A robust alternative to the ellipse area when the trace has excursions: a single lean
    that lasts two seconds changes the ellipse considerably and this hardly at all.
    """
    p = np.asarray(xy, float)
    p = p[np.isfinite(p).all(axis=1)][:, :2]
    if len(p) < 3:
        return float("nan")
    c = p - p.mean(axis=0)
    return float(np.quantile(np.linalg.norm(c, axis=1), quantile))

Circular statistics

micromotion.circular

Statistics for angles.

Sway direction, the phase of a breath, the time of day a session happened: all are circular, and ordinary statistics quietly give wrong answers on them. The mean of 350 degrees and 10 degrees is 0, not 180.

Sway direction needs a further distinction. A body swaying forwards and backwards along one line has no preferred direction, only a preferred axis — 10 degrees and 190 degrees are the same posture. Statistics that treat those as opposite will report a person with a strong front-back sway as having no directional preference at all. The axial variants here double the angles before averaging, which is the standard fix, and they are the ones to use for sway.

circ_mean

circ_mean(angles, weights=None) -> dict

Mean direction and concentration.

Returns the mean angle in radians and R, the resultant length, which runs from 0 for angles spread uniformly around the circle to 1 for angles all identical. R is the circular answer to "how consistent is this", and there is no separate standard deviation worth reporting beside it.

Source code in src/micromotion/circular.py

def circ_mean(angles, weights=None) -> dict:
    """Mean direction and concentration.

    Returns the mean angle in radians and ``R``, the resultant length, which runs from 0 for
    angles spread uniformly around the circle to 1 for angles all identical. ``R`` is the
    circular answer to "how consistent is this", and there is no separate standard deviation
    worth reporting beside it.
    """
    a = np.asarray(angles, float)
    a = a[np.isfinite(a)]
    if not len(a):
        return {"mean": float("nan"), "R": float("nan"), "n": 0}
    w = np.ones_like(a) if weights is None else np.asarray(weights, float)[: len(a)]
    z = np.sum(w * np.exp(1j * a)) / np.sum(w)
    return {"mean": float(np.angle(z)), "R": float(np.abs(z)), "n": len(a)}

rayleigh

rayleigh(angles) -> dict

Test whether angles are spread uniformly around the circle.

The null hypothesis is no preferred direction. A small p rejects it. Uses the standard approximation, which is accurate for n above about 10.

Note what this cannot detect: a perfectly bidirectional distribution, with half the angles at 0 and half at 180, has a resultant length of zero and passes as uniform. Use :func:rayleigh_axial when the data are axial, which sway is.

Source code in src/micromotion/circular.py

def rayleigh(angles) -> dict:
    """Test whether angles are spread uniformly around the circle.

    The null hypothesis is no preferred direction. A small p rejects it. Uses the standard
    approximation, which is accurate for n above about 10.

    Note what this cannot detect: a perfectly bidirectional distribution, with half the
    angles at 0 and half at 180, has a resultant length of zero and passes as uniform. Use
    :func:`rayleigh_axial` when the data are axial, which sway is.
    """
    a = np.asarray(angles, float)
    a = a[np.isfinite(a)]
    n = len(a)
    if n < 3:
        return {"R": float("nan"), "z": float("nan"), "p": float("nan"), "n": n}
    R = np.abs(np.sum(np.exp(1j * a))) / n
    z = n * R**2
    p = np.exp(np.sqrt(1 + 4 * n + 4 * (n**2 - (n * R) ** 2)) - (1 + 2 * n))
    return {"R": float(R), "z": float(z), "p": float(min(1.0, p)), "n": n}

rayleigh_axial

rayleigh_axial(angles) -> dict

Rayleigh test for axial data, where an angle and its opposite are the same.

Doubles the angles, tests those, and halves the resulting mean direction back. This is the correct test for sway direction, marker orientation, and anything else defined on a line rather than a ray.

Source code in src/micromotion/circular.py

def rayleigh_axial(angles) -> dict:
    """Rayleigh test for axial data, where an angle and its opposite are the same.

    Doubles the angles, tests those, and halves the resulting mean direction back. This is
    the correct test for sway direction, marker orientation, and anything else defined on a
    line rather than a ray.
    """
    a = np.asarray(angles, float)
    a = a[np.isfinite(a)]
    out = rayleigh(2 * a)
    m = circ_mean(2 * a)
    out["mean_axis"] = float(m["mean"] / 2)
    return out

axial_dispersion

axial_dispersion(angles) -> float

Circular standard deviation of axial data, in radians.

Zero for a body swaying along one fixed line, rising towards the isotropic case.

Source code in src/micromotion/circular.py

def axial_dispersion(angles) -> float:
    """Circular standard deviation of axial data, in radians.

    Zero for a body swaying along one fixed line, rising towards the isotropic case.
    """
    a = np.asarray(angles, float)
    a = a[np.isfinite(a)]
    if len(a) < 2:
        return float("nan")
    R = np.abs(np.sum(np.exp(2j * a))) / len(a)
    return float(np.sqrt(-2 * np.log(max(R, 1e-12))) / 2)

circ_corr

circ_corr(a, b) -> dict

Correlation between two circular variables, after Jammalamadaka and Sengupta.

Source code in src/micromotion/circular.py

def circ_corr(a, b) -> dict:
    """Correlation between two circular variables, after Jammalamadaka and Sengupta."""
    a = np.asarray(a, float)
    b = np.asarray(b, float)
    m = np.isfinite(a) & np.isfinite(b)
    a, b = a[m], b[m]
    n = len(a)
    if n < 4:
        return {"r": float("nan"), "p": float("nan"), "n": n}
    sa = np.sin(a - circ_mean(a)["mean"])
    sb = np.sin(b - circ_mean(b)["mean"])
    r = np.sum(sa * sb) / np.sqrt(np.sum(sa**2) * np.sum(sb**2))
    l20, l02 = np.mean(sa**2), np.mean(sb**2)
    l22 = np.mean(sa**2 * sb**2)
    z = np.sqrt(n * l20 * l02 / l22) * r
    return {"r": float(r), "p": float(2 * (1 - _stats.norm.cdf(abs(z)))), "n": n}

circ_corr_linear

circ_corr_linear(angles, x) -> dict

Correlation between a circular variable and a linear one.

The question behind "does quantity of motion depend on the time of day", or on the day of the year, where the predictor wraps and the outcome does not.

Source code in src/micromotion/circular.py

def circ_corr_linear(angles, x) -> dict:
    """Correlation between a circular variable and a linear one.

    The question behind "does quantity of motion depend on the time of day", or on the day of
    the year, where the predictor wraps and the outcome does not.
    """
    a = np.asarray(angles, float)
    x = np.asarray(x, float)
    m = np.isfinite(a) & np.isfinite(x)
    a, x = a[m], x[m]
    n = len(a)
    if n < 4:
        return {"r": float("nan"), "p": float("nan"), "n": n}
    rxs = np.corrcoef(x, np.sin(a))[0, 1]
    rxc = np.corrcoef(x, np.cos(a))[0, 1]
    rcs = np.corrcoef(np.sin(a), np.cos(a))[0, 1]
    r2 = (rxc**2 + rxs**2 - 2 * rxc * rxs * rcs) / (1 - rcs**2)
    r2 = float(np.clip(r2, 0, 1))
    p = 1 - _stats.chi2.cdf(n * r2, 2)
    return {"r": float(np.sqrt(r2)), "r2": r2, "p": float(p), "n": n}

vtest

vtest(angles, mu: float) -> dict

Test for a preferred direction at a specified angle.

More powerful than :func:rayleigh when there is a prior expectation of where the preference lies — a hypothesis that sway is front-to-back, for instance, rather than a search for whichever direction happens to win.

Source code in src/micromotion/circular.py

def vtest(angles, mu: float) -> dict:
    """Test for a preferred direction at a specified angle.

    More powerful than :func:`rayleigh` when there is a prior expectation of where the
    preference lies — a hypothesis that sway is front-to-back, for instance, rather than a
    search for whichever direction happens to win.
    """
    a = np.asarray(angles, float)
    a = a[np.isfinite(a)]
    n = len(a)
    if n < 3:
        return {"V": float("nan"), "p": float("nan"), "n": n}
    R = np.abs(np.sum(np.exp(1j * a))) / n
    mean = circ_mean(a)["mean"]
    v = R * np.cos(mean - mu)
    u = v * np.sqrt(2 * n)
    return {"V": float(v), "u": float(u),
            "p": float(1 - _stats.norm.cdf(u)), "n": n}

Alignment

micromotion.align

Aligning recordings that share no clock.

The recurring problem in this corpus: two instruments recorded the same person at the same time and neither knows what time it was. A phone counts seconds from when its app started, an fNIRS headband from when its acquisition began, an audio recorder from when someone pressed the button. Putting them on one timeline means recovering the offset from the signals themselves.

Two methods, for two situations. When both signals carry the same physiological rhythm, track that rhythm in each and cross-correlate the resulting curves: :func:instantaneous_rate then :func:xcorr_lag. When the signals are of different kinds, different lengths, or sampled far apart, search integer offsets directly with :func:search_lag, which tolerates all three.

Both report how confident they are, and neither should be used without looking. A cross-correlation always has a maximum; that the maximum means anything is a separate claim.

instantaneous_rate

instantaneous_rate(x, fs: float, band: tuple[float, float] = CARDIAC_BAND, win_s: float = 20.0, step_s: float = 1.0, per_minute: bool = True)

Track a rhythm's frequency over time, by sliding a Welch estimate along it.

Returns (times, rate). With per_minute the rate is in beats or breaths per minute; otherwise in Hz.

This turns a raw signal into something two instruments can be compared on even when they measure entirely different quantities. A haemoglobin trace and a chest accelerometer have no common units, but both carry a heartbeat, and the way that heart rate wanders over ten minutes is a signature specific enough to align them.

Source code in src/micromotion/align.py

def instantaneous_rate(x, fs: float, band: tuple[float, float] = CARDIAC_BAND,
                       win_s: float = 20.0, step_s: float = 1.0,
                       per_minute: bool = True):
    """Track a rhythm's frequency over time, by sliding a Welch estimate along it.

    Returns ``(times, rate)``. With ``per_minute`` the rate is in beats or breaths per
    minute; otherwise in Hz.

    This turns a raw signal into something two instruments can be compared on even when they
    measure entirely different quantities. A haemoglobin trace and a chest accelerometer have
    no common units, but both carry a heartbeat, and the way that heart rate wanders over ten
    minutes is a signature specific enough to align them.
    """
    x = np.asarray(x, float)
    n, s = int(win_s * fs), int(step_s * fs)
    if n <= 0 or s <= 0 or len(x) < n:
        return np.array([]), np.array([])
    times, out = [], []
    nper = int(min(n, fs * 10))
    for i in range(0, len(x) - n, s):
        seg = _signal.detrend(x[i:i + n])
        f, p = _signal.welch(seg, fs, nperseg=nper)
        m = (f >= band[0]) & (f <= band[1])
        if m.any():
            out.append(f[m][np.argmax(p[m])])
            times.append((i + n / 2) / fs)
    rate = np.array(out) * (60.0 if per_minute else 1.0)
    return np.array(times), rate

xcorr_lag

xcorr_lag(a, b, fs: float = 1.0, max_lag_s: float | None = None, min_r: float = 0.5, difference: bool = True) -> dict

Offset between two equally sampled signals, by cross-correlation.

Returns the lag in seconds and samples, the correlation at that lag, and confident, which is whether that correlation reaches min_r. A positive lag means b starts later than a.

Both series are differenced first, and that default is load-bearing rather than a stylistic choice. Correlating two drifting series and taking the best lag is the classic spurious-regression trap: over 200 pairs of independent random walks, the best-lag correlation had a median of 0.47, exceeded 0.5 forty per cent of the time, and reached 0.98 at worst. No threshold can survive that. Differencing the same pairs gave a median of 0.11 and never exceeded 0.16. Differencing is shift-equivariant, so the lag itself is unaffected — it removes the drift, not the alignment. Pass difference=False only if the inputs are already stationary and you know it.

The confidence test is a threshold on the correlation and deliberately not a measure of how sharply the peak stands out from the rest of the curve. Two such measures were tried and both get it backwards: a drifting series correlates highly with itself at every nearby lag, so its true peak looks unremarkable against the curve, while white noise gives a flat curve whose maximum looks strikingly sharp. A genuine alignment scored 0.98 and pure noise 0.07, but the noise had the higher peak-to-background ratio of the two.

Source code in src/micromotion/align.py

def xcorr_lag(a, b, fs: float = 1.0, max_lag_s: float | None = None,
              min_r: float = 0.5, difference: bool = True) -> dict:
    """Offset between two equally sampled signals, by cross-correlation.

    Returns the lag in seconds and samples, the correlation at that lag, and ``confident``,
    which is whether that correlation reaches ``min_r``. A positive lag means ``b`` starts
    later than ``a``.

    Both series are differenced first, and that default is load-bearing rather than a
    stylistic choice. Correlating two drifting series and taking the best lag is the classic
    spurious-regression trap: over 200 pairs of independent random walks, the best-lag
    correlation had a median of 0.47, exceeded 0.5 forty per cent of the time, and reached
    0.98 at worst. No threshold can survive that. Differencing the same pairs gave a median
    of 0.11 and never exceeded 0.16. Differencing is shift-equivariant, so the lag itself is
    unaffected — it removes the drift, not the alignment. Pass ``difference=False`` only if
    the inputs are already stationary and you know it.

    The confidence test is a threshold on the correlation and deliberately not a measure of
    how sharply the peak stands out from the rest of the curve. Two such measures were tried
    and both get it backwards: a drifting series correlates highly with itself at every
    nearby lag, so its true peak looks unremarkable against the curve, while white noise
    gives a flat curve whose maximum looks strikingly sharp. A genuine alignment scored 0.98
    and pure noise 0.07, but the noise had the higher peak-to-background ratio of the two.
    """
    a = np.asarray(a, float)
    b = np.asarray(b, float)
    if difference and len(a) > 1 and len(b) > 1:
        a, b = np.diff(a), np.diff(b)
    a = (a - a.mean()) / (a.std() + 1e-12)
    b = (b - b.mean()) / (b.std() + 1e-12)

    c = _signal.correlate(a, b, mode="full") / min(len(a), len(b))
    lags = _signal.correlation_lags(len(a), len(b), mode="full")
    if max_lag_s is not None:
        keep = np.abs(lags) <= int(max_lag_s * fs)
        c, lags = c[keep], lags[keep]
    if not len(c):
        return {"lag_s": float("nan"), "r": float("nan"), "confident": False}

    k = int(np.argmax(c))
    peak = float(c[k])
    return {
        "lag_samples": int(lags[k]),
        "lag_s": float(lags[k]) / fs,
        "r": peak,
        "confident": bool(peak >= min_r),
    }

search_lag

search_lag(t_a, x_a, t_b, x_b, max_lag_s: float = 300.0, step_s: float = 1.0, min_overlap_s: float = 120.0, min_r: float = 0.5) -> dict

Offset between two irregular or unequal-length signals, by direct search.

Both series are put on a common grid, then every integer offset within max_lag_s is scored by Pearson correlation over whatever the two share at that offset. Slower than :func:xcorr_lag and far more tolerant: the series may differ in length, in sampling, and in how much of the session they cover.

confident is True only when the best correlation reaches min_r and the overlap reaches min_overlap_s. This recovered the audio time origin for 331 of 356 StillStanding365 days; the remaining 25 failed this test and were left blank rather than given a number nobody could trust.

Source code in src/micromotion/align.py

def search_lag(t_a, x_a, t_b, x_b, max_lag_s: float = 300.0, step_s: float = 1.0,
               min_overlap_s: float = 120.0, min_r: float = 0.5) -> dict:
    """Offset between two irregular or unequal-length signals, by direct search.

    Both series are put on a common grid, then every integer offset within ``max_lag_s`` is
    scored by Pearson correlation over whatever the two share at that offset. Slower than
    :func:`xcorr_lag` and far more tolerant: the series may differ in length, in sampling,
    and in how much of the session they cover.

    ``confident`` is True only when the best correlation reaches ``min_r`` and the overlap
    reaches ``min_overlap_s``. This recovered the audio time origin for 331 of 356
    StillStanding365 days; the remaining 25 failed this test and were left blank rather than
    given a number nobody could trust.
    """
    ga = np.arange(np.min(t_a), np.max(t_a), step_s)
    gb = np.arange(np.min(t_b), np.max(t_b), step_s)
    a = np.interp(ga, t_a, x_a)
    b = np.interp(gb, t_b, x_b)
    a = (a - a.mean()) / (a.std() + 1e-12)
    b = (b - b.mean()) / (b.std() + 1e-12)

    # The two grids start at different absolute times, so an offset in array index is not
    # the offset in seconds. Carrying the origins back in is the whole point: without it the
    # function silently reports every alignment as zero, whatever the true offset.
    origin = float(ga[0] - gb[0])

    n_min = int(min_overlap_s / step_s)
    k_max = int(max_lag_s / step_s)
    best = {"lag_s": float("nan"), "r": -np.inf, "n_overlap": 0, "confident": False}
    for k in range(-k_max, k_max + 1):
        if k >= 0:
            aa, bb = a[k:], b[: len(a) - k]
        else:
            aa, bb = a[: len(b) + k], b[-k:]
        n = min(len(aa), len(bb))
        if n < n_min:
            continue
        r = float(np.corrcoef(aa[:n], bb[:n])[0, 1])
        if np.isfinite(r) and r > best["r"]:
            best = {"lag_s": origin + k * step_s, "r": r, "n_overlap": n,
                    "confident": False}
    if not np.isfinite(best["r"]):
        return {"lag_s": float("nan"), "r": float("nan"), "n_overlap": 0,
                "confident": False}
    best["confident"] = bool(best["r"] >= min_r and best["n_overlap"] >= n_min)
    return best

find_transient

find_transient(x, fs: float, threshold: float = 8.0, search_s: float | None = None, min_separation_s: float = 1.0)

Locate impulsive events: a hand clap, a tap on a sensor, a heel strike.

Returns their times in seconds. threshold is in robust standard deviations above the median of the signal envelope, using the median absolute deviation so that the events themselves do not inflate the scale they are measured against.

The StillStanding365 sessions open and close with a synchronisation clap, and the deposited pipeline handles it by trimming a fixed twelve seconds from each end. Detecting it instead gives the offset rather than discarding it, which is what turns a clap into an alignment rather than a nuisance.

Source code in src/micromotion/align.py

def find_transient(x, fs: float, threshold: float = 8.0, search_s: float | None = None,
                   min_separation_s: float = 1.0):
    """Locate impulsive events: a hand clap, a tap on a sensor, a heel strike.

    Returns their times in seconds. ``threshold`` is in robust standard deviations above the
    median of the signal envelope, using the median absolute deviation so that the events
    themselves do not inflate the scale they are measured against.

    The StillStanding365 sessions open and close with a synchronisation clap, and the
    deposited pipeline handles it by trimming a fixed twelve seconds from each end. Detecting
    it instead gives the offset rather than discarding it, which is what turns a clap into an
    alignment rather than a nuisance.
    """
    x = np.asarray(x, float)
    if x.ndim > 1:
        x = np.linalg.norm(x, axis=1)
    env = np.abs(_signal.hilbert(x - np.median(x)))

    if search_s is not None:
        n = int(search_s * fs)
        mask = np.zeros(len(env), bool)
        mask[:n] = True
        mask[-n:] = True
    else:
        mask = np.ones(len(env), bool)

    med = np.median(env)
    mad = np.median(np.abs(env - med)) * 1.4826
    if mad <= 0:
        return np.array([])
    peaks, _ = _signal.find_peaks(
        np.where(mask, env, -np.inf),
        height=med + threshold * mad,
        distance=max(1, int(min_separation_s * fs)),
    )
    return peaks / fs

apply_lag

apply_lag(t, lag_s: float)

Shift a timebase onto another recording's origin.

Source code in src/micromotion/align.py

def apply_lag(t, lag_s: float):
    """Shift a timebase onto another recording's origin."""
    return np.asarray(t, float) + lag_s