From 841f229b73736f8f7f6cc5bc12ca3a3acbfcc23b Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 09:14:49 +0200 Subject: [PATCH 1/4] Add _posture: posturography and sway-dynamics metrics Ports the still standing study's posturography stack into pure numpy/scipy: CoP sway metrics, confidence-ellipse and convex-hull area, Collins-De Luca stabilogram diffusion, DFA, sample entropy, spectral edges, sway texture, sway orientation, axial Rayleigh, spatial extent, and principal-axis projection. From-scratch SDA/DFA/SampEn are validated against known-answer synthetic signals. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/_posture.py | 631 ++++++++++++++++++++++++++++++++++++ tests/test_posture.py | 245 ++++++++++++++ 2 files changed, 876 insertions(+) create mode 100644 musicalgestures/_posture.py create mode 100644 tests/test_posture.py diff --git a/musicalgestures/_posture.py b/musicalgestures/_posture.py new file mode 100644 index 0000000..2ecaad7 --- /dev/null +++ b/musicalgestures/_posture.py @@ -0,0 +1,631 @@ +""" +Posturography and standstill-sway metrics for centre-of-pressure (CoP) and +head/marker position signals. + +This module ports the "still standing" study's posturography stack into +pure numpy/scipy surfaces that operate on plain arrays -- no study-specific +loaders, axis conventions, or marker loops. Three families of measures are +provided: + +* **Sway amount / geometry** -- :func:`cop_sway_metrics`, + :func:`confidence_ellipse_area`, :func:`convex_hull_area`. +* **Control dynamics / complexity** -- :func:`stabilogram_diffusion` + (Collins-De Luca SDA), :func:`dfa` (detrended fluctuation analysis), + :func:`sample_entropy`, :func:`spectral_edges`, :func:`sway_texture`, + :func:`principal_axis_projection`. +* **Direction / extent** -- :func:`sway_orientation`, :func:`axial_rayleigh`, + :func:`spatial_extent`. + +The from-scratch SDA / DFA / sample-entropy implementations are validated in +the test-suite against known-answer synthetic signals (white noise -> +DFA alpha ~= 0.5 and SDA Hurst ~= 0.5; a sine -> low sample entropy relative +to its shuffle). + +Source: still standing study (Jensenius) -- posturography and micromotion +analyses of the international "standstill" championships and related datasets. +""" + +import numpy as np + + +# --------------------------------------------------------------------------- +# Sway amount / geometry +# --------------------------------------------------------------------------- +def confidence_ellipse_area(xy, conf=0.95): + """ + Area of the confidence ellipse of a 2-D point cloud (e.g. a + centre-of-pressure trace). + + The ellipse is the standard bivariate-Gaussian confidence region + ``area = pi * chi2_conf,2df * sqrt(det Cov)`` where ``Cov`` is the + 2x2 covariance of the (mean-removed) points. For a CoP sway path this + is the classic 95% "sway-ellipse area". + + Source: still standing study (Jensenius), HpSp balance analysis. + + Args: + xy (np.ndarray): Point cloud of shape ``(T, 2)``. + conf (float, optional): Confidence level in ``(0, 1)``. Defaults to + 0.95. + + Returns: + float: Ellipse area in squared position units (e.g. mm^2), or + ``nan`` if fewer than three finite points are available. + """ + from scipy.stats import chi2 + + xy = np.asarray(xy, dtype=float) + xy = xy[np.isfinite(xy).all(axis=1)] + if len(xy) < 3: + return np.nan + cov = np.cov(xy.T) + return float(np.pi * chi2.ppf(conf, 2) * np.sqrt(max(np.linalg.det(cov), 0.0))) + + +def convex_hull_area(xy): + """ + Area of the 2-D convex hull of a point cloud. + + A non-parametric alternative to :func:`confidence_ellipse_area` for the + region occupied by a sway path: it makes no Gaussian assumption and is + driven by the outermost excursions. + + Source: still standing study (Jensenius); complements the confidence + ellipse used in the balance reports. + + Args: + xy (np.ndarray): Point cloud of shape ``(T, 2)``. + + Returns: + float: Convex-hull area in squared position units, or ``nan`` if + fewer than three non-collinear finite points are available. + """ + from scipy.spatial import ConvexHull, QhullError + + xy = np.asarray(xy, dtype=float) + xy = xy[np.isfinite(xy).all(axis=1)] + if len(xy) < 3: + return np.nan + try: + return float(ConvexHull(xy).volume) # 2-D "volume" == area + except QhullError: + return np.nan + + +def cop_sway_metrics(xy, t=None, fs=None, *, freq_band=(0.1, 5.0), + resample_fs=50.0): + """ + Standard centre-of-pressure (CoP) sway metrics from a 2-D sway path. + + Computes the classic posturographic descriptors: CoP path length and + path rate, the 95% confidence-ellipse area, medio-lateral (ML) and + antero-posterior (AP) ranges and standard deviations, the AP/ML range + and SD ratios, and the mean sway frequency of each axis (the + power-weighted mean frequency of a Welch spectrum inside ``freq_band``, + computed on a uniform grid at ``resample_fs``). + + The first column of ``xy`` is treated as ML and the second as AP, + matching the study convention. Sampling time may be given either as an + explicit time vector ``t`` (seconds; may be irregular) or a constant + rate ``fs`` (Hz); if neither is supplied a rate of 1 Hz is assumed. + + Source: still standing study (Jensenius), HpSp balance analysis + (``analyze_balance``). + + Args: + xy (np.ndarray): CoP path of shape ``(T, 2)`` as ``[ML, AP]`` in + position units (e.g. mm). + t (np.ndarray, optional): Per-sample timestamps in seconds. May be + irregular. Defaults to None. + fs (float, optional): Constant sampling rate in Hz, used when ``t`` + is not given. Defaults to None (interpreted as 1 Hz). + freq_band (tuple, optional): ``(low, high)`` band in Hz for the mean + sway frequency. Defaults to ``(0.1, 5.0)``. + resample_fs (float, optional): Uniform rate in Hz onto which the + path is interpolated before the spectral estimate. Defaults to + 50.0. + + Returns: + dict: Metrics with keys ``n``, ``dur``, ``fs_mean``, ``path_len``, + ``path_rate``, ``area95``, ``ml_range``, ``ap_range``, + ``ml_sd``, ``ap_sd``, ``ap_ml_range_ratio``, + ``ap_ml_sd_ratio``, ``mf_ml``, ``mf_ap`` and ``mf_mean``. + """ + from scipy import signal as _sig + + xy = np.asarray(xy, dtype=float) + if xy.ndim != 2 or xy.shape[1] != 2: + raise ValueError("xy must have shape (T, 2)") + ml = xy[:, 0] + ap = xy[:, 1] + n = len(xy) + + if t is not None: + t = np.asarray(t, dtype=float) + t = t - t[0] + elif fs is not None: + t = np.arange(n) / float(fs) + else: + t = np.arange(n, dtype=float) + dur = float(t[-1]) if n > 1 else 0.0 + fs_mean = (n - 1) / dur if dur > 0 else np.nan + + ml_c = ml - np.nanmean(ml) + ap_c = ap - np.nanmean(ap) + + dpath = np.sqrt(np.diff(ml) ** 2 + np.diff(ap) ** 2) + path_len = float(np.nansum(dpath)) + path_rate = path_len / dur if dur > 0 else np.nan + + ml_range = float(np.nanmax(ml) - np.nanmin(ml)) + ap_range = float(np.nanmax(ap) - np.nanmin(ap)) + ml_sd = float(np.nanstd(ml, ddof=1)) + ap_sd = float(np.nanstd(ap, ddof=1)) + + area95 = confidence_ellipse_area(np.column_stack([ml_c, ap_c]), conf=0.95) + + lo, hi = freq_band + tu = np.arange(0, dur, 1.0 / resample_fs) if dur > 0 else np.array([]) + + def _mean_freq(x): + if len(tu) < 4: + return np.nan + xu = np.interp(tu, t, x) + xu = _sig.detrend(xu) + nperseg = min(len(xu), int(resample_fs * 20)) + if nperseg < 4: + return np.nan + f, P = _sig.welch(xu, fs=resample_fs, nperseg=nperseg) + band = (f >= lo) & (f <= hi) + if P[band].sum() <= 0: + return np.nan + return float(np.sum(f[band] * P[band]) / np.sum(P[band])) + + mf_ml = _mean_freq(ml_c) + mf_ap = _mean_freq(ap_c) + + return dict( + n=n, dur=dur, fs_mean=fs_mean, + path_len=path_len, path_rate=path_rate, area95=area95, + ml_range=ml_range, ap_range=ap_range, ml_sd=ml_sd, ap_sd=ap_sd, + ap_ml_range_ratio=ap_range / ml_range if ml_range else np.nan, + ap_ml_sd_ratio=ap_sd / ml_sd if ml_sd else np.nan, + mf_ml=mf_ml, mf_ap=mf_ap, mf_mean=float(np.nanmean([mf_ml, mf_ap])), + ) + + +# --------------------------------------------------------------------------- +# Control dynamics / complexity +# --------------------------------------------------------------------------- +def principal_axis_projection(xy): + """ + Project a 2-D (or N-D) point cloud onto its principal axis. + + Runs a PCA on the mean-removed points and returns the 1-D coordinate + along the direction of greatest variance -- the natural 1-D reduction of + a sway path used by the dynamics/complexity measures. + + Source: still standing study (Jensenius), sway-dynamics analysis. + + Args: + xy (np.ndarray): Point cloud of shape ``(T, D)`` (typically + ``D == 2``). + + Returns: + np.ndarray: 1-D projection of shape ``(T,)``. + """ + xy = np.asarray(xy, dtype=float) + xy = xy - xy.mean(axis=0) + cov = np.cov(xy.T) + w, v = np.linalg.eigh(cov) + return xy @ v[:, -1] # eigvecs ascending -> last is major axis + + +def stabilogram_diffusion(xy, fs, *, short_max_s=0.6, long_min_s=1.5, + n_lags=40): + """ + Collins-De Luca stabilogram-diffusion analysis (SDA) of a sway path. + + Fits the mean-square-displacement (MSD) curve + ``<[r(t+dt) - r(t)]^2>`` versus time-lag ``dt`` in log-log space and + reports a short-term and a long-term Hurst exponent (each ``slope / 2``) + plus the critical crossover time where the two regression lines + intersect. Persistent (open-loop) drift gives a short-term Hurst above + 0.5; anti-persistent (closed-loop correction) gives a long-term Hurst + below 0.5. + + Source: still standing study (Jensenius), sway-dynamics analysis; + method of Collins & De Luca (1993). + + Args: + xy (np.ndarray): Sway path of shape ``(T, D)`` (``D >= 1``). A 1-D + input of shape ``(T,)`` is accepted and treated as a single + axis. + fs (float): Sampling rate in Hz. + short_max_s (float, optional): Upper bound (s) of the short-term + fitting window. Defaults to 0.6. + long_min_s (float, optional): Lower bound (s) of the long-term + fitting window. Defaults to 1.5. + n_lags (int, optional): Number of log-spaced lags at which the MSD + is evaluated. Defaults to 40. + + Returns: + dict: ``{"H_short", "H_long", "crossover_s"}``. Entries are ``nan`` + when a window contains fewer than three usable lags. + """ + xy = np.asarray(xy, dtype=float) + if xy.ndim == 1: + xy = xy[:, None] + n = len(xy) + if n < 8: + return dict(H_short=np.nan, H_long=np.nan, crossover_s=np.nan) + + lags = np.unique( + np.round(np.logspace(0, np.log10(max(n // 4, 2)), n_lags)).astype(int)) + lags = lags[lags >= 1] + msd = np.array([np.mean(np.sum((xy[l:] - xy[:-l]) ** 2, axis=1)) + for l in lags]) + t = lags / fs + ok = msd > 0 + t, msd = t[ok], msd[ok] + short = t < short_max_s + long = t > long_min_s + + def _slope_intercept(mask): + if mask.sum() < 3: + return np.nan, np.nan + a, b = np.polyfit(np.log(t[mask]), np.log(msd[mask]), 1) + return a, b + + a_s, b_s = _slope_intercept(short) + a_l, b_l = _slope_intercept(long) + H_short = a_s / 2.0 if np.isfinite(a_s) else np.nan + H_long = a_l / 2.0 if np.isfinite(a_l) else np.nan + crossover = np.nan + if np.isfinite(a_s) and np.isfinite(a_l) and abs(a_s - a_l) > 1e-9: + crossover = float(np.exp((b_l - b_s) / (a_s - a_l))) + return dict(H_short=float(H_short), H_long=float(H_long), + crossover_s=crossover) + + +def dfa(x, *, n_scales=18, min_scale=10): + """ + Detrended fluctuation analysis (DFA) scaling exponent. + + Integrates the mean-removed signal, then measures the RMS of the + linearly-detrended integrated profile within non-overlapping windows of + increasing size; the slope of ``log F(n)`` versus ``log n`` is the DFA + exponent ``alpha``. White noise gives ``alpha ~= 0.5``; a random walk + (Brownian) gives ``alpha ~= 1.5``; ``alpha == 1`` is 1/f noise. + + Source: still standing study (Jensenius), sway-complexity analysis; + method of Peng et al. (1994). + + Args: + x (np.ndarray): 1-D input signal. + n_scales (int, optional): Number of log-spaced window sizes. + Defaults to 18. + min_scale (int, optional): Smallest window size in samples. Defaults + to 10. + + Returns: + float: The DFA exponent ``alpha`` (``nan`` if too short). + """ + x = np.asarray(x, dtype=float) + x = x[np.isfinite(x)] + n = len(x) + if n < 4 * min_scale: + return np.nan + y = np.cumsum(x - x.mean()) + scales = np.unique( + np.round(np.logspace(np.log10(min_scale), + np.log10(n // 4), n_scales)).astype(int)) + F = [] + used = [] + for m in scales: + k = n // m + if k < 1: + continue + seg = y[:k * m].reshape(k, m) + tt = np.arange(m) + rms = [np.sqrt(np.mean((s - np.polyval(np.polyfit(tt, s, 1), tt)) ** 2)) + for s in seg] + fm = np.mean(rms) + if fm > 0: + F.append(fm) + used.append(m) + if len(F) < 2: + return np.nan + return float(np.polyfit(np.log(used), np.log(F), 1)[0]) + + +def sample_entropy(x, m=2, r=0.2): + """ + Sample entropy (SampEn) of a 1-D signal. + + Measures regularity/predictability: the negative log conditional + probability that sequences close (within tolerance ``r``) for ``m`` + samples remain close for ``m + 1`` samples, self-matches excluded. + Lower values mean more repetitive/predictable signals. The signal is + z-scored internally so ``r`` is expressed as a fraction of its standard + deviation. Neighbour counts use a Chebyshev (max-norm) KD-tree. + + Source: still standing study (Jensenius), sway-complexity analysis; + method of Richman & Moorman (2000). + + Args: + x (np.ndarray): 1-D input signal. + m (int, optional): Embedding (template) length. Defaults to 2. + r (float, optional): Tolerance as a fraction of the signal SD. + Defaults to 0.2. + + Returns: + float: Sample entropy (``nan`` if undefined, e.g. no matches). + """ + from scipy.spatial import cKDTree + + x = np.asarray(x, dtype=float) + x = x[np.isfinite(x)] + N = len(x) + if N < m + 2: + return np.nan + x = (x - x.mean()) / (x.std() + 1e-12) + + def _count(mm): + emb = np.array([x[i:i + mm] for i in range(N - mm + 1)]) + tree = cKDTree(emb) + pairs = tree.query_ball_point(emb, r, p=np.inf) + return sum(len(p) - 1 for p in pairs) # exclude self-match + + B = _count(m) + A = _count(m + 1) + if B == 0 or A == 0: + return np.nan + return float(-np.log(A / B)) + + +def spectral_edges(x, fs, *, edges=(0.5, 0.95), nperseg=None): + """ + Spectral-edge frequencies of a signal. + + Returns the frequencies below which a given cumulative fraction of the + Welch power spectrum lies. With the default ``edges`` the first value is + the median frequency (50% edge) and the second the 95% spectral-edge + frequency, two standard descriptors of sway spectral shape. + + Source: still standing study (Jensenius), sway-dynamics analysis. + + Args: + x (np.ndarray): 1-D input signal. + fs (float): Sampling rate in Hz. + edges (tuple, optional): Cumulative-power fractions in ``(0, 1)``. + Defaults to ``(0.5, 0.95)``. + nperseg (int, optional): Welch segment length in samples. Defaults + to ``min(2048, len(x))``. + + Returns: + dict: Mapping ``"f"`` (e.g. ``"f50"``, ``"f95"``) to the edge + frequency in Hz. + """ + from scipy import signal as _sig + + x = np.asarray(x, dtype=float) + x = x[np.isfinite(x)] + if len(x) < 8: + return {f"f{int(round(e * 100))}": np.nan for e in edges} + if nperseg is None: + nperseg = min(2048, len(x)) + f, P = _sig.welch(x - x.mean(), fs, nperseg=nperseg) + cP = np.cumsum(P) + if cP[-1] <= 0: + return {f"f{int(round(e * 100))}": np.nan for e in edges} + cP = cP / cP[-1] + return {f"f{int(round(e * 100))}": float(np.interp(e, cP, f)) + for e in edges} + + +def sway_texture(speed, fs, *, frozen_threshold=2.0): + """ + Micro-texture of a sway speed signal: frozen fraction and burst rate. + + Distinguishes a smooth wander from intermittent ballistic corrections. + The frozen fraction is the share of time the speed is below + ``frozen_threshold``; the burst rate is the number of upward threshold + crossings (onset of a velocity burst) per minute. + + Source: still standing study (Jensenius), sway-texture analysis. + + Args: + speed (np.ndarray): 1-D speed signal (e.g. mm/s). + fs (float): Sampling rate in Hz. + frozen_threshold (float, optional): Speed below which the signal is + considered "frozen", in the units of ``speed``. Defaults to 2.0. + + Returns: + dict: ``{"frozen_fraction", "burst_rate"}`` where ``burst_rate`` is + in bursts per minute. + """ + speed = np.asarray(speed, dtype=float) + speed = speed[np.isfinite(speed)] + if len(speed) < 2: + return dict(frozen_fraction=np.nan, burst_rate=np.nan) + frozen = float((speed < frozen_threshold).mean()) + above = (speed >= frozen_threshold).astype(int) + bursts = int(np.sum(np.diff(above) == 1)) + minutes = len(speed) / fs / 60.0 + burst_rate = bursts / minutes if minutes > 0 else np.nan + return dict(frozen_fraction=frozen, burst_rate=burst_rate) + + +# --------------------------------------------------------------------------- +# Direction / extent +# --------------------------------------------------------------------------- +def sway_orientation(xy): + """ + Principal sway-axis orientation and anisotropy of a 2-D point cloud. + + A PCA of the mean-removed horizontal positions gives the orientation of + the major axis as an axial angle in ``[0, 180)`` degrees (undirected -- + a line, not an arrow) and the anisotropy ``sqrt(lambda_max / lambda_min)`` + of the sway ellipse. Anisotropy 1.0 is isotropic/circular; values above + ~1.3 indicate clearly directional sway. + + Source: still standing study (Jensenius), sway-direction analysis. + + Args: + xy (np.ndarray): Point cloud of shape ``(T, 2)``. + + Returns: + dict: ``{"angle_deg", "anisotropy"}``. Both are ``nan`` if the + covariance is degenerate or there are too few finite points. + """ + xy = np.asarray(xy, dtype=float) + xy = xy[np.isfinite(xy).all(axis=1)] + if len(xy) < 3: + return dict(angle_deg=np.nan, anisotropy=np.nan) + x = xy[:, 0] - xy[:, 0].mean() + y = xy[:, 1] - xy[:, 1].mean() + cov = np.cov(x, y) + w, v = np.linalg.eigh(cov) # ascending eigenvalues + if w[0] <= 0: + return dict(angle_deg=np.nan, anisotropy=np.nan) + angle = float(np.degrees(np.arctan2(v[1, 1], v[0, 1])) % 180.0) + anisotropy = float(np.sqrt(w[1] / w[0])) + return dict(angle_deg=angle, anisotropy=anisotropy) + + +def axial_rayleigh(angles_deg): + """ + Axial Rayleigh test for a preferred orientation among axial angles. + + Tests whether a sample of axial (undirected, ``[0, 180)`` deg) angles -- + e.g. per-session principal sway axes -- clusters around a common + orientation. Angles are doubled to map the axial circle onto the full + circle before computing the mean resultant length ``R`` and the Rayleigh + p-value (small ``p`` with large ``R`` means a shared preferred axis). + + Source: still standing study (Jensenius), sway-direction analysis. + + Args: + angles_deg (np.ndarray): Axial angles in degrees. + + Returns: + dict: ``{"R", "p", "mean_axis_deg", "n"}``. + """ + a = 2 * np.radians(np.asarray(angles_deg, dtype=float)) + a = a[np.isfinite(a)] + n = len(a) + if n < 2: + return dict(R=np.nan, p=np.nan, mean_axis_deg=np.nan, n=n) + C = np.mean(np.cos(a)) + S = np.mean(np.sin(a)) + R = float(np.hypot(C, S)) + Z = n * R * R + p = float(np.exp(-Z) * (1 + (2 * Z - Z * Z) / (4 * n))) + mean_axis = float(np.degrees(0.5 * np.arctan2(S, C)) % 180.0) + return dict(R=R, p=p, mean_axis_deg=mean_axis, n=n) + + +def spatial_extent(pos, fs, *, ellipse_conf=0.95, window_s=20.0, + vertical_axis=None): + """ + Spatial extent / occupied volume of a 3-D (or 2-D) position trace. + + Complements sway magnitude (QoM) by describing *how large a region* a + marker occupies. Reports the RMS dispersion radius about the session + centroid, the Gaussian confidence-ellipsoid volume and its cube-root + radius, the mean within-window dispersion (which removes slow drift), + and a drift ratio ``full_dispersion / within_window_dispersion`` (``> 1`` + when slow drift enlarges the occupied region over the session). When a + ``vertical_axis`` is given the drift is additionally split into + horizontal and vertical components. + + Source: still standing study (Jensenius), spatial-range analysis + (``session_metrics``). + + Args: + pos (np.ndarray): Position trace of shape ``(T, D)`` with ``D`` 2 or + 3, in position units (e.g. mm). + fs (float): Sampling rate in Hz. + ellipse_conf (float, optional): Confidence level for the ellipsoid + volume. Defaults to 0.95. + window_s (float, optional): Window length in seconds for the + within-window dispersion. Defaults to 20.0. + vertical_axis (int, optional): Index of the vertical axis (``0``, + ``1`` or ``2``); when given, drift is decomposed into horizontal + and vertical parts. Defaults to None. + + Returns: + dict: ``dispersion``, ``ellipsoid_volume``, ``ellipsoid_radius``, + ``within_window_dispersion``, ``drift_ratio`` and (when + ``vertical_axis`` is set) ``drift_horizontal`` and + ``drift_vertical``. Returns ``None`` if fewer than one window of + finite samples is available. + """ + from scipy.stats import chi2 + + pos = np.asarray(pos, dtype=float) + if pos.ndim != 2: + raise ValueError("pos must have shape (T, D)") + D = pos.shape[1] + mask = np.isfinite(pos).all(axis=1) + P = pos[mask] + w = int(fs * window_s) + if len(P) < max(w, D + 2): + return None + + centroid = P.mean(axis=0) + disp_vec = P - centroid + d = np.sqrt((disp_vec * disp_vec).sum(axis=1)) + dispersion = float(np.sqrt((d * d).mean())) + + cov = np.cov(P.T) + det = np.linalg.det(cov) + chi_d = chi2.ppf(ellipse_conf, D) + # volume of a D-dim Gaussian confidence ellipsoid + if D == 2: + vol = np.pi * chi_d * np.sqrt(max(det, 0.0)) + else: # D == 3 + vol = (4.0 / 3.0) * np.pi * (chi_d ** 1.5) * np.sqrt(max(det, 0.0)) + ellipsoid_volume = float(vol) + ellipsoid_radius = float(vol ** (1.0 / D)) + + def _within_window(Q): + cols = Q.shape[1] if Q.ndim == 2 else 1 + Q = Q.reshape(len(Q), cols) + wr = [] + for s in range(0, len(Q) - w + 1, w): + seg = Q[s:s + w] + cc = seg.mean(axis=0) + e = seg - cc + wr.append(np.sqrt((e * e).sum(axis=1).mean())) + return np.mean(wr) if wr else np.nan + + def _full_dispersion(Q): + cols = Q.shape[1] if Q.ndim == 2 else 1 + Q = Q.reshape(len(Q), cols) + cc = Q.mean(axis=0) + e = Q - cc + return np.sqrt((e * e).sum(axis=1).mean()) + + win_disp = float(_within_window(P)) + drift_ratio = (dispersion / win_disp + if win_disp and win_disp > 0 else np.nan) + + out = dict(dispersion=dispersion, ellipsoid_volume=ellipsoid_volume, + ellipsoid_radius=ellipsoid_radius, + within_window_dispersion=win_disp, drift_ratio=drift_ratio) + + if vertical_axis is not None: + hax = [k for k in range(D) if k != vertical_axis] + + def _comp_drift(cols): + Q = P[:, cols] + wm = _within_window(Q) + full = _full_dispersion(Q) + return full / wm if wm and wm > 0 else np.nan + + out["drift_horizontal"] = float(_comp_drift(hax)) + out["drift_vertical"] = float(_comp_drift([vertical_axis])) + + return out diff --git a/tests/test_posture.py b/tests/test_posture.py new file mode 100644 index 0000000..e34989d --- /dev/null +++ b/tests/test_posture.py @@ -0,0 +1,245 @@ +"""Tests for musicalgestures._posture (posturography / sway metrics). + +All ground truth is synthetic -- no data files. Circular CoP paths give +closed-form path length and ellipse area; anisotropic Gaussian clouds give a +known orientation and anisotropy; drifting vs stationary random walks give a +drift-ratio separation; white noise validates the from-scratch DFA and SDA +implementations (alpha ~= 0.5, H ~= 0.5); a sine validates that sample entropy +is low relative to its shuffle. +""" +import numpy as np +import pytest + +from musicalgestures import ( + cop_sway_metrics, + confidence_ellipse_area, + convex_hull_area, + stabilogram_diffusion, + dfa, + sample_entropy, + spectral_edges, + sway_texture, + sway_orientation, + axial_rayleigh, + spatial_extent, + principal_axis_projection, +) + + +class TestCopSwayMetrics: + def test_circular_path_length_and_ellipse(self): + # A circle of radius R sampled densely: path length -> 2*pi*R, and + # the covariance of a uniform circle is (R^2/2) I, so the 95% + # ellipse area is pi * chi2_.95,2 * R^2/2. + from scipy.stats import chi2 + R = 10.0 + fs = 100.0 + n = 4000 + theta = np.linspace(0, 2 * np.pi, n, endpoint=False) + xy = np.column_stack([R * np.cos(theta), R * np.sin(theta)]) + m = cop_sway_metrics(xy, fs=fs) + assert m["path_len"] == pytest.approx(2 * np.pi * R, rel=1e-3) + expected_area = np.pi * chi2.ppf(0.95, 2) * (R ** 2 / 2.0) + assert m["area95"] == pytest.approx(expected_area, rel=0.02) + + def test_path_rate_uses_duration(self): + R = 5.0 + fs = 50.0 + n = 2000 + theta = np.linspace(0, 2 * np.pi, n, endpoint=False) + xy = np.column_stack([R * np.cos(theta), R * np.sin(theta)]) + m = cop_sway_metrics(xy, fs=fs) + dur = (n - 1) / fs + assert m["path_rate"] == pytest.approx(2 * np.pi * R / dur, rel=1e-2) + + def test_ap_dominant_sway_ratio(self): + # AP (col 1) has larger amplitude than ML (col 0). + rng = np.random.default_rng(0) + fs = 50.0 + n = 5000 + ml = 1.0 * rng.standard_normal(n) + ap = 3.0 * rng.standard_normal(n) + m = cop_sway_metrics(np.column_stack([ml, ap]), fs=fs) + assert m["ap_ml_sd_ratio"] == pytest.approx(3.0, rel=0.1) + + def test_mean_sway_frequency(self): + # A pure 0.5 Hz oscillation on both axes -> mean freq ~ 0.5 Hz. + fs = 50.0 + t = np.arange(0, 120, 1 / fs) + f0 = 0.5 + xy = np.column_stack([np.sin(2 * np.pi * f0 * t), + np.cos(2 * np.pi * f0 * t)]) + m = cop_sway_metrics(xy, fs=fs) + assert m["mf_mean"] == pytest.approx(f0, abs=0.05) + + def test_accepts_time_vector(self): + R = 4.0 + n = 1000 + t = np.linspace(0, 20, n) + theta = np.linspace(0, 2 * np.pi, n, endpoint=False) + xy = np.column_stack([R * np.cos(theta), R * np.sin(theta)]) + m = cop_sway_metrics(xy, t=t) + assert m["dur"] == pytest.approx(20.0, rel=1e-6) + + +class TestEllipseAndHull: + def test_ellipse_area_matches_formula(self): + from scipy.stats import chi2 + rng = np.random.default_rng(1) + sx, sy = 2.0, 5.0 + xy = np.column_stack([sx * rng.standard_normal(20000), + sy * rng.standard_normal(20000)]) + area = confidence_ellipse_area(xy, conf=0.95) + expected = np.pi * chi2.ppf(0.95, 2) * sx * sy + assert area == pytest.approx(expected, rel=0.05) + + def test_convex_hull_of_square(self): + # Dense fill of a 4x6 rectangle -> hull area ~ 24. + rng = np.random.default_rng(2) + xy = np.column_stack([rng.uniform(0, 4, 5000), + rng.uniform(0, 6, 5000)]) + assert convex_hull_area(xy) == pytest.approx(24.0, rel=0.05) + + def test_degenerate_inputs(self): + assert np.isnan(confidence_ellipse_area(np.zeros((2, 2)))) + assert np.isnan(convex_hull_area(np.zeros((2, 2)))) + + +class TestSwayOrientation: + def test_known_angle_and_anisotropy(self): + # Anisotropic Gaussian cloud rotated by 30 deg: major axis 5, minor 1. + rng = np.random.default_rng(3) + n = 20000 + major = 5.0 * rng.standard_normal(n) + minor = 1.0 * rng.standard_normal(n) + phi = np.radians(30.0) + Rm = np.array([[np.cos(phi), -np.sin(phi)], + [np.sin(phi), np.cos(phi)]]) + xy = (Rm @ np.vstack([major, minor])).T + out = sway_orientation(xy) + assert out["angle_deg"] == pytest.approx(30.0, abs=2.0) + assert out["anisotropy"] == pytest.approx(5.0, rel=0.1) + + def test_isotropic_is_low_anisotropy(self): + rng = np.random.default_rng(4) + xy = rng.standard_normal((20000, 2)) + assert sway_orientation(xy)["anisotropy"] < 1.1 + + +class TestAxialRayleigh: + def test_clustered_axes_significant(self): + rng = np.random.default_rng(5) + angles = (45.0 + 3.0 * rng.standard_normal(40)) % 180.0 + out = axial_rayleigh(angles) + assert out["R"] > 0.9 + assert out["p"] < 1e-6 + assert out["mean_axis_deg"] == pytest.approx(45.0, abs=3.0) + + def test_uniform_axes_not_significant(self): + angles = np.linspace(0, 180, 60, endpoint=False) + out = axial_rayleigh(angles) + assert out["R"] < 0.1 + assert out["p"] > 0.2 + + +class TestSpatialExtent: + def test_drift_ratio_separates_drift_from_stationary(self): + rng = np.random.default_rng(6) + fs = 50.0 + n = 30000 # 600 s + # stationary: white jitter about a fixed point (no slow drift) + stat = 2.0 * rng.standard_normal((n, 3)) + # drifting: random walk that wanders over the session + steps = 0.5 * rng.standard_normal((n, 3)) + drift = np.cumsum(steps, axis=0) + m_stat = spatial_extent(stat, fs, window_s=20.0) + m_drift = spatial_extent(drift, fs, window_s=20.0) + assert m_stat["drift_ratio"] == pytest.approx(1.0, abs=0.2) + assert m_drift["drift_ratio"] > 2.0 + + def test_vertical_split(self): + rng = np.random.default_rng(7) + fs = 50.0 + n = 20000 + pos = rng.standard_normal((n, 3)) + # add a slow vertical (axis 2) drift only + pos[:, 2] += np.linspace(0, 200, n) + m = spatial_extent(pos, fs, window_s=20.0, vertical_axis=2) + assert m["drift_vertical"] > m["drift_horizontal"] + + def test_too_short_returns_none(self): + assert spatial_extent(np.zeros((10, 3)), fs=50.0) is None + + +class TestDynamicsKnownAnswers: + def test_dfa_white_noise(self): + rng = np.random.default_rng(10) + x = rng.standard_normal(20000) + alpha = dfa(x) + assert alpha == pytest.approx(0.5, abs=0.1) + + def test_dfa_random_walk(self): + rng = np.random.default_rng(11) + x = np.cumsum(rng.standard_normal(20000)) + alpha = dfa(x) + # Brownian motion -> alpha ~ 1.5 + assert alpha == pytest.approx(1.5, abs=0.15) + + def test_sda_random_walk_hurst(self): + # Collins-De Luca SDA operates on the CoP *trajectory*, modelled as a + # diffusion process. A pure random walk (Brownian trajectory) has + # MSD ~ dt^1, so the Hurst exponent (slope/2) is ~0.5 on all scales. + rng = np.random.default_rng(12) + fs = 50.0 + xy = np.cumsum(rng.standard_normal((20000, 2)), axis=0) + out = stabilogram_diffusion(xy, fs) + assert out["H_short"] == pytest.approx(0.5, abs=0.15) + assert out["H_long"] == pytest.approx(0.5, abs=0.2) + + def test_sda_white_noise_positions_flat(self): + # White-noise *positions* (no diffusion): MSD is flat vs lag, so the + # short-scale Hurst collapses toward 0 -- a useful sanity contrast to + # the random-walk case above. + rng = np.random.default_rng(15) + fs = 50.0 + xy = rng.standard_normal((20000, 2)) + out = stabilogram_diffusion(xy, fs) + assert out["H_short"] < 0.15 + + def test_sample_entropy_sine_low_vs_shuffled_high(self): + t = np.linspace(0, 40 * np.pi, 3000) + sine = np.sin(t) + rng = np.random.default_rng(13) + shuffled = sine.copy() + rng.shuffle(shuffled) + se_sine = sample_entropy(sine, m=2, r=0.2) + se_shuf = sample_entropy(shuffled, m=2, r=0.2) + assert se_sine < se_shuf + assert se_sine < 0.5 + + def test_spectral_edges_of_lowpass_signal(self): + # A 0.5 Hz sine: essentially all power at 0.5 Hz -> both edges ~0.5. + fs = 50.0 + t = np.arange(0, 200, 1 / fs) + x = np.sin(2 * np.pi * 0.5 * t) + edges = spectral_edges(x, fs) + assert edges["f50"] == pytest.approx(0.5, abs=0.1) + assert edges["f95"] == pytest.approx(0.5, abs=0.2) + + def test_sway_texture_frozen_fraction(self): + # Speed below threshold 80% of the time. + fs = 50.0 + n = 10000 + speed = np.concatenate([np.full(8000, 0.5), np.full(2000, 10.0)]) + out = sway_texture(speed, fs, frozen_threshold=2.0) + assert out["frozen_fraction"] == pytest.approx(0.8, abs=0.01) + + def test_principal_axis_projection_recovers_major_axis(self): + rng = np.random.default_rng(14) + n = 5000 + major = 10.0 * rng.standard_normal(n) + minor = 0.5 * rng.standard_normal(n) + xy = np.column_stack([major, minor]) + proj = principal_axis_projection(xy) + # projection should be almost perfectly correlated with the major axis + assert abs(np.corrcoef(proj, major)[0, 1]) > 0.99 From 730df4d0302ec0b59c71c69cfc80e8af938a7fe0 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 09:15:26 +0200 Subject: [PATCH 2/4] Add _physio: respiration rate and spectral band fractions Ports the still standing study's physiology features: windowed respiration rate (band-pass + Welch peak per window) and caller-supplied spectral band fractions (the cardiorespiratory spectral-composition diagnostic, with bands supplied by the caller so there is no sensor dependency). Registers the _posture and _physio exports in the package namespace. Validated against a 0.25 Hz breathing sinusoid (~15 br/min) and a two-tone band split. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/__init__.py | 18 +++++ musicalgestures/_physio.py | 138 ++++++++++++++++++++++++++++++++++++ tests/test_physio.py | 66 +++++++++++++++++ 3 files changed, 222 insertions(+) create mode 100644 musicalgestures/_physio.py create mode 100644 tests/test_physio.py diff --git a/musicalgestures/__init__.py b/musicalgestures/__init__.py index 03e4fe8..8739f64 100644 --- a/musicalgestures/__init__.py +++ b/musicalgestures/__init__.py @@ -69,3 +69,21 @@ def __init__(self): rayleigh_test, synchrony, ) +from musicalgestures._posture import ( + cop_sway_metrics, + confidence_ellipse_area, + convex_hull_area, + stabilogram_diffusion, + dfa, + sample_entropy, + spectral_edges, + sway_texture, + sway_orientation, + axial_rayleigh, + spatial_extent, + principal_axis_projection, +) +from musicalgestures._physio import ( + respiration_rate, + spectral_band_fractions, +) diff --git a/musicalgestures/_physio.py b/musicalgestures/_physio.py new file mode 100644 index 0000000..5b03a07 --- /dev/null +++ b/musicalgestures/_physio.py @@ -0,0 +1,138 @@ +""" +Physiology signal features for standstill / micromotion studies. + +Two pure numpy/scipy surfaces ported from the "still standing" study: + +* :func:`respiration_rate` -- windowed breathing rate (breaths per minute) + from a respiration waveform, via band-pass filtering and a Welch spectral + peak per window. +* :func:`spectral_band_fractions` -- the fraction of a signal's Welch power + falling in each of a set of caller-supplied named frequency bands. This is + the generic "cardiorespiratory QoM" spectral-composition diagnostic with + the heart-rate/respiration bands supplied by the caller, so the function + carries no dependency on any particular physiological sensor. + +Source: still standing study (Jensenius) -- Deichman / Equivital physiology +analyses. +""" + +import numpy as np + + +def respiration_rate(waveform, fs, *, band=(0.1, 0.6), window_s=30, + step_s=30): + """ + Windowed respiration rate (breaths per minute) from a breathing waveform. + + Each analysis window is band-pass filtered to the respiration band and + its dominant frequency is taken as the Welch spectral peak inside that + band; the rate is that frequency times 60. Windows advance by ``step_s`` + seconds. The default band ``(0.1, 0.6)`` Hz corresponds to about + 6-36 breaths/min. + + Source: still standing study (Jensenius), Deichman respiration analysis + (``compute_qom_resp``). + + Args: + waveform (np.ndarray): 1-D respiration/breathing waveform. + fs (float): Sampling rate in Hz. + band (tuple, optional): ``(low, high)`` respiration band in Hz. + Defaults to ``(0.1, 0.6)``. + window_s (float, optional): Window length in seconds. Defaults to 30. + step_s (float, optional): Hop between windows in seconds. Defaults to + 30. + + Returns: + dict: ``{"rate_bpm", "times_s", "median_bpm"}`` where ``rate_bpm`` is + the per-window rate (breaths/min, ``nan`` for windows without a + clear peak), ``times_s`` the window centre times in seconds, and + ``median_bpm`` the median across valid windows. + """ + from scipy.signal import butter, filtfilt, welch + + x = np.asarray(waveform, dtype=float) + x = x[np.isfinite(x)] + lo, hi = band + nyq = fs / 2.0 + if len(x) < int(fs * window_s): + # single short window: still attempt one estimate + window_s = max(1.0, len(x) / fs) + + b, a = butter(2, [lo / nyq, hi / nyq], btype="band") + xf = filtfilt(b, a, x - x.mean()) + + win = int(fs * window_s) + step = int(fs * step_s) + win = max(win, 1) + step = max(step, 1) + + rates = [] + times = [] + for start in range(0, max(len(xf) - win + 1, 1), step): + seg = xf[start:start + win] + if len(seg) < int(fs * 5): # need at least ~5 s + rates.append(np.nan) + times.append((start + win / 2) / fs) + continue + nperseg = min(len(seg), int(fs * window_s)) + f, P = welch(seg, fs, nperseg=nperseg) + mask = (f >= lo) & (f <= hi) + if mask.any() and P[mask].sum() > 0: + fpk = f[mask][np.argmax(P[mask])] + rates.append(float(fpk * 60.0)) + else: + rates.append(np.nan) + times.append((start + win / 2) / fs) + + rate_bpm = np.array(rates, dtype=float) + return dict(rate_bpm=rate_bpm, times_s=np.array(times, dtype=float), + median_bpm=float(np.nanmedian(rate_bpm)) + if np.isfinite(rate_bpm).any() else np.nan) + + +def spectral_band_fractions(signal, fs, bands, *, total_band=(0.1, 8.0), + nperseg_s=20): + """ + Fraction of a signal's power in each of a set of named frequency bands. + + Estimates the Welch power spectrum and, for each named band in ``bands``, + returns that band's summed power divided by the summed power in + ``total_band``. This is the generic spectral-composition diagnostic used + for the "cardiorespiratory QoM artifact" analysis (e.g. how much of a + chest-accelerometer QoM signal sits in a cardiac vs a respiration band), + with the bands supplied by the caller so there is no built-in dependence + on a heart-rate or respiration sensor. + + Source: still standing study (Jensenius), Deichman chest-QoM + cardiorespiratory spectral-composition analysis (``deichman_full``). + + Args: + signal (np.ndarray): 1-D input signal. + fs (float): Sampling rate in Hz. + bands (dict): Mapping of band name to ``(low, high)`` in Hz, e.g. + ``{"cardiac": (0.9, 1.3), "resp": (0.12, 0.5)}``. + total_band (tuple, optional): ``(low, high)`` reference band whose + power is the denominator. Defaults to ``(0.1, 8.0)``. + nperseg_s (float, optional): Welch segment length in seconds. + Defaults to 20. + + Returns: + dict: Mapping of each band name to its power fraction in ``[0, 1]`` + (``nan`` if the total band contains no power). + """ + from scipy.signal import welch + + x = np.asarray(signal, dtype=float) + x = x[np.isfinite(x)] + if len(x) < 8: + return {name: np.nan for name in bands} + nperseg = min(len(x), max(8, int(fs * nperseg_s))) + f, P = welch(x - x.mean(), fs, nperseg=nperseg) + tlo, thi = total_band + total = P[(f >= tlo) & (f <= thi)].sum() + if total <= 0: + return {name: np.nan for name in bands} + out = {} + for name, (lo, hi) in bands.items(): + out[name] = float(P[(f >= lo) & (f <= hi)].sum() / total) + return out diff --git a/tests/test_physio.py b/tests/test_physio.py new file mode 100644 index 0000000..1a897cb --- /dev/null +++ b/tests/test_physio.py @@ -0,0 +1,66 @@ +"""Tests for musicalgestures._physio (respiration + spectral band fractions). + +Synthetic ground truth: a 0.25 Hz breathing sinusoid must give a respiration +rate of ~15 breaths/min; a two-tone signal must split its Welch power between +two named bands in the expected proportion. +""" +import numpy as np +import pytest + +from musicalgestures import respiration_rate, spectral_band_fractions + + +class TestRespirationRate: + def test_quarter_hz_breathing_is_15_bpm(self): + fs = 25.0 + dur = 120.0 + t = np.arange(0, dur, 1 / fs) + wave = np.sin(2 * np.pi * 0.25 * t) # 0.25 Hz -> 15 br/min + out = respiration_rate(wave, fs, band=(0.1, 0.6), + window_s=30, step_s=30) + assert out["median_bpm"] == pytest.approx(15.0, abs=1.0) + assert np.isfinite(out["rate_bpm"]).all() + + def test_windows_and_times(self): + fs = 20.0 + dur = 120.0 + t = np.arange(0, dur, 1 / fs) + wave = np.sin(2 * np.pi * 0.2 * t) # 12 br/min + out = respiration_rate(wave, fs, window_s=30, step_s=30) + assert len(out["rate_bpm"]) == len(out["times_s"]) + assert out["median_bpm"] == pytest.approx(12.0, abs=1.0) + + def test_added_noise_still_recovers_rate(self): + rng = np.random.default_rng(0) + fs = 25.0 + t = np.arange(0, 120, 1 / fs) + wave = np.sin(2 * np.pi * 0.25 * t) + 0.3 * rng.standard_normal(len(t)) + out = respiration_rate(wave, fs) + assert out["median_bpm"] == pytest.approx(15.0, abs=1.5) + + +class TestSpectralBandFractions: + def test_two_tone_split(self): + # Equal-amplitude tones at 0.3 Hz and 1.2 Hz; each named band should + # capture close to half the in-band power. + fs = 50.0 + t = np.arange(0, 200, 1 / fs) + x = np.sin(2 * np.pi * 0.3 * t) + np.sin(2 * np.pi * 1.2 * t) + bands = {"low": (0.2, 0.5), "high": (1.0, 1.5)} + frac = spectral_band_fractions(x, fs, bands, total_band=(0.1, 4.0)) + assert frac["low"] == pytest.approx(0.5, abs=0.1) + assert frac["high"] == pytest.approx(0.5, abs=0.1) + assert frac["low"] + frac["high"] == pytest.approx(1.0, abs=0.1) + + def test_dominant_band_gets_most_power(self): + fs = 50.0 + t = np.arange(0, 200, 1 / fs) + # a strong 1.2 Hz tone and a weak 0.3 Hz tone + x = 0.2 * np.sin(2 * np.pi * 0.3 * t) + np.sin(2 * np.pi * 1.2 * t) + bands = {"low": (0.2, 0.5), "high": (1.0, 1.5)} + frac = spectral_band_fractions(x, fs, bands) + assert frac["high"] > frac["low"] + + def test_empty_signal(self): + frac = spectral_band_fractions(np.array([]), 50.0, {"a": (0.1, 0.5)}) + assert np.isnan(frac["a"]) From f241e8b843e55c669b8420040c2f875f57957896 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 09:15:43 +0200 Subject: [PATCH 3/4] Add _mocap: robust QTM TSV loader and cross-modality utilities Consolidates the several near-duplicate QTM/TSV loaders into one robust read_qtm_tsv (MARKER_NAMES header, autodetected data start, zeros->NaN, encoding fallback, returns marker names / [T,M,3] array / fs). Adds compare_modality_envelopes (per-second resample + Pearson r, decoupled from QoM computation) and a Welch-based dominant_frequency (not re-exported at package level to avoid colliding with _analysis.dominant_frequency). Validated with an in-test QTM TSV fixture and known-lag envelopes. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/__init__.py | 8 ++ musicalgestures/_mocap.py | 211 ++++++++++++++++++++++++++++++++++++ tests/test_mocap.py | 112 +++++++++++++++++++ 3 files changed, 331 insertions(+) create mode 100644 musicalgestures/_mocap.py create mode 100644 tests/test_mocap.py diff --git a/musicalgestures/__init__.py b/musicalgestures/__init__.py index 8739f64..19278b1 100644 --- a/musicalgestures/__init__.py +++ b/musicalgestures/__init__.py @@ -87,3 +87,11 @@ def __init__(self): respiration_rate, spectral_band_fractions, ) +# NOTE: musicalgestures._mocap also defines `dominant_frequency` (a Welch-peak +# variant from the Westney study). To avoid shadowing the pre-existing +# `_analysis.dominant_frequency` exported above, it is intentionally NOT +# re-exported here; reach it as `musicalgestures._mocap.dominant_frequency`. +from musicalgestures._mocap import ( + read_qtm_tsv, + compare_modality_envelopes, +) diff --git a/musicalgestures/_mocap.py b/musicalgestures/_mocap.py new file mode 100644 index 0000000..70d4659 --- /dev/null +++ b/musicalgestures/_mocap.py @@ -0,0 +1,211 @@ +""" +Motion-capture I/O and cross-modality utilities. + +Pure numpy/scipy helpers ported from the "still standing" and +Westney-comparisons studies: + +* :func:`read_qtm_tsv` -- a single robust reader for Qualisys Track Manager + (QTM) TSV exports, consolidating the four/five near-duplicate loaders that + were copy-pasted across the study scripts. +* :func:`compare_modality_envelopes` -- resample two motion envelopes onto a + common per-second grid and correlate them (e.g. video-pose vs mocap + validation). +* :func:`dominant_frequency` -- the dominant spectral peak of a signal within + a band, via Welch. + +.. note:: + :func:`compare_modality_envelopes` deliberately takes *precomputed* 1-D + motion envelopes rather than computing quantity-of-motion internally, so + this module stays independent of the QoM machinery. The natural producer + of such envelopes is ``musicalgestures._qom.band_limited_qom`` followed by + a per-second binning (``envelope`` / ``bin_series``), arriving in a + sibling PR. + +Source: still standing study and Westney-comparisons study (Jensenius). +""" + +import numpy as np + + +def read_qtm_tsv(path): + """ + Read a Qualisys Track Manager (QTM) TSV motion-capture export. + + Consolidates the several near-duplicate loaders used across the studies + into one robust reader. It locates the ``MARKER_NAMES`` header row to + recover marker labels, autodetects where the numeric data block starts + (the first row whose first field parses as a float), drops a trailing + all-empty column produced by a trailing tab, converts exact-zero XYZ + triples (Qualisys gap fills) to ``NaN``, and falls back from UTF-8 to + latin-1 encoding. When a ``FREQUENCY`` header field is present the frame + rate is returned as well. + + Source: still standing study and Westney-comparisons study (Jensenius) -- + unifies the ``load_qtm`` variants in the balance/dynamics/circular/ + spatial-range reports and the latin-1 variant in ``compare_mp_mocap``. + + Args: + path (str): Path to the ``.tsv`` file. + + Returns: + tuple: ``(marker_names, data, fs)`` where ``marker_names`` is a list + of ``M`` strings (empty if no header was found), ``data`` is a + float array of shape ``(T, M, 3)`` with gaps as ``NaN``, and + ``fs`` is the frame rate in Hz or ``None`` if not derivable from + the header. + """ + import io + + def _read_lines(enc): + with io.open(path, encoding=enc) as fh: + return fh.readlines() + + try: + lines = _read_lines("utf-8") + except (UnicodeDecodeError, UnicodeError): + lines = _read_lines("latin-1") + + marker_names = [] + fs = None + data_start = None + for i, ln in enumerate(lines): + parts = ln.rstrip("\n").split("\t") + key = parts[0].strip().upper() + if key == "MARKER_NAMES": + marker_names = [p.strip() for p in parts[1:] if p.strip()] + continue + if key == "FREQUENCY" and len(parts) > 1: + try: + fs = float(parts[1]) + except ValueError: + pass + continue + # first row whose leading field is numeric marks the data block + try: + float(parts[0]) + data_start = i + break + except ValueError: + continue + + if data_start is None: + raise ValueError(f"No numeric data block found in {path!r}") + + rows = [] + for ln in lines[data_start:]: + parts = ln.rstrip("\n").split("\t") + try: + rows.append([float(p) for p in parts if p != ""]) + except ValueError: + continue + if not rows: + raise ValueError(f"No parseable data rows in {path!r}") + + ncol = min(len(r) for r in rows) + arr = np.array([r[:ncol] for r in rows], dtype=float) + + # A QTM data row may carry a leading time/frame column(s); the marker + # block is the trailing 3*M columns. Prefer the marker count from the + # header; otherwise infer the largest multiple of 3. + if marker_names: + M = len(marker_names) + else: + M = ncol // 3 + offset = arr.shape[1] - 3 * M + if offset < 0: + # header lists more markers than columns present; fall back + M = arr.shape[1] // 3 + offset = arr.shape[1] - 3 * M + marker_names = marker_names[:M] + block = arr[:, offset:offset + 3 * M] + data = block.reshape(block.shape[0], M, 3) + + # exact-zero triples are Qualisys gap fills -> NaN + gap = np.all(data == 0, axis=2) + data[gap] = np.nan + + return marker_names, data, fs + + +def dominant_frequency(x, fs, band=(0.3, 4.0)): + """ + Dominant frequency of a signal within a band, via a Welch spectrum. + + Returns the frequency of the largest Welch power-spectral-density peak + inside ``band`` -- e.g. the dominant oscillation rate of a body-part + speed or vertical-position signal. + + Source: Westney-comparisons study (Jensenius), extended motion-feature + analysis (motion dominant frequency of vertical trunk position). + + Args: + x (np.ndarray): 1-D input signal. + fs (float): Sampling rate in Hz. + band (tuple, optional): ``(low, high)`` search band in Hz. Defaults + to ``(0.3, 4.0)``. + + Returns: + float: The dominant frequency in Hz, or ``nan`` if the band is empty. + """ + from scipy.signal import welch + + x = np.asarray(x, dtype=float) + x = x[np.isfinite(x)] + if len(x) < 8: + return np.nan + lo, hi = band + f, P = welch(x - x.mean(), fs, nperseg=min(2048, len(x))) + mask = (f >= lo) & (f <= hi) + if not mask.any(): + return np.nan + return float(f[mask][np.argmax(P[mask])]) + + +def compare_modality_envelopes(env_a, env_b, fs_a, fs_b): + """ + Correlate two motion envelopes after resampling to a common grid. + + Resamples both 1-D envelopes onto a shared one-sample-per-second grid + (by averaging within each second), truncates to the common length, and + returns the Pearson correlation -- the video-vs-mocap (or view-vs-view) + agreement measure. Both inputs are treated as already-computed motion + envelopes (e.g. per-frame band-limited quantity-of-motion), keeping this + function decoupled from the QoM computation itself. + + Source: still standing / Westney-comparisons study (Jensenius), + MediaPipe-vs-mocap validation (``compare_mp_mocap``). + + Args: + env_a (np.ndarray): First 1-D motion envelope. + env_b (np.ndarray): Second 1-D motion envelope. + fs_a (float): Sampling rate of ``env_a`` in Hz. + fs_b (float): Sampling rate of ``env_b`` in Hz. + + Returns: + dict: ``{"r", "n"}`` where ``r`` is the Pearson correlation of the + two per-second envelopes and ``n`` the number of common seconds + (``r`` is ``nan`` if fewer than three overlapping seconds or if + either resampled envelope is constant). + """ + def _per_second(env, fs): + env = np.asarray(env, dtype=float) + nsec = int(len(env) / fs) + if nsec < 1: + return np.array([]) + step = int(round(fs)) + step = max(step, 1) + return np.array([np.nanmean(env[i * step:(i + 1) * step]) + for i in range(nsec)]) + + a = _per_second(env_a, fs_a) + b = _per_second(env_b, fs_b) + n = min(len(a), len(b)) + if n < 3: + return dict(r=np.nan, n=n) + a = a[:n] + b = b[:n] + ok = np.isfinite(a) & np.isfinite(b) + if ok.sum() < 3 or a[ok].std() == 0 or b[ok].std() == 0: + return dict(r=np.nan, n=int(ok.sum())) + r = float(np.corrcoef(a[ok], b[ok])[0, 1]) + return dict(r=r, n=int(ok.sum())) diff --git a/tests/test_mocap.py b/tests/test_mocap.py new file mode 100644 index 0000000..9c9ba89 --- /dev/null +++ b/tests/test_mocap.py @@ -0,0 +1,112 @@ +"""Tests for musicalgestures._mocap (QTM TSV loader + cross-modality utils). + +Synthetic ground truth: an in-test QTM TSV string fixture exercises the +loader (marker names, autodetected data start, zeros->NaN); known-lag +envelopes exercise the per-second correlation; a synthetic oscillation +exercises the dominant-frequency estimator. +""" +import numpy as np +import pytest + +from musicalgestures import read_qtm_tsv, compare_modality_envelopes +from musicalgestures._mocap import dominant_frequency + + +QTM_TSV = "\t".join([ + "NO_OF_FRAMES", "4"]) + "\n" + \ + "\t".join(["NO_OF_MARKERS", "2"]) + "\n" + \ + "\t".join(["FREQUENCY", "100"]) + "\n" + \ + "\t".join(["MARKER_NAMES", "HEAD", "HAND"]) + "\n" + \ + "0\t0.0\t1.0\t2.0\t3.0\t10.0\t11.0\t12.0\n" + \ + "1\t0.01\t1.1\t2.1\t3.1\t10.1\t11.1\t12.1\n" + \ + "2\t0.02\t0.0\t0.0\t0.0\t10.2\t11.2\t12.2\n" + \ + "3\t0.03\t1.3\t2.3\t3.3\t10.3\t11.3\t12.3\n" + + +class TestReadQtmTsv: + def test_parses_fixture(self, tmp_path): + p = tmp_path / "trial.tsv" + p.write_text(QTM_TSV) + names, data, fs = read_qtm_tsv(str(p)) + assert names == ["HEAD", "HAND"] + assert fs == 100.0 + assert data.shape == (4, 2, 3) + # first frame, HEAD marker + assert np.allclose(data[0, 0], [1.0, 2.0, 3.0]) + assert np.allclose(data[0, 1], [10.0, 11.0, 12.0]) + + def test_zero_triples_become_nan(self, tmp_path): + p = tmp_path / "trial.tsv" + p.write_text(QTM_TSV) + _, data, _ = read_qtm_tsv(str(p)) + # frame 2, HEAD marker was all zeros -> NaN + assert np.all(np.isnan(data[2, 0])) + # HAND marker at frame 2 is intact + assert np.allclose(data[2, 1], [10.2, 11.2, 12.2]) + + def test_latin1_fallback(self, tmp_path): + p = tmp_path / "trial.tsv" + # a latin-1 byte (0xB5 = micro sign) in a comment-ish header line + content = QTM_TSV.replace("HEAD", "HE\xb5D") + p.write_bytes(content.encode("latin-1")) + names, data, fs = read_qtm_tsv(str(p)) + assert names[0] == "HE\xb5D" + assert data.shape == (4, 2, 3) + + def test_no_header_infers_markers(self, tmp_path): + p = tmp_path / "bare.tsv" + # no MARKER_NAMES, pure 6-column (2 marker) numeric block + p.write_text("1.0\t2.0\t3.0\t4.0\t5.0\t6.0\n" + "1.1\t2.1\t3.1\t4.1\t5.1\t6.1\n" + "1.2\t2.2\t3.2\t4.2\t5.2\t6.2\n") + names, data, fs = read_qtm_tsv(str(p)) + assert names == [] + assert data.shape == (3, 2, 3) + assert fs is None + + +class TestCompareModalityEnvelopes: + def test_identical_envelopes_correlate_perfectly(self): + rng = np.random.default_rng(0) + # 60 s of a slow envelope, sampled at two different rates + secs = 60 + base = np.abs(np.sin(np.linspace(0, 6 * np.pi, secs))) + \ + 0.1 * rng.standard_normal(secs) + fs_a, fs_b = 30.0, 100.0 + env_a = np.repeat(base, int(fs_a)) + env_b = np.repeat(base, int(fs_b)) + out = compare_modality_envelopes(env_a, env_b, fs_a, fs_b) + assert out["r"] == pytest.approx(1.0, abs=1e-6) + assert out["n"] == secs + + def test_uncorrelated_envelopes_low_r(self): + rng = np.random.default_rng(1) + fs = 30.0 + env_a = np.repeat(rng.standard_normal(60), int(fs)) + env_b = np.repeat(rng.standard_normal(60), int(fs)) + out = compare_modality_envelopes(env_a, env_b, fs, fs) + assert abs(out["r"]) < 0.4 + + def test_constant_envelope_returns_nan(self): + env_a = np.ones(600) + env_b = np.repeat(np.arange(20.0), 30) + out = compare_modality_envelopes(env_a, env_b, 30.0, 30.0) + assert np.isnan(out["r"]) + + +class TestDominantFrequency: + def test_recovers_oscillation(self): + fs = 50.0 + t = np.arange(0, 60, 1 / fs) + x = np.sin(2 * np.pi * 1.5 * t) + assert dominant_frequency(x, fs, band=(0.3, 4.0)) == \ + pytest.approx(1.5, abs=0.1) + + def test_band_excludes_out_of_band_peak(self): + fs = 50.0 + t = np.arange(0, 60, 1 / fs) + # dominant tone at 6 Hz, weak tone at 1 Hz; band (0.3,4) should + # return ~1 Hz, not 6 Hz. + x = np.sin(2 * np.pi * 6.0 * t) + 0.3 * np.sin(2 * np.pi * 1.0 * t) + assert dominant_frequency(x, fs, band=(0.3, 4.0)) == \ + pytest.approx(1.0, abs=0.2) From cbab7a3b2b1a1f13ebe7f3390f9437e5e88e7ccc Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 09:30:18 +0200 Subject: [PATCH 4/4] fix: PR2 review nits (half-open bands, 15s respiration window, doc notes) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/_mocap.py | 5 ++++- musicalgestures/_physio.py | 13 ++++++++----- musicalgestures/_posture.py | 4 +++- 3 files changed, 15 insertions(+), 7 deletions(-) diff --git a/musicalgestures/_mocap.py b/musicalgestures/_mocap.py index 70d4659..5085fe3 100644 --- a/musicalgestures/_mocap.py +++ b/musicalgestures/_mocap.py @@ -170,7 +170,10 @@ def compare_modality_envelopes(env_a, env_b, fs_a, fs_b): returns the Pearson correlation -- the video-vs-mocap (or view-vs-view) agreement measure. Both inputs are treated as already-computed motion envelopes (e.g. per-frame band-limited quantity-of-motion), keeping this - function decoupled from the QoM computation itself. + function decoupled from the QoM computation itself. The per-second binning + uses an integer-rounded step, so non-integer frame rates (e.g. 29.97 fps) + drift slightly over long signals; this function is intended for validation + rather than precise alignment. Source: still standing / Westney-comparisons study (Jensenius), MediaPipe-vs-mocap validation (``compare_mp_mocap``). diff --git a/musicalgestures/_physio.py b/musicalgestures/_physio.py index 5b03a07..79115d4 100644 --- a/musicalgestures/_physio.py +++ b/musicalgestures/_physio.py @@ -28,7 +28,8 @@ def respiration_rate(waveform, fs, *, band=(0.1, 0.6), window_s=30, its dominant frequency is taken as the Welch spectral peak inside that band; the rate is that frequency times 60. Windows advance by ``step_s`` seconds. The default band ``(0.1, 0.6)`` Hz corresponds to about - 6-36 breaths/min. + 6-36 breaths/min. Each window must contain at least 15 seconds of valid + samples for spectral estimation. Source: still standing study (Jensenius), Deichman respiration analysis (``compute_qom_resp``). @@ -70,7 +71,7 @@ def respiration_rate(waveform, fs, *, band=(0.1, 0.6), window_s=30, times = [] for start in range(0, max(len(xf) - win + 1, 1), step): seg = xf[start:start + win] - if len(seg) < int(fs * 5): # need at least ~5 s + if len(seg) < int(fs * 15): # need at least 15 s for spectral estimation rates.append(np.nan) times.append((start + win / 2) / fs) continue @@ -101,7 +102,9 @@ def spectral_band_fractions(signal, fs, bands, *, total_band=(0.1, 8.0), for the "cardiorespiratory QoM artifact" analysis (e.g. how much of a chest-accelerometer QoM signal sits in a cardiac vs a respiration band), with the bands supplied by the caller so there is no built-in dependence - on a heart-rate or respiration sensor. + on a heart-rate or respiration sensor. Power is bin-summed on the Welch + grid; the study source integrated with trapz, which yields nearly + identical results on the uniform frequency spacing of Welch. Source: still standing study (Jensenius), Deichman chest-QoM cardiorespiratory spectral-composition analysis (``deichman_full``). @@ -129,10 +132,10 @@ def spectral_band_fractions(signal, fs, bands, *, total_band=(0.1, 8.0), nperseg = min(len(x), max(8, int(fs * nperseg_s))) f, P = welch(x - x.mean(), fs, nperseg=nperseg) tlo, thi = total_band - total = P[(f >= tlo) & (f <= thi)].sum() + total = P[(f >= tlo) & (f < thi)].sum() if total <= 0: return {name: np.nan for name in bands} out = {} for name, (lo, hi) in bands.items(): - out[name] = float(P[(f >= lo) & (f <= hi)].sum() / total) + out[name] = float(P[(f >= lo) & (f < hi)].sum() / total) return out diff --git a/musicalgestures/_posture.py b/musicalgestures/_posture.py index 2ecaad7..4751bc0 100644 --- a/musicalgestures/_posture.py +++ b/musicalgestures/_posture.py @@ -203,7 +203,9 @@ def principal_axis_projection(xy): Runs a PCA on the mean-removed points and returns the 1-D coordinate along the direction of greatest variance -- the natural 1-D reduction of - a sway path used by the dynamics/complexity measures. + a sway path used by the dynamics/complexity measures. The PCA eigenvector + sign is arbitrary; the projection may be globally flipped across calls or + datasets. Source: still standing study (Jensenius), sway-dynamics analysis.