From 51a5b64e22734bc2f56841b00a1e356885531952 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 08:51:12 +0200 Subject: [PATCH 1/6] feat: canonical peak-picker (_peaks) and pulse/cycle segmentation (_pulse) Add two plain-numpy signal modules following the _analysis.py convention (standalone functions, explicit __init__ exports, lazy scipy imports): - _peaks.pick_peaks: the ONE canonical adaptive peak-picker (moving-average smoothing, relative/absolute threshold, minimum inter-peak interval, optional prominence gate) shared by all sound-motion modules in this series. Reimplemented from the cymbal-comparison paper's method prose; its default constants are documented as provisional because the paper and its deposited JSON summaries disagree on some values. - _pulse: DP-based grouping of stroke onsets into per-cycle groups (group_strokes, ported faithfully from the ro study, including the honest caveat that the carried EMA estimate is not part of the DP state), segment_cycles/Cycle/cycle_table, exponential accelerando fitting (fit_accelerando, array-based), and motion_onsets (steepest sustained rises of a motion signal), which routes through pick_peaks. Tests use synthetic ground truth only, adapting the ro study's make_ro generator (tests/_synth.py); the deep-accelerando case (25 cycles, shrinking stroke gap) requires >=95% double-group recovery. Source studies: ro; cymbal comparison (Jensenius). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/__init__.py | 11 ++ musicalgestures/_peaks.py | 88 +++++++++++ musicalgestures/_pulse.py | 297 ++++++++++++++++++++++++++++++++++++ tests/_synth.py | 83 ++++++++++ tests/test_peaks.py | 69 +++++++++ tests/test_pulse.py | 124 +++++++++++++++ 6 files changed, 672 insertions(+) create mode 100644 musicalgestures/_peaks.py create mode 100644 musicalgestures/_pulse.py create mode 100644 tests/_synth.py create mode 100644 tests/test_peaks.py create mode 100644 tests/test_pulse.py diff --git a/musicalgestures/__init__.py b/musicalgestures/__init__.py index 03e4fe8..5a8d0d3 100644 --- a/musicalgestures/__init__.py +++ b/musicalgestures/__init__.py @@ -69,3 +69,14 @@ def __init__(self): rayleigh_test, synchrony, ) + +# --- Sound--motion signal methods (ro / stillstanding / cymbal / Westney studies) --- +from musicalgestures._peaks import pick_peaks +from musicalgestures._pulse import ( + Cycle, + group_strokes, + segment_cycles, + cycle_table, + fit_accelerando, + motion_onsets, +) diff --git a/musicalgestures/_peaks.py b/musicalgestures/_peaks.py new file mode 100644 index 0000000..57805ba --- /dev/null +++ b/musicalgestures/_peaks.py @@ -0,0 +1,88 @@ +""" +Canonical adaptive peak-picking for sound and motion signals. + +This module provides the ONE peak-picker shared by the pulse, alignment, +quantity-of-motion and audio-feature modules (`_pulse`, `_alignment`, +`_qom`, `_audiofeatures`), so that every event-detection step in the +toolbox uses the same, well-tested convention: optional moving-average +smoothing, a relative (or absolute) amplitude threshold, a minimum +inter-peak interval, and an optional prominence gate. + +The function is independent of the MgVideo/MgAudio classes and operates on +any 1-D numpy signal (audio onset-detection functions, quantity-of-motion +curves, wrist-speed signals, acceleration magnitudes, ...). +""" + +import numpy as np + + +def pick_peaks(x, fs=1.0, smooth=3, rel_threshold=0.5, min_interval=0.3, + rel_prominence=0.2, threshold=None, prominence=None): + """ + Adaptive peak-picker: smoothing, relative threshold, minimum inter-peak + interval, and an optional prominence gate. + + The processing chain is: (1) an optional short moving-average smoothing + (`smooth` taps); (2) discard candidate maxima below an amplitude + threshold, expressed as a fraction of the signal's peak + (`rel_threshold`) or absolutely (`threshold`); (3) enforce a minimum + inter-peak interval of `min_interval` seconds (stronger peaks win); + (4) optionally require each peak to exceed its flanking local minima by + a prominence, again expressed as a fraction of the signal's peak + (`rel_prominence`) or absolutely (`prominence`). + + The default constants (3-tap smoothing, 0.50 x peak threshold, 0.30 s + minimum interval, 0.20 x peak prominence) are the "selective" video + quantity-of-motion settings from the cymbal-comparison study and are + PROVISIONAL defaults: that study's prose and deposited JSON summaries + disagree on some values (e.g. 0.25 x peak with a 0.10 s interval in one + deposit), so tune the parameters to the signal at hand rather than + relying on the defaults. For reference, the same study used + 0.12 x peak / 0.10 s for hand-acceleration impacts, 0.15 x peak / + 0.06 s for audio energy onsets, and 0.40 x peak / 0.20 s for + wrist-speed peaks. + + Source: cymbal-comparison study (Jensenius), reimplemented from the + paper's method description; also subsumes the peak-picking conventions + of the Westney-comparisons and ro studies. + + Args: + x (np.ndarray): Input 1-D signal. + fs (float, optional): Sampling rate of the signal (Hz). Defaults to 1.0 + (i.e. `min_interval` is then expressed in samples). + smooth (int, optional): Length of the moving-average smoothing window in + samples (taps). None, 0 or 1 disables smoothing. Defaults to 3. + rel_threshold (float, optional): Amplitude threshold as a fraction of the + (smoothed) signal's maximum. None disables the threshold. Defaults to 0.5. + min_interval (float, optional): Minimum inter-peak interval in seconds + (given `fs`). Defaults to 0.3. + rel_prominence (float, optional): Required peak prominence as a fraction of + the (smoothed) signal's maximum. None disables the gate. Defaults to 0.2. + threshold (float, optional): Absolute amplitude threshold. Overrides + `rel_threshold` when given. Defaults to None. + prominence (float, optional): Absolute prominence requirement. Overrides + `rel_prominence` when given. Defaults to None. + + Returns: + np.ndarray: Integer sample indices of the detected peaks (divide by `fs` + for times in seconds). + """ + from scipy.signal import find_peaks + from scipy.ndimage import uniform_filter1d + + x = np.asarray(x, dtype=float) + if len(x) < 3: + return np.array([], dtype=int) + + if smooth is not None and smooth > 1: + x = uniform_filter1d(x, size=int(smooth)) + + peak = float(np.max(x)) + height = threshold if threshold is not None else ( + rel_threshold * peak if rel_threshold is not None else None) + prom = prominence if prominence is not None else ( + rel_prominence * peak if rel_prominence is not None else None) + distance = max(1, int(round(min_interval * fs))) if min_interval else 1 + + idx, _ = find_peaks(x, height=height, distance=distance, prominence=prom) + return idx diff --git a/musicalgestures/_pulse.py b/musicalgestures/_pulse.py new file mode 100644 index 0000000..3cb984f --- /dev/null +++ b/musicalgestures/_pulse.py @@ -0,0 +1,297 @@ +""" +Pulse and cycle segmentation for accelerating rhythmic sequences. + +Tools for grouping event onsets (e.g. drum strokes) into per-cycle stroke +groups, tabulating per-cycle metrics, fitting an exponential accelerando +model, and detecting motion onsets from a quantity-of-motion signal. + +These helpers are independent of the MgVideo/MgAudio classes and operate on +plain numpy arrays of onset times or 1-D motion signals. + +Source: ro study (Jensenius) -- analysis of the accelerating ro ritual's +drum-stroke cycles and their coupling to body motion. +""" + +from dataclasses import dataclass + +import numpy as np + +from musicalgestures._peaks import pick_peaks + + +@dataclass +class Cycle: + """ + One rhythmic cycle: a group of stroke onsets plus an optional secondary + event (e.g. a shout) that falls inside the cycle. + + Source: ro study (Jensenius). + + Attributes: + index (int): Zero-based cycle index. + strokes (list): Onset times (s) of the strokes in this cycle. + event (float | None): Time (s) of the cycle's secondary event, or None. + """ + index: int + strokes: list + event: float | None = None + + @property + def t_start(self): + """float: Time (s) of the cycle's first stroke.""" + return self.strokes[0] + + @property + def n_strokes(self): + """int: Number of strokes in the cycle.""" + return len(self.strokes) + + @property + def stroke_gap(self): + """float: Gap (s) between the first two strokes (NaN if fewer than 2).""" + return self.strokes[1] - self.strokes[0] if len(self.strokes) >= 2 else np.nan + + +def group_strokes(onset_times, max_strokes=4, gap_lo=0.10, gap_hi=0.60, + w_abs=6.0, w_within=4.0, tol_within=0.15, + w_order=2.5, w_trend=2.0, tol_trend=0.25, + size_costs=(1.0, 0.0, 1.5, 6.0)): + """ + Segment stroke onsets into per-cycle stroke groups by dynamic + programming over candidate group boundaries (Viterbi over segmentations), + with a greedily carried stroke-gap estimate (EMA) that is exact given + that carried estimate but is not itself part of the DP state key, so + Bellman optimality does not strictly hold over the full segmentation. + + The cost of a segmentation encodes structural priors for an accelerating + cyclic pattern (developed for the ro ritual's double drum strokes): + + * size prior -- `size_costs[k-1]` per k-stroke group: with the defaults, + double strokes are free, singles/triples carry a small penalty, >=4 a + steep one; + * within-gap plausibility -- gaps inside a group should fall in + `[gap_lo, gap_hi]` seconds (`w_abs` x log-excess outside) and stay + close to a running stroke-gap estimate (EMA, weight 0.5): `w_within` + x |log-ratio| beyond `tol_within`; + * ordering prior -- a between-group gap shorter than the current + stroke-gap estimate costs `w_order` x the log-ratio shortfall; + * accelerando prior -- successive between-group gaps should not grow: + an increase beyond `tol_trend` (log) costs `w_trend` x the excess + (decreases are free, so a climax's shrinking gaps cost nothing). + + Unlike a single running threshold, the trend and ordering terms let the + decision boundary between stroke gaps and cycle gaps shrink with the + accelerando, so the climax stays resolved even when the cycle gap drops + to (or just below) the stroke gap in the final cycles. + + Source: ro study (Jensenius). + + Args: + onset_times (np.ndarray): Stroke onset times in seconds (any order). + max_strokes (int, optional): Maximum strokes per group. Defaults to 4. + gap_lo (float, optional): Lower bound (s) of plausible within-group gaps. Defaults to 0.10. + gap_hi (float, optional): Upper bound (s) of plausible within-group gaps. Defaults to 0.60. + w_abs (float, optional): Weight of the absolute within-gap plausibility term. + Defaults to 6.0. + w_within (float, optional): Weight of the within-gap-vs-EMA term. Defaults to 4.0. + tol_within (float, optional): Log-ratio tolerance of the within-gap-vs-EMA term. + Defaults to 0.15. + w_order (float, optional): Weight of the ordering prior. Defaults to 2.5. + w_trend (float, optional): Weight of the accelerando (non-increasing gaps) prior. + Defaults to 2.0. + tol_trend (float, optional): Log-ratio tolerance of the accelerando prior. Defaults to 0.25. + size_costs (tuple, optional): Cost per group of size 1, 2, 3, >=4. + Defaults to (1.0, 0.0, 1.5, 6.0). + + Returns: + list: A list of groups, each a list of stroke onset times (floats, seconds). + """ + t = np.sort(np.asarray(onset_times, float)) + n = len(t) + if n == 0: + return [] + if n == 1: + return [[float(t[0])]] + log = np.log + gaps = np.maximum(np.diff(t), 1e-3) + + def size_cost(k): + return size_costs[min(k, len(size_costs)) - 1] + + def within_cost(w, ema): + c = w_abs * (max(0.0, log(w / gap_hi)) + max(0.0, log(gap_lo / w))) + if ema is not None: + c += w_within * max(0.0, abs(log(w / ema)) - tol_within) + return c + + # State (i, j): a completed group spans onsets j..i. Value: cost, the + # between-gap that preceded the group, the stroke-gap EMA, backpointer. + best = {} + + def push(key, cost, prev_g, ema, back): + cur = best.get(key) + if cur is None or cost < cur[0]: + best[key] = (cost, prev_g, ema, back) + + for i in range(min(n, max_strokes)): # groups starting at onset 0 + cost, ema = size_cost(i + 1), None + for k in range(i): + cost += within_cost(gaps[k], ema) + ema = gaps[k] if ema is None else 0.5 * ema + 0.5 * gaps[k] + push((i, 0), cost, None, ema, None) + + for i in range(n - 1): + for j in range(max(0, i - max_strokes + 1), i + 1): + state = best.get((i, j)) + if state is None: + continue + cost0, prev_g, ema0, _ = state + g = gaps[i] # between-group gap + cost0 += w_trend * max(0.0, log(g / prev_g) - tol_trend) \ + if prev_g is not None else 0.0 + if ema0 is not None: + cost0 += w_order * max(0.0, log(ema0 / g)) + for i2 in range(i + 1, min(n, i + 1 + max_strokes)): + size = i2 - i + cost, ema = cost0 + size_cost(size), ema0 + for k in range(i + 1, i2): + cost += within_cost(gaps[k], ema) + ema = gaps[k] if ema is None else 0.5 * ema + 0.5 * gaps[k] + push((i2, i + 1), cost, g, ema, (i, j)) + + key = min((k for k in best if k[0] == n - 1), key=lambda k: best[k][0]) + starts = [] + while key is not None: + starts.append(key[1]) + key = best[key][3] + starts.reverse() + return [list(map(float, t[a:b])) + for a, b in zip(starts, starts[1:] + [n])] + + +def segment_cycles(onset_times, event_times=None, **kwargs): + """ + Segment stroke onsets into `Cycle` objects: one Cycle per stroke group + (see `group_strokes`); the cycle's secondary event is the first event + onset in [group start, next group start). + + Source: ro study (Jensenius) -- the secondary events were the ritual's + shouts. + + Args: + onset_times (np.ndarray): Stroke onset times in seconds. + event_times (np.ndarray, optional): Onset times (s) of a secondary event + stream (e.g. shouts) to assign to cycles. Defaults to None. + **kwargs: Passed on to `group_strokes`. + + Returns: + list: A list of `Cycle` objects. + """ + groups = group_strokes(onset_times, **kwargs) + event_times = np.sort(np.asarray( + [] if event_times is None else event_times, float)) + starts = [g[0] for g in groups] + [np.inf] + out = [] + for i, g in enumerate(groups): + ev = event_times[(event_times >= starts[i]) & (event_times < starts[i + 1])] + out.append(Cycle(i, list(g), float(ev[0]) if len(ev) else None)) + return out + + +def cycle_table(cycles, clip_id="", context=""): + """ + Tabulate per-cycle metrics from a list of `Cycle` objects. + + Source: ro study (Jensenius). + + Args: + cycles (list): A list of `Cycle` objects (see `segment_cycles`). + clip_id (str, optional): Identifier written to the `clip` column. Defaults to "". + context (str, optional): Label written to the `context` column. Defaults to "". + + Returns: + pd.DataFrame: One row per cycle with columns `clip`, `context`, `cycle`, + `t` (cycle start, s), `ioi` (to next cycle start, s), `n_strokes`, + `stroke_gap` (s), `event` (secondary-event time, s or NaN) and + `stroke_event_ioi` (event minus cycle start, s or NaN). + """ + import pandas as pd + columns = ['clip', 'context', 'cycle', 't', 'ioi', 'n_strokes', + 'stroke_gap', 'event', 'stroke_event_ioi'] + rows = [] + for i, c in enumerate(cycles): + nxt = cycles[i + 1].t_start if i + 1 < len(cycles) else np.nan + rows.append(dict( + clip=clip_id, context=context, cycle=c.index, t=c.t_start, + ioi=nxt - c.t_start, n_strokes=c.n_strokes, + stroke_gap=c.stroke_gap, event=c.event, + stroke_event_ioi=(c.event - c.t_start) if c.event is not None else np.nan, + )) + if not rows: + return pd.DataFrame(columns=columns) + return pd.DataFrame(rows) + + +def fit_accelerando(times, iois): + """ + Fit an exponential accelerando model IOI(t) = ioi0 * 2**(-t / t_double) + via least squares on log2(IOI): `t_double` is the time (s) it takes the + inter-onset interval to halve (i.e. the tempo to double). + + Source: ro study (Jensenius). + + Args: + times (np.ndarray): Cycle start times in seconds. + iois (np.ndarray): Inter-onset intervals (s) at those times. Non-finite + or non-positive entries are ignored. + + Returns: + tuple: `(ioi0, t_double, r2)` where `ioi0` is the fitted IOI at t=0 (s), + `t_double` is the tempo-doubling time (s; `np.inf` if the sequence is + not accelerating), and `r2` is the fit's coefficient of determination + on log2(IOI). + """ + t = np.asarray(times, float) + ioi = np.asarray(iois, float) + m = np.isfinite(t) & np.isfinite(ioi) & (ioi > 0) + t, ioi = t[m], np.log2(ioi[m]) + A = np.vstack([t, np.ones_like(t)]).T + (slope, intercept), *_ = np.linalg.lstsq(A, ioi, rcond=None) + pred = A @ [slope, intercept] + ss_tot = ((ioi - ioi.mean()) ** 2).sum() + r2 = 1 - ((ioi - pred) ** 2).sum() / ss_tot if ss_tot > 0 else np.nan + return 2.0 ** intercept, (-1.0 / slope if slope < 0 else np.inf), r2 + + +def motion_onsets(motion, fs, min_interval=0.25, smooth_cutoff=8.0): + """ + Times of the steepest sustained rises in a motion signal (e.g. a + quantity-of-motion curve): the signal is low-pass filtered, its positive + time-derivative is formed, and peaks of that derivative are picked with + the canonical peak-picker (`musicalgestures.pick_peaks`) using a robust + prominence gate of 0.25 x (99th percentile - median) of the derivative. + + Source: ro study (Jensenius) -- motion onsets of the rowing gesture, + related to the drum cycles via `per_cycle_motion_delta`. + + Args: + motion (np.ndarray): 1-D motion signal (e.g. mean absolute frame difference). + fs (float): Sampling rate of the signal (Hz, e.g. video frames per second). + min_interval (float, optional): Minimum interval between onsets (s). Defaults to 0.25. + smooth_cutoff (float, optional): Low-pass cutoff (Hz) applied before + differentiation. Defaults to 8.0. + + Returns: + np.ndarray: Onset times in seconds. + """ + from scipy.signal import butter, filtfilt + motion = np.asarray(motion, float) + cutoff = min(smooth_cutoff, 0.45 * fs) + b, a = butter(3, cutoff / (fs / 2)) + s = filtfilt(b, a, motion) + d = np.clip(np.gradient(s) * fs, 0, None) + prom = (np.percentile(d, 99) - np.median(d)) * 0.25 + idx = pick_peaks(d, fs=fs, smooth=None, rel_threshold=None, + min_interval=min_interval, rel_prominence=None, + prominence=prom) + return idx / fs diff --git a/tests/_synth.py b/tests/_synth.py new file mode 100644 index 0000000..0eeba50 --- /dev/null +++ b/tests/_synth.py @@ -0,0 +1,83 @@ +"""Synthetic ground-truth generators for the sound--motion analysis tests. + +Adapted from the ro study's test synthesizer (Jensenius): accelerating +double-stroke pulse trains with known cycle starts, click/burst audio +rendering, and simple helpers for click trains and decaying tones. +No media fixtures are needed: every test signal is generated here. +""" +import numpy as np + +SR = 22050 + + +def click(sr=SR, f=180.0, dur=0.03): + """A short decaying sine click.""" + t = np.arange(int(dur * sr)) / sr + return (np.sin(2 * np.pi * f * t) * np.exp(-t * 60)).astype("float32") + + +def shout_burst(sr=SR, dur=0.3, f0=300.0): + """Vowel-like harmonic stack (300-2700 Hz), like a crowd shout.""" + t = np.arange(int(dur * sr)) / sr + y = sum(np.sin(2 * np.pi * f0 * k * t) / k for k in range(1, 10)) + y = y / np.abs(y).max() + return (0.8 * y * np.hanning(len(t))).astype("float32") + + +def ro_times(ioi0=2.0, t_double=12.0, n_cycles=15, stroke_gap=0.25, + shout_frac=0.5, gap_shrink=0.0): + """Ground-truth event times (no audio) for an accelerating ro sequence: + IOI(t) = ioi0 * 2**(-(t - 1.0) / t_double); each cycle is a double drum + stroke plus one shout at shout_frac * current IOI. gap_shrink linearly + shrinks the stroke gap to (1 - gap_shrink) * stroke_gap by the last cycle. + gt['stroke_gap'] stays the initial gap for backward compatibility.""" + starts, t = [], 1.0 + for _ in range(n_cycles): + starts.append(t) + t += ioi0 * 2 ** (-(t - 1.0) / t_double) + gt = {"starts": starts, "strokes": [], "shouts": [], + "stroke_gap": stroke_gap} + for k, s in enumerate(starts): + ioi_here = ioi0 * 2 ** (-(s - 1.0) / t_double) + frac = k / max(1, n_cycles - 1) + gap_here = stroke_gap * (1.0 - gap_shrink * frac) + gt["strokes"] += [s, s + gap_here] + gt["shouts"].append(s + shout_frac * ioi_here) + return gt + + +def make_ro(ioi0=2.0, t_double=12.0, n_cycles=15, stroke_gap=0.25, + shout_frac=0.5, sr=SR, gap_shrink=0.0): + """Render ro_times() to audio: click per stroke, vowel burst per shout.""" + gt = ro_times(ioi0, t_double, n_cycles, stroke_gap, shout_frac, gap_shrink) + total = gt["starts"][-1] + 2.0 + y = np.zeros(int(total * sr), "float32") + c, b = click(sr), shout_burst(sr) + for st in gt["strokes"]: + i0 = int(st * sr) + y[i0:i0 + len(c)] += c + for sh in gt["shouts"]: + i0 = int(sh * sr) + y[i0:i0 + len(b)] += b + return y / max(1e-9, np.abs(y).max()), gt + + +def click_train(times, sr=SR, tail=0.5, **click_kw): + """Render a click at each time (s); returns the waveform.""" + times = np.asarray(times, float) + c = click(sr, **click_kw) + y = np.zeros(int((times.max() + tail) * sr), "float32") + for t in times: + i0 = int(t * sr) + y[i0:i0 + len(c)] += c + return y / max(1e-9, np.abs(y).max()) + + +def decaying_tone(t60, sr=SR, f=440.0, dur=None, onset=0.05): + """A tone with an exact exponential decay of the given T60 (s).""" + if dur is None: + dur = onset + 0.8 * t60 + n = int(dur * sr) + t = np.arange(n) / sr + envelope = np.where(t < onset, t / onset, 10 ** (-3 * (t - onset) / t60)) + return (np.sin(2 * np.pi * f * t) * envelope).astype("float32") diff --git a/tests/test_peaks.py b/tests/test_peaks.py new file mode 100644 index 0000000..1634a0e --- /dev/null +++ b/tests/test_peaks.py @@ -0,0 +1,69 @@ +"""Tests for musicalgestures._peaks (canonical adaptive peak-picker).""" +import numpy as np + +from musicalgestures import pick_peaks + + +def bumpy_signal(fs=100.0, peak_times=(1.0, 2.0, 3.5), amps=(1.0, 0.6, 0.9), + dur=5.0, width=0.05): + t = np.arange(int(dur * fs)) / fs + x = np.zeros_like(t) + for pt, a in zip(peak_times, amps): + x += a * np.exp(-0.5 * ((t - pt) / width) ** 2) + return t, x + + +class TestPickPeaks: + def test_recovers_known_peaks(self): + fs = 100.0 + t, x = bumpy_signal(fs) + idx = pick_peaks(x, fs=fs, rel_threshold=0.3, min_interval=0.3, + rel_prominence=0.2) + assert len(idx) == 3 + assert np.allclose(idx / fs, [1.0, 2.0, 3.5], atol=0.05) + + def test_relative_threshold_gates_small_peaks(self): + fs = 100.0 + t, x = bumpy_signal(fs, amps=(1.0, 0.2, 0.9)) + idx = pick_peaks(x, fs=fs, rel_threshold=0.5, min_interval=0.1, + rel_prominence=None) + assert np.allclose(idx / fs, [1.0, 3.5], atol=0.05) + + def test_min_interval_keeps_stronger_peak(self): + fs = 100.0 + t, x = bumpy_signal(fs, peak_times=(1.0, 1.15), amps=(1.0, 0.8)) + idx = pick_peaks(x, fs=fs, rel_threshold=0.3, min_interval=0.5, + rel_prominence=None, smooth=None) + assert len(idx) == 1 + assert abs(idx[0] / fs - 1.0) < 0.05 + + def test_prominence_gate_rejects_shoulder(self): + fs = 100.0 + t = np.arange(int(5 * fs)) / fs + # A big slow hump with a tiny ripple riding on it: the ripple's peak + # is high in absolute terms but has almost no prominence. + x = np.exp(-0.5 * ((t - 2.5) / 1.0) ** 2) + x += 0.05 * np.exp(-0.5 * ((t - 3.5) / 0.02) ** 2) + idx_gated = pick_peaks(x, fs=fs, smooth=None, rel_threshold=0.3, + min_interval=0.2, rel_prominence=0.2) + idx_open = pick_peaks(x, fs=fs, smooth=None, rel_threshold=0.3, + min_interval=0.2, rel_prominence=None) + assert len(idx_gated) == 1 + assert len(idx_open) == 2 + + def test_absolute_overrides(self): + fs = 100.0 + t, x = bumpy_signal(fs, amps=(1.0, 0.6, 0.9)) + idx = pick_peaks(x, fs=fs, rel_threshold=None, threshold=0.7, + min_interval=0.1, rel_prominence=None) + assert np.allclose(idx / fs, [1.0, 3.5], atol=0.05) + + def test_empty_and_short_input(self): + assert len(pick_peaks(np.array([]))) == 0 + assert len(pick_peaks(np.array([1.0, 2.0]))) == 0 + + def test_returns_integer_indices(self): + fs = 100.0 + t, x = bumpy_signal(fs) + idx = pick_peaks(x, fs=fs) + assert idx.dtype.kind == "i" diff --git a/tests/test_pulse.py b/tests/test_pulse.py new file mode 100644 index 0000000..8c80f75 --- /dev/null +++ b/tests/test_pulse.py @@ -0,0 +1,124 @@ +"""Tests for musicalgestures._pulse (cycle segmentation, accelerando, motion onsets). + +All ground truth is synthetic (see tests/_synth.py, adapted from the ro study). +""" +import numpy as np +import pytest + +from musicalgestures import ( + Cycle, + group_strokes, + segment_cycles, + cycle_table, + fit_accelerando, + motion_onsets, +) +from _synth import ro_times + + +class TestGroupStrokes: + def test_basic_double_strokes(self): + gt = ro_times(n_cycles=15) + groups = group_strokes(gt["strokes"]) + assert len(groups) == 15 + assert all(len(g) == 2 for g in groups) + assert np.allclose([g[0] for g in groups], gt["starts"]) + + def test_deep_accelerando_25_cycles(self): + # Deep-accelerando case: 25 cycles with a shrinking stroke gap; the + # final cycle gaps approach the stroke gap. At least 95% of the + # groups must be recovered as double strokes at the right starts. + gt = ro_times(ioi0=2.0, t_double=10.0, n_cycles=25, + stroke_gap=0.25, gap_shrink=0.4) + groups = group_strokes(gt["strokes"]) + doubles = [g for g in groups if len(g) == 2] + assert len(groups) == 25 + assert len(doubles) / len(groups) >= 0.95 + recovered = sorted(g[0] for g in groups) + assert np.allclose(recovered, gt["starts"], atol=1e-9) + + def test_empty_and_single(self): + assert group_strokes([]) == [] + assert group_strokes([3.0]) == [[3.0]] + + def test_unsorted_input(self): + gt = ro_times(n_cycles=8) + shuffled = np.asarray(gt["strokes"])[::-1] + groups = group_strokes(shuffled) + assert len(groups) == 8 + + +class TestSegmentCycles: + def test_assigns_first_event_in_cycle(self): + gt = ro_times(n_cycles=10) + cycles = segment_cycles(gt["strokes"], gt["shouts"]) + assert len(cycles) == 10 + assert all(isinstance(c, Cycle) for c in cycles) + got = [c.event for c in cycles] + assert np.allclose(got, gt["shouts"]) + + def test_no_events(self): + gt = ro_times(n_cycles=5) + cycles = segment_cycles(gt["strokes"]) + assert all(c.event is None for c in cycles) + + def test_cycle_properties(self): + c = Cycle(0, [1.0, 1.25], 1.8) + assert c.t_start == 1.0 + assert c.n_strokes == 2 + assert c.stroke_gap == pytest.approx(0.25) + assert np.isnan(Cycle(0, [1.0]).stroke_gap) + + +class TestCycleTable: + def test_columns_and_values(self): + gt = ro_times(n_cycles=6) + df = cycle_table(segment_cycles(gt["strokes"], gt["shouts"]), + clip_id="clip1", context="test") + assert list(df.columns) == ['clip', 'context', 'cycle', 't', 'ioi', + 'n_strokes', 'stroke_gap', 'event', + 'stroke_event_ioi'] + assert len(df) == 6 + assert np.allclose(df.t.to_numpy(), gt["starts"]) + assert np.allclose(df.ioi.to_numpy()[:-1], np.diff(gt["starts"])) + assert np.isnan(df.ioi.to_numpy()[-1]) + assert (df.n_strokes == 2).all() + + def test_empty(self): + df = cycle_table([]) + assert len(df) == 0 + assert 'ioi' in df.columns + + +class TestFitAccelerando: + def test_recovers_known_parameters(self): + gt = ro_times(ioi0=2.0, t_double=12.0, n_cycles=20) + t = np.asarray(gt["starts"]) + ioi = np.append(np.diff(t), np.nan) + ioi0, t_double, r2 = fit_accelerando(t, ioi) + # ro_times references its exponential to t=1.0 s, so the fitted + # ioi0 (at t=0) is ioi0 * 2**(1/t_double). + assert t_double == pytest.approx(12.0, rel=0.05) + assert ioi0 == pytest.approx(2.0 * 2 ** (1 / 12.0), rel=0.05) + assert r2 > 0.99 + + def test_steady_pulse_gives_inf(self): + t = np.arange(10, dtype=float) + ioi = np.ones(10) + _, t_double, _ = fit_accelerando(t, ioi) + assert np.isinf(t_double) + + +class TestMotionOnsets: + def test_recovers_ramp_onsets(self): + fs = 50.0 + rng = np.random.default_rng(0) + t = np.arange(int(20 * fs)) / fs + truth = np.array([3.0, 8.0, 13.0, 17.0]) + sig = 0.005 * rng.standard_normal(len(t)) + for on in truth: + sig += 1.0 / (1.0 + np.exp(-(t - on) / 0.1)) * np.exp(-(t - on).clip(0) / 1.5) + onsets = motion_onsets(sig, fs) + matched = [np.min(np.abs(onsets - x)) for x in truth] + assert len(onsets) == len(truth) + assert max(matched) < 0.25 From de9001bbf554d50ada5d688877338c284794b577 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 08:51:36 +0200 Subject: [PATCH 2/6] feat: unified quantity-of-motion signal cores (_qom) One QoM surface unifying the band-limited QoM cores that were duplicated across the stillstanding and Westney-comparisons studies (mocap markers, MediaPipe landmarks, slow postural sway): - band_limited_qom: zero-phase band-pass + speed, with an automatic decimate+SOS regime for very low bands (the "slow sway" case) where a direct high-order band-pass is numerically fragile - accel_to_speed: high-pass, integrate, high-pass accelerometer-to-speed (the stillstanding "corpus method"), with optional gravity normalisation - group_qom / pose_qom: per-body-part QoM over marker/landmark groups - body_scale / normalized_qom: median torso length and framing-invariant, dimensionless QoM (body-lengths/s) - grid_qom: spatial grid QoM heatmap + per-cell series - envelope / bin_series: smoothed z-scored envelopes and binned means Tests use synthetic trajectories with analytically known speeds and band content, including a 2x zoom simulation that must leave normalized_qom invariant, and a decimation-regime check for the low band. Source studies: stillstanding; Westney comparisons (Jensenius). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/__init__.py | 11 ++ musicalgestures/_qom.py | 353 ++++++++++++++++++++++++++++++++++++ tests/test_qom.py | 188 +++++++++++++++++++ 3 files changed, 552 insertions(+) create mode 100644 musicalgestures/_qom.py create mode 100644 tests/test_qom.py diff --git a/musicalgestures/__init__.py b/musicalgestures/__init__.py index 5a8d0d3..e394884 100644 --- a/musicalgestures/__init__.py +++ b/musicalgestures/__init__.py @@ -80,3 +80,14 @@ def __init__(self): fit_accelerando, motion_onsets, ) +from musicalgestures._qom import ( + band_limited_qom, + accel_to_speed, + group_qom, + pose_qom, + body_scale, + normalized_qom, + grid_qom, + envelope, + bin_series, +) diff --git a/musicalgestures/_qom.py b/musicalgestures/_qom.py new file mode 100644 index 0000000..170da7a --- /dev/null +++ b/musicalgestures/_qom.py @@ -0,0 +1,353 @@ +""" +Quantity-of-motion (QoM) signal cores for position, pose and accelerometer +data. + +Band-limited QoM (with an automatic decimate+SOS regime for very low +frequency bands), accelerometer-to-speed integration, per-landmark-group +pose QoM, body-scale normalisation for framing-invariant comparisons, +spatial grid QoM, and small envelope/binning helpers. + +These functions are independent of the MgVideo/MgAudio classes and operate +on plain numpy arrays (marker/landmark trajectories, accelerometer data, +grayscale frame stacks, 1-D signals). + +Sources: stillstanding study and Westney-comparisons study (Jensenius). +""" + +import numpy as np + + +def _interp_nans(x): + """Linearly interpolate non-finite entries of a 1-D array (returns a copy).""" + x = np.asarray(x, float).copy() + m = np.isfinite(x) + if m.sum() > 2 and not m.all(): + x = np.interp(np.arange(len(x)), np.flatnonzero(m), x[m]) + return x + + +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 + + +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) + + +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; empty if the input is too + short) and `fs_out` is its sampling rate (equal to `fs` unless + decimated). + """ + 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])]) + if len(pos) < int(fs) + 5 or not np.isfinite(pos).all(): + return np.array([]), fs + + hi_eff = min(hi, 0.9 * fs / 2) + 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 + 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 + + +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) + + +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 + + +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) + + +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 + + +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. + """ + 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) + return qom / scale, speed / scale, fs_out + + +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. + """ + frames = np.asarray(frames, dtype=np.float32) + 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 diff --git a/tests/test_qom.py b/tests/test_qom.py new file mode 100644 index 0000000..e586a5c --- /dev/null +++ b/tests/test_qom.py @@ -0,0 +1,188 @@ +"""Tests for musicalgestures._qom (band-limited QoM, normalisation, grid QoM). + +All ground truth is synthetic: sinusoidal trajectories with known band +content and analytically known speeds, a simulated 2x camera zoom, and a +moving spot on a synthetic frame stack. +""" +import numpy as np +import pytest + +from musicalgestures import ( + band_limited_qom, + accel_to_speed, + group_qom, + pose_qom, + body_scale, + normalized_qom, + grid_qom, + envelope, + bin_series, +) + + +def circular_motion(fs=100.0, dur=30.0, f=1.0, radius=10.0): + """Circular trajectory: constant speed 2*pi*f*radius, all energy at f Hz.""" + t = np.arange(int(dur * fs)) / fs + return np.column_stack([radius * np.cos(2 * np.pi * f * t), + radius * np.sin(2 * np.pi * f * t)]), t + + +class TestBandLimitedQom: + def test_in_band_speed_close_to_truth(self): + fs = 100.0 + pos, _ = circular_motion(fs, f=1.0, radius=10.0) + speed, fs_out = band_limited_qom(pos, fs, lo=0.3, hi=15.0) + assert fs_out == fs + true_speed = 2 * np.pi * 1.0 * 10.0 + mid = speed[len(speed) // 4: -len(speed) // 4] + assert np.median(mid) == pytest.approx(true_speed, rel=0.05) + + def test_out_of_band_motion_rejected(self): + fs = 100.0 + pos_slow, _ = circular_motion(fs, f=0.05, radius=10.0) # below lo + pos_in, _ = circular_motion(fs, f=1.0, radius=10.0) + s_slow, _ = band_limited_qom(pos_slow, fs, lo=0.3, hi=5.0, + auto_decimate=False) + s_in, _ = band_limited_qom(pos_in, fs, lo=0.3, hi=5.0) + assert s_slow.mean() < 0.05 * s_in.mean() + + def test_low_band_uses_decimation(self): + fs = 100.0 + pos, _ = circular_motion(fs, dur=60.0, f=0.3, radius=10.0) + speed, fs_out = band_limited_qom(pos, fs, lo=0.1, hi=0.5) + assert fs_out < fs # decimated regime + true_speed = 2 * np.pi * 0.3 * 10.0 + mid = speed[len(speed) // 4: -len(speed) // 4] + assert np.median(mid) == pytest.approx(true_speed, rel=0.1) + + def test_nan_interpolation(self): + fs = 100.0 + pos, _ = circular_motion(fs) + pos[200:205] = np.nan + speed, _ = band_limited_qom(pos, fs) + assert len(speed) == len(pos) - 1 + assert np.isfinite(speed).all() + + def test_too_short_input(self): + speed, fs_out = band_limited_qom(np.zeros((10, 2)), 100.0) + assert len(speed) == 0 + + +class TestAccelToSpeed: + def test_recovers_sinusoidal_speed(self): + # Position x(t) = A sin(2 pi f t) => acceleration is its second + # derivative; the integrated speed magnitude should peak near + # A * 2 pi f. + fs = 256.0 + f, A = 1.0, 0.05 + t = np.arange(int(60 * fs)) / fs + w = 2 * np.pi * f + acc = np.column_stack([-A * w ** 2 * np.sin(w * t), + np.zeros_like(t), np.zeros_like(t)]) + speed = accel_to_speed(acc, fs) + mid = speed[len(speed) // 4: -len(speed) // 4] + assert np.max(mid) == pytest.approx(A * w, rel=0.1) + + def test_gravity_normalisation(self): + fs = 256.0 + t = np.arange(int(30 * fs)) / fs + counts_per_g = 340.0 + w = 2 * np.pi * 1.0 + A = 0.05 + # raw counts: gravity on z plus motion on x + acc = np.column_stack([ + -A * w ** 2 * np.sin(w * t) / 9.80665 * counts_per_g, + np.zeros_like(t), + counts_per_g * np.ones_like(t)]) + speed = accel_to_speed(acc, fs, normalize_gravity=True) + mid = speed[len(speed) // 4: -len(speed) // 4] + assert np.max(mid) == pytest.approx(A * w, rel=0.15) + + +class TestGroupAndPoseQom: + def test_group_mean_of_identical_markers(self): + fs = 50.0 + pos, _ = circular_motion(fs, f=0.7) + points = np.stack([pos, pos, pos], axis=1) + qom, speed, fs_out = group_qom(points, fs, lo=0.3, hi=5.0) + single, _ = band_limited_qom(pos, fs, lo=0.3, hi=5.0) + assert qom == pytest.approx(single.mean(), rel=1e-6) + assert np.allclose(speed, single[: len(speed)]) + + def test_pose_qom_single_landmark(self): + fs = 25.0 + pos, _ = circular_motion(fs, f=0.7) + qom, speed, _ = pose_qom(pos, fs) + assert np.isfinite(qom) and qom > 0 + + +class TestBodyScaleNormalisation: + def make_person(self, fs=25.0, dur=30.0, scale=1.0, f=0.6): + """A minimal MediaPipe-indexed skeleton with oscillating wrists.""" + t = np.arange(int(dur * fs)) / fs + L = 33 + lm = np.zeros((len(t), L, 2)) + lm[:, 11] = [100, 100] + lm[:, 12] = [140, 100] + lm[:, 23] = [105, 200] + lm[:, 24] = [135, 200] + wrist = 30 * np.sin(2 * np.pi * f * t) + lm[:, 15, 0] = 80 + wrist + lm[:, 15, 1] = 150 + lm[:, 16, 0] = 160 - wrist + lm[:, 16, 1] = 150 + return lm * scale + + def test_body_scale_torso_length(self): + lm = self.make_person() + # torso: shoulder-mid (120, 100) to hip-mid (120, 200) => 100 px + assert body_scale(lm) == pytest.approx(100.0) + + def test_zoom_invariance(self): + """A simulated 2x camera zoom must not change the normalised QoM.""" + fs = 25.0 + lm1 = self.make_person(fs, scale=1.0) + lm2 = self.make_person(fs, scale=2.0) + q1, _, _ = normalized_qom(lm1, fs) + q2, _, _ = normalized_qom(lm2, fs) + raw1, _, _ = pose_qom(lm1, fs) + raw2, _, _ = pose_qom(lm2, fs) + assert raw2 == pytest.approx(2 * raw1, rel=1e-6) # raw QoM scales + assert q2 == pytest.approx(q1, rel=1e-6) # normalised is invariant + + +class TestGridQom: + def test_moving_spot_lights_up_its_cell(self): + T, H, W = 40, 40, 60 + frames = np.zeros((T, H, W), dtype=np.float32) + # A flickering spot inside cell (row 0, col 0) region + frames[::2, 2:8, 2:8] = 200.0 + series, heat = grid_qom(frames, grid=(6, 4), threshold=8.0) + assert series.shape == (T - 1, 24) + assert heat.shape == (4, 6) + assert np.argmax(heat) == 0 + assert heat[0, 0] > 10 * (heat.sum() - heat[0, 0] + 1e-9) + + def test_region_of_interest(self): + T, H, W = 10, 40, 60 + frames = np.zeros((T, H, W), dtype=np.float32) + frames[::2, :H // 2, :] = 100.0 # motion only in the top half + _, heat_bottom = grid_qom(frames, grid=(2, 2), region=(0, 1, 0.5, 1)) + assert heat_bottom.max() == 0 + + +class TestEnvelopeBinSeries: + def test_envelope_zscored(self): + fs = 25.0 + x = np.sin(2 * np.pi * 0.2 * np.arange(int(20 * fs)) / fs) + e = envelope(x, fs) + assert len(e) == len(x) + assert abs(e.mean()) < 1e-6 + assert e.std() == pytest.approx(1.0, rel=1e-3) + + def test_bin_series_means(self): + fs = 10.0 + x = np.repeat([1.0, 2.0, 3.0], int(fs)) + b = bin_series(x, fs, bin_s=1.0) + assert np.allclose(b, [1.0, 2.0, 3.0]) + assert len(bin_series(np.ones(5), fs)) == 0 From ef26c9ce09f3b1f5414f082358689cf1ad588fc2 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 08:52:02 +0200 Subject: [PATCH 3/6] feat: audio-motion alignment tools (_alignment) Stream-to-stream alignment methods for signals and event lists: - xcorr_lag: canonical vectorized lead/lag estimate (scipy.signal.correlate/correlation_lags), normalized, preferring the smallest-|lag| candidate among near-ties so periodic envelopes do not alias; envelope_lag is kept as a thin wrapper preserving the ro study's interface, unifying the two studies' implementations - per_cycle_motion_delta: nearest-to-start motion-onset assignment with descending-order dedup over overlapping lookback windows (ported faithfully from the ro study, semantics documented in the docstring) - anchor_and_match + offset_stats: per-take relative event alignment for streams with independent clocks (anchor the strongest events at t=0, match the rest within a window; no absolute sync claim); the 0.15 s window is a provisional default reimplemented from the cymbal paper - sliding_correlation: windowed local coupling profile - envelope_agreement: N-source envelope correlation matrix (cross-view agreement), built on _qom.envelope Tests: known 0.3 s envelope shift recovery, periodic smallest-lag preference, synthetic offset events for anchor_and_match, and a built-in coupling change profiled by sliding_correlation. Source studies: ro; cymbal comparison; Westney comparisons (Jensenius). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/__init__.py | 9 + musicalgestures/_alignment.py | 313 ++++++++++++++++++++++++++++++++++ tests/test_alignment.py | 165 ++++++++++++++++++ 3 files changed, 487 insertions(+) create mode 100644 musicalgestures/_alignment.py create mode 100644 tests/test_alignment.py diff --git a/musicalgestures/__init__.py b/musicalgestures/__init__.py index e394884..bd90a2a 100644 --- a/musicalgestures/__init__.py +++ b/musicalgestures/__init__.py @@ -91,3 +91,12 @@ def __init__(self): envelope, bin_series, ) +from musicalgestures._alignment import ( + xcorr_lag, + envelope_lag, + per_cycle_motion_delta, + anchor_and_match, + offset_stats, + sliding_correlation, + envelope_agreement, +) diff --git a/musicalgestures/_alignment.py b/musicalgestures/_alignment.py new file mode 100644 index 0000000..d8b8cc0 --- /dev/null +++ b/musicalgestures/_alignment.py @@ -0,0 +1,313 @@ +""" +Audio--motion (and stream-to-stream) alignment tools. + +Lead/lag estimation by cross-correlation, per-cycle motion-onset deltas, +anchor-and-match relative event alignment with offset statistics, +sliding-window coupling, and N-source envelope agreement. + +These functions are independent of the MgVideo/MgAudio classes and operate +on plain numpy arrays (envelopes, onset/event time lists). + +Sources: ro study, cymbal-comparison study and Westney-comparisons study +(Jensenius). +""" + +import numpy as np + +from musicalgestures._qom import envelope + + +def xcorr_lag(x, y, fs, max_lag=1.5): + """ + Canonical lead/lag estimate between two signals by vectorized + cross-correlation: the lag of `y` relative to `x` that maximizes their + correlation, searched within +/- `max_lag` seconds. Positive lag means + `y` happens after `x`. + + Both signals are mean-removed and the correlation is normalized to a + Pearson-like coefficient over the full window. Among near-tied maxima + (common for periodic envelopes, where peaks recur at +/- one period), + the smallest-magnitude lag is returned rather than an arbitrary aliased + one. + + Source: Westney-comparisons study (Jensenius) -- camera/audio residual + sync by onset-envelope cross-correlation; unified with the ro study's + `envelope_lag`. + + Args: + x (np.ndarray): Reference 1-D signal. + y (np.ndarray): Comparison 1-D signal (same sampling rate). + fs (float): Sampling rate of both signals (Hz). + max_lag (float, optional): Maximum absolute lag to search (s). Defaults to 1.5. + + Returns: + tuple: `(lag, r)` where `lag` is the lag of `y` relative to `x` in + seconds (positive = `y` later) and `r` is the normalized correlation + at that lag. + """ + from scipy.signal import correlate, correlation_lags + x = np.asarray(x, float) + y = np.asarray(y, float) + n = min(len(x), len(y)) + x, y = x[:n] - x[:n].mean(), y[:n] - y[:n].mean() + cc = correlate(y, x, mode="full") + lags = correlation_lags(n, n, mode="full") + m = np.abs(lags) <= max(1, int(round(max_lag * fs))) + cc, lags = cc[m], lags[m] + cc = cc / (n * x.std() * y.std() + 1e-12) + # Near-exact ties are floating-point noise; prefer the smallest |lag|. + cand = np.flatnonzero(cc >= cc.max() - 1e-9) + best = cand[np.argmin(np.abs(lags[cand]))] + return lags[best] / fs, float(cc[best]) + + +def envelope_lag(x, y, rate, max_lag_s=1.5): + """ + Lag (s) of `y` relative to `x` maximizing their correlation. Positive + lag = `y` happens after `x`. Thin wrapper around `xcorr_lag`, kept as + the ro study's interface for envelope-to-envelope lags (e.g. voice + envelope vs motion envelope). + + Source: ro study (Jensenius). + + Args: + x (np.ndarray): Reference envelope. + y (np.ndarray): Comparison envelope (same sampling rate). + rate (float): Sampling rate of both envelopes (Hz). + max_lag_s (float, optional): Maximum absolute lag to search (s). Defaults to 1.5. + + Returns: + tuple: `(lag, r)` as in `xcorr_lag`. + """ + return xcorr_lag(x, y, rate, max_lag=max_lag_s) + + +def per_cycle_motion_delta(cycle_starts, motion_times, lookback=0.3): + """ + For each cycle, the time of its assigned motion onset minus the + cycle start (NaN if none). Each cycle's natural window is + [start-lookback, next_start): the lookback lets a motion onset that + anticipates this cycle's stroke be credited to it, while the + unshifted upper bound at next_start lets genuine post-onset motion + (motion after the cycle's own stroke) be credited too. + + Because windows overlap by up to `lookback`, a single motion onset + can fall in two consecutive cycles' windows. To avoid double- + counting (and to stop a slow cycle from "stealing" a fast cycle's + only motion onset when IOI < lookback), cycles are processed in + DESCENDING time order and each claims the unclaimed motion onset + still inside its window that is NEAREST to its own cycle start + (ties broken toward the earlier onset), marking it claimed. + A cyclic gesture can plausibly produce two motion peaks per slow + cycle: a forward-stroke peak close to the sound onset and a later + return-stroke peak near the tail of the cycle's window. The + forward-stroke response to the sound is the one of interest, so + nearest-to-start assignment credits that one to the cycle instead + of the return-stroke peak that a latest-wins rule would grab. + Processing cycles in descending order still ensures a fast climax + cycle can secure its own motion onset before an earlier, slower + cycle's wider (lookback-only) window claims it instead. + + Source: ro study (Jensenius) -- rowing-gesture onsets vs drum cycles. + + Args: + cycle_starts (np.ndarray): Cycle start times (s), ascending. + motion_times (np.ndarray): Motion onset times (s), any order. + lookback (float, optional): How far before its start a cycle may claim + a motion onset (s). Defaults to 0.3. + + Returns: + np.ndarray: One delta (s) per cycle: assigned motion onset time minus + cycle start, NaN where no onset was assigned. + """ + cycle_starts = np.asarray(cycle_starts, float) + motion_times = np.sort(np.asarray(motion_times, float)) + ends = np.append(cycle_starts[1:], np.inf) + claimed = np.zeros(len(motion_times), dtype=bool) + out = np.full(len(cycle_starts), np.nan) + for i in range(len(cycle_starts) - 1, -1, -1): + s, e = cycle_starts[i], ends[i] + in_window = (motion_times >= s - lookback) & (motion_times < e) & ~claimed + idx = np.flatnonzero(in_window) + if len(idx): + j = idx[np.argmin(np.abs(motion_times[idx] - s))] + claimed[j] = True + out[i] = motion_times[j] - s + return out + + +def anchor_and_match(times_a, times_b, anchor_a=None, anchor_b=None, + weights_a=None, weights_b=None, window=0.15): + """ + Per-take relative event alignment between two streams with independent + clocks: both streams are shifted so that their anchor events -- by + default the strongest event of each stream, physically the same moment + (e.g. the hardest strike) -- sit at t = 0; then every non-anchor event + of stream `a` is matched to the nearest event of stream `b` within + +/- `window` seconds, and the signed offset (`b` minus `a`; positive = + `b` later) is recorded. The anchor pair contributes an offset of 0 by + construction and is excluded. No absolute cross-stream synchronisation + is claimed: the result measures whether the two streams agree on the + RELATIVE timing of the remaining events. + + The default matching window of 0.15 s is a PROVISIONAL default, + reimplemented from the cymbal-comparison paper's method description. + + Source: cymbal-comparison study (Jensenius) -- audio-onset vs kinematic- + impact (and video/pose) timing agreement without a common clock. + + Args: + times_a (np.ndarray): Event times (s) of stream a (e.g. kinematic impacts). + times_b (np.ndarray): Event times (s) of stream b (e.g. audio onsets). + anchor_a (float, optional): Anchor time in stream a. Defaults to None. + anchor_b (float, optional): Anchor time in stream b. Defaults to None. + weights_a (np.ndarray, optional): Event strengths for stream a, used to + pick the anchor (argmax) when `anchor_a` is None. Defaults to None. + weights_b (np.ndarray, optional): Event strengths for stream b, used to + pick the anchor (argmax) when `anchor_b` is None. Defaults to None. + window (float, optional): Maximum absolute offset (s) for a match. + Defaults to 0.15. + + Returns: + np.ndarray: Signed offsets (s), one per matched non-anchor event of + stream a (`b` minus `a`). + + Raises: + ValueError: If an anchor can be determined for neither stream (no + anchor time and no weights given). + """ + times_a = np.asarray(times_a, float) + times_b = np.asarray(times_b, float) + + def pick_anchor(times, anchor, weights, name): + if anchor is not None: + return float(anchor) + if weights is not None: + return float(times[int(np.argmax(weights))]) + raise ValueError( + f"anchor_and_match: provide anchor_{name} or weights_{name} " + f"to determine the anchor event of stream {name}") + + a0 = pick_anchor(times_a, anchor_a, weights_a, "a") + b0 = pick_anchor(times_b, anchor_b, weights_b, "b") + a = times_a - a0 + b = times_b - b0 + a = a[np.abs(a) > 1e-12] # exclude the anchor itself + if len(a) == 0 or len(b) == 0: + return np.array([]) + offsets = [] + for t in a: + j = int(np.argmin(np.abs(b - t))) + if abs(b[j] - t) <= window: + offsets.append(b[j] - t) + return np.asarray(offsets, float) + + +def offset_stats(offsets): + """ + Summary statistics of a signed-offset distribution (s), as reported for + anchor-and-match timing agreement. + + Source: cymbal-comparison study (Jensenius). + + Args: + offsets (np.ndarray): Signed offsets (s), e.g. from `anchor_and_match`. + + Returns: + dict: `n`, `median`, `mean`, `std`, `iqr` (interquartile range), + `abs_median` (median absolute offset), `min` and `max`; statistics + are NaN when `n` is 0. + """ + offsets = np.asarray(offsets, float) + offsets = offsets[np.isfinite(offsets)] + if len(offsets) == 0: + nan = float("nan") + return dict(n=0, median=nan, mean=nan, std=nan, iqr=nan, + abs_median=nan, min=nan, max=nan) + q1, q3 = np.percentile(offsets, [25, 75]) + return dict( + n=int(len(offsets)), + median=float(np.median(offsets)), + mean=float(np.mean(offsets)), + std=float(np.std(offsets)), + iqr=float(q3 - q1), + abs_median=float(np.median(np.abs(offsets))), + min=float(np.min(offsets)), + max=float(np.max(offsets)), + ) + + +def sliding_correlation(x, y, fs, window=10.0, step=2.0, min_std=1e-6): + """ + Local (windowed) Pearson correlation profile between two signals: + correlation within sliding windows of `window` seconds, hopped by + `step` seconds. Windows where either signal is (near-)constant yield + NaN. Useful for testing whether two streams that are uncorrelated + globally couple locally (e.g. motion effort vs loudness). + + Source: Westney-comparisons study (Jensenius) -- local motion-loudness + coupling in 10 s windows. + + Args: + x (np.ndarray): First 1-D signal. + y (np.ndarray): Second 1-D signal (same sampling rate). + fs (float): Sampling rate of both signals (Hz). + window (float, optional): Window length (s). Defaults to 10.0. + step (float, optional): Hop between windows (s). Defaults to 2.0. + min_std (float, optional): Minimum in-window standard deviation for a + valid correlation. Defaults to 1e-6. + + Returns: + tuple: `(times, r)` where `times` are the window centres (s) and `r` + the windowed correlations (NaN where undefined). + """ + x = np.asarray(x, float) + y = np.asarray(y, float) + n = min(len(x), len(y)) + W = max(2, int(round(window * fs))) + S = max(1, int(round(step * fs))) + times, r = [], [] + for s0 in range(0, n - W + 1, S): + a, b = x[s0:s0 + W], y[s0:s0 + W] + times.append((s0 + W / 2) / fs) + if a.std() > min_std and b.std() > min_std: + r.append(float(np.corrcoef(a, b)[0, 1])) + else: + r.append(np.nan) + return np.asarray(times), np.asarray(r) + + +def envelope_agreement(signals, fs, smooth=1.0): + """ + Agreement among N parallel envelopes (e.g. per-camera-angle quantity-of- + motion curves of the same performance): the signals are truncated to + their common length, smoothed and z-scored (see + `musicalgestures.envelope`), and their Pearson correlation matrix is + computed. The mean off-diagonal correlation summarises how well the + sources agree. + + Source: Westney-comparisons study (Jensenius) -- cross-view agreement of + five camera angles' motion envelopes. + + Args: + signals (sequence): Sequence of N 1-D arrays (possibly of different + lengths), or a 2-D array of shape (N, T). + fs (float): Sampling rate of the signals (Hz). + smooth (float, optional): Envelope smoothing window (s); None or 0 + disables smoothing. Defaults to 1.0. + + Returns: + tuple: `(C, mean_r)` where `C` is the (N, N) correlation matrix and + `mean_r` is the mean of its upper off-diagonal entries (NaN for + fewer than two signals). + """ + signals = [np.asarray(s, float) for s in signals] + if len(signals) == 0: + return np.zeros((0, 0)), float("nan") + L = min(len(s) for s in signals) + M = np.array([envelope(s[:L], fs, smooth=smooth, normalize=True) + for s in signals]) + C = np.corrcoef(M) if len(signals) > 1 else np.ones((1, 1)) + iu = np.triu_indices(len(signals), k=1) + mean_r = float(np.mean(C[iu])) if len(iu[0]) else float("nan") + return C, mean_r diff --git a/tests/test_alignment.py b/tests/test_alignment.py new file mode 100644 index 0000000..d0cfd6e --- /dev/null +++ b/tests/test_alignment.py @@ -0,0 +1,165 @@ +"""Tests for musicalgestures._alignment (lead/lag, matching, coupling profiles). + +All ground truth is synthetic: lagged envelopes with known shifts, event +lists with known offsets, and signals with a built-in coupling change. +""" +import numpy as np +import pytest + +from musicalgestures import ( + xcorr_lag, + envelope_lag, + per_cycle_motion_delta, + anchor_and_match, + offset_stats, + sliding_correlation, + envelope_agreement, +) + + +def noisy_envelope(fs=50.0, dur=30.0, seed=1): + rng = np.random.default_rng(seed) + t = np.arange(int(dur * fs)) / fs + x = np.zeros_like(t) + for a in (1.0, 0.7, 0.4): + f = rng.uniform(0.2, 2.0) + x += a * np.sin(2 * np.pi * f * t + rng.uniform(0, 2 * np.pi)) + return x + 0.05 * rng.standard_normal(len(t)) + + +class TestXcorrLag: + def test_recovers_300ms_shift(self): + fs = 50.0 + x = noisy_envelope(fs) + shift = int(0.3 * fs) + y = np.roll(x, shift) # y happens 0.3 s after x + lag, r = xcorr_lag(x[shift:-shift], y[shift:-shift], fs, max_lag=1.5) + assert lag == pytest.approx(0.3, abs=1.5 / fs) + assert r > 0.9 + + def test_negative_lag(self): + fs = 50.0 + x = noisy_envelope(fs, seed=2) + shift = int(0.4 * fs) + y = np.roll(x, -shift) # y precedes x + lag, _ = xcorr_lag(x[shift:-shift], y[shift:-shift], fs, max_lag=1.5) + assert lag == pytest.approx(-0.4, abs=1.5 / fs) + + def test_periodic_prefers_smallest_lag(self): + # A pure sinusoid has equally good peaks at lag +/- one period; + # the smallest-|lag| candidate (0) must win. + fs = 100.0 + t = np.arange(int(20 * fs)) / fs + x = np.sin(2 * np.pi * 1.0 * t) + lag, r = xcorr_lag(x, x, fs, max_lag=1.5) + assert lag == 0.0 + assert r > 0.99 + + def test_envelope_lag_wrapper_agrees(self): + fs = 50.0 + x = noisy_envelope(fs, seed=3) + y = np.roll(x, int(0.2 * fs)) + assert envelope_lag(x, y, fs) == xcorr_lag(x, y, fs) + + +class TestPerCycleMotionDelta: + def test_nearest_to_start_and_dedup(self): + starts = np.array([0.0, 2.0, 4.0, 5.0]) + # one anticipating onset, one late return-stroke peak in cycle 0, + # one onset shared between cycle 2's window and cycle 3's lookback + motion = np.array([-0.1, 1.5, 2.05, 4.9]) + out = per_cycle_motion_delta(starts, motion, lookback=0.3) + assert out[0] == pytest.approx(-0.1) + assert out[1] == pytest.approx(0.05) + # 4.9 is in both cycle 2's window [3.7, 5.0) and cycle 3's [4.7, inf); + # descending order lets cycle 3 claim it first, leaving cycle 2 empty. + assert np.isnan(out[2]) + assert out[3] == pytest.approx(-0.1) + + def test_no_motion(self): + out = per_cycle_motion_delta([0.0, 1.0], []) + assert np.all(np.isnan(out)) + + +class TestAnchorAndMatch: + def test_recovers_known_offsets(self): + # Stream b = stream a + per-event offsets, on a shifted clock. + a = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + true_off = np.array([0.0, 0.05, -0.04, 0.02, 0.10]) + b = a + true_off + 7.3 # independent clock + wa = np.array([1, 5, 1, 1, 1]) # anchor: a[1] / b[1] + wb = wa.copy() + off = anchor_and_match(a, b, weights_a=wa, weights_b=wb, window=0.15) + # relative to the anchor pair (offset 0.05), excluded by construction + expected = np.array([0.0, -0.04, 0.02, 0.10]) - 0.05 + assert np.allclose(np.sort(off), np.sort(expected), atol=1e-9) + + def test_window_excludes_far_events(self): + a = np.array([0.0, 1.0, 2.0]) + b = np.array([0.0, 1.5]) # the 1.0 event has no match within 0.15 + off = anchor_and_match(a, b, anchor_a=0.0, anchor_b=0.0, window=0.15) + assert len(off) == 0 + + def test_requires_anchor_or_weights(self): + with pytest.raises(ValueError): + anchor_and_match([0.0, 1.0], [0.0, 1.0]) + + def test_offset_stats(self): + off = np.array([-0.02, 0.0, 0.02, 0.04]) + st = offset_stats(off) + assert st["n"] == 4 + assert st["median"] == pytest.approx(0.01) + assert st["mean"] == pytest.approx(0.01) + assert st["min"] == pytest.approx(-0.02) + assert st["max"] == pytest.approx(0.04) + assert offset_stats([])["n"] == 0 + assert np.isnan(offset_stats([])["median"]) + + +class TestSlidingCorrelation: + def test_profiles_coupling_change(self): + # First half: y coupled to x; second half: y independent. + fs = 25.0 + dur = 60.0 + rng = np.random.default_rng(4) + x = noisy_envelope(fs, dur, seed=5) + y = x.copy() + half = len(x) // 2 + y[half:] = noisy_envelope(fs, dur, seed=6)[half:] + times, r = sliding_correlation(x, y, fs, window=10.0, step=2.0) + first = r[times < dur / 2 - 5] + second = r[times > dur / 2 + 5] + assert np.nanmedian(first) > 0.95 + assert np.nanmedian(second) < 0.5 + + def test_constant_window_is_nan(self): + fs = 10.0 + x = np.ones(int(30 * fs)) + y = noisy_envelope(fs, 30.0) + _, r = sliding_correlation(x, y, fs, window=5.0, step=5.0) + assert np.all(np.isnan(r)) + + +class TestEnvelopeAgreement: + def test_agreeing_sources(self): + fs = 25.0 + base = noisy_envelope(fs, seed=7) + rng = np.random.default_rng(8) + views = [base + 0.1 * rng.standard_normal(len(base)) for _ in range(3)] + C, mean_r = envelope_agreement(views, fs) + assert C.shape == (3, 3) + assert np.allclose(np.diag(C), 1.0) + assert mean_r > 0.9 + + def test_disagreeing_source(self): + fs = 25.0 + a = noisy_envelope(fs, seed=9) + b = noisy_envelope(fs, seed=10) + C, mean_r = envelope_agreement([a, b], fs) + assert abs(mean_r) < 0.5 + + def test_unequal_lengths_truncated(self): + fs = 25.0 + a = noisy_envelope(fs, dur=30.0, seed=11) + C, _ = envelope_agreement([a, a[: len(a) - 40]], fs) + assert C.shape == (2, 2) From b676b5ef23316c9c05aeb98e29cc35c502d275da Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 08:52:33 +0200 Subject: [PATCH 4/6] feat: scipy-only audio feature extraction (_audiofeatures) Plain-numpy audio descriptors complementing the librosa-based, figure-producing MgAudio methods with lightweight numeric outputs (a new module rather than MgAudio methods, following the grain of the _analysis.py plain-function convention; MgAudio methods return MgFigures, these return arrays/floats): - rms_envelope: windowed RMS energy envelope - spectral_flux / spectral_flux_onsets: positive-STFT-difference onset detection function, peak-picked with the canonical pick_peaks at an adaptive mean+std threshold - energy_onsets: RMS-envelope half-wave-rectified-difference onsets (0.15 x peak / 0.06 s provisional defaults from the cymbal paper) - t60_backward_decay: ISO-3382-style backward-decay reverberation time (-5..-35 dB span, -5..-25 fallback, 6 dB re-rise stop; provisional defaults reimplemented from the cymbal paper's prose) - attack_spectral_centroid: energy-weighted STFT centroid over the first 120 ms after the envelope peak Tests: click trains with known onset times, synthetic exponential decays with exact T60 (including the T20-fallback path), and bright-vs-dark attack ordering. Source studies: cymbal comparison; Westney comparisons (Jensenius). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/__init__.py | 8 + musicalgestures/_audiofeatures.py | 246 ++++++++++++++++++++++++++++++ tests/test_audiofeatures.py | 99 ++++++++++++ 3 files changed, 353 insertions(+) create mode 100644 musicalgestures/_audiofeatures.py create mode 100644 tests/test_audiofeatures.py diff --git a/musicalgestures/__init__.py b/musicalgestures/__init__.py index bd90a2a..c7d5d10 100644 --- a/musicalgestures/__init__.py +++ b/musicalgestures/__init__.py @@ -100,3 +100,11 @@ def __init__(self): sliding_correlation, envelope_agreement, ) +from musicalgestures._audiofeatures import ( + rms_envelope, + spectral_flux, + spectral_flux_onsets, + energy_onsets, + t60_backward_decay, + attack_spectral_centroid, +) diff --git a/musicalgestures/_audiofeatures.py b/musicalgestures/_audiofeatures.py new file mode 100644 index 0000000..313478f --- /dev/null +++ b/musicalgestures/_audiofeatures.py @@ -0,0 +1,246 @@ +""" +Scipy-only audio feature extraction for sound--motion analysis. + +RMS envelopes, spectral-flux and energy-based onset detection, T60-style +backward-decay reverberation time, and attack spectral centroid. + +These functions are independent of the MgAudio class and operate on plain +numpy waveforms (`y`, `sr`), so they can be used on any mono audio array. +They complement the librosa-based, figure-producing methods of `MgAudio` +with lightweight numeric outputs. All onset detectors use the canonical +peak-picker (`musicalgestures.pick_peaks`). + +Sources: cymbal-comparison study and Westney-comparisons study (Jensenius). +""" + +import numpy as np + +from musicalgestures._peaks import pick_peaks + + +def rms_envelope(y, sr, window=0.02): + """ + RMS energy envelope over consecutive, non-overlapping windows. + + Source: cymbal-comparison study (Jensenius); also the + Westney-comparisons study's high-rate sync envelope. + + Args: + y (np.ndarray): Mono waveform. + sr (int): Sampling rate (Hz). + window (float, optional): Window length in seconds. Defaults to 0.02. + + Returns: + tuple: `(env, rate)` where `env` is the RMS envelope (one value per + window) and `rate` its sampling rate (1 / `window` Hz). + """ + y = np.asarray(y, float) + win = max(1, int(round(window * sr))) + n = len(y) // win + env = np.sqrt((y[:n * win].reshape(n, win) ** 2).mean(axis=1) + 1e-12) + return env, sr / win + + +def spectral_flux(y, sr, nperseg=2048, noverlap=1536): + """ + Spectral-flux onset-detection function: the positive first difference of + the STFT magnitude, summed over frequency and normalized to a maximum of + 1. Rises sharply at note/percussion onsets. + + Source: Westney-comparisons study (Jensenius). + + Args: + y (np.ndarray): Mono waveform. + sr (int): Sampling rate (Hz). + nperseg (int, optional): STFT window length in samples. Defaults to 2048. + noverlap (int, optional): STFT window overlap in samples. Defaults to 1536. + + Returns: + tuple: `(flux, times)` where `flux` is the onset-detection function and + `times` are its frame times in seconds (frame rate + `sr / (nperseg - noverlap)`). + """ + from scipy.signal import stft + y = np.asarray(y, float) + f, t, Z = stft(y, sr, nperseg=nperseg, noverlap=noverlap) + mag = np.abs(Z) + flux = np.maximum(0, np.diff(mag, axis=1)).sum(axis=0) + flux = flux / (flux.max() + 1e-9) + return flux, t[1:] + + +def spectral_flux_onsets(y, sr, nperseg=2048, noverlap=1536, + threshold=None, min_interval=0.05): + """ + Onset times from the spectral-flux onset-detection function (see + `spectral_flux`), peak-picked with the canonical peak-picker. By default + the threshold adapts to the signal as mean + 1 standard deviation of the + flux, following the source study. + + Source: Westney-comparisons study (Jensenius). + + Args: + y (np.ndarray): Mono waveform. + sr (int): Sampling rate (Hz). + nperseg (int, optional): STFT window length in samples. Defaults to 2048. + noverlap (int, optional): STFT window overlap in samples. Defaults to 1536. + threshold (float, optional): Absolute flux threshold (the flux has + maximum 1). Defaults to None, which uses mean + std of the flux. + min_interval (float, optional): Minimum inter-onset interval (s). + Defaults to 0.05. + + Returns: + np.ndarray: Onset times in seconds. + """ + flux, times = spectral_flux(y, sr, nperseg=nperseg, noverlap=noverlap) + if len(flux) == 0: + return np.array([]) + if threshold is None: + threshold = float(flux.mean() + flux.std()) + rate = sr / (nperseg - noverlap) + idx = pick_peaks(flux, fs=rate, smooth=None, rel_threshold=None, + threshold=threshold, min_interval=min_interval, + rel_prominence=None) + return times[idx] + + +def energy_onsets(y, sr, window=0.02, rel_threshold=0.15, min_interval=0.06): + """ + Energy-based onset times: the RMS envelope's half-wave-rectified first + difference is peak-picked (canonical peak-picker) at a threshold + relative to the per-file peak, with a minimum inter-onset interval. + + Reliable for discrete strokes; over-fragments sustained rolls/tremolo + and can trigger on near-noise material. The default constants + (0.15 x peak, 0.06 s) are PROVISIONAL defaults reimplemented from the + cymbal-comparison paper's method description. + + Source: cymbal-comparison study (Jensenius). + + Args: + y (np.ndarray): Mono waveform. + sr (int): Sampling rate (Hz). + window (float, optional): RMS window length (s). Defaults to 0.02. + rel_threshold (float, optional): Threshold as a fraction of the + onset-function's peak. Defaults to 0.15. + min_interval (float, optional): Minimum inter-onset interval (s). + Defaults to 0.06. + + Returns: + np.ndarray: Onset times in seconds. + """ + env, rate = rms_envelope(y, sr, window=window) + odf = np.maximum(np.diff(env), 0) + idx = pick_peaks(odf, fs=rate, smooth=None, rel_threshold=rel_threshold, + min_interval=min_interval, rel_prominence=None) + return (idx + 1) / rate + + +def t60_backward_decay(y, sr, window=0.02, spans=((-5, -35), (-5, -25)), + rerise_db=6.0): + """ + T60-style reverberation time by an ISO-3382-style backward-decay method: + the RMS envelope is converted to dB relative to its peak; from the peak + the decay is followed (stopping if the envelope re-rises by more than + `rerise_db` dB, marking a new onset), and T60 is estimated by linear + regression of the dB curve over the first available level span -- + by default -5 to -35 dB (a T30 measure, extrapolated x2), falling back + to -5 to -25 dB (T20, x3) when the deeper level is not reached. + + The constants (20 ms window, -5/-35 with -5/-25 fallback, 6 dB re-rise + stop) are PROVISIONAL defaults reimplemented from the cymbal-comparison + paper's method description. + + Source: cymbal-comparison study (Jensenius) -- instrument decay of + damped vs undamped cymbal strokes. + + Args: + y (np.ndarray): Mono waveform of a decaying event (analysis starts at + the envelope peak). + sr (int): Sampling rate (Hz). + window (float, optional): RMS window length (s). Defaults to 0.02. + spans (tuple, optional): Level spans (dB re. peak) to try in order, + each a `(top, bottom)` pair. Defaults to ((-5, -35), (-5, -25)). + rerise_db (float, optional): Stop following the decay when the envelope + re-rises this many dB above its running minimum. Defaults to 6.0. + + Returns: + tuple: `(t60, span)` where `t60` is the estimated reverberation time in + seconds (NaN if no span was usable) and `span` is the `(top, bottom)` + pair actually used (None if none). + """ + env, rate = rms_envelope(y, sr, window=window) + if len(env) < 3: + return float("nan"), None + db = 20 * np.log10(env / (env.max() + 1e-12) + 1e-12) + p = int(np.argmax(env)) + + # Follow the decay from the peak until the envelope re-rises. + stop = len(db) + runmin = db[p] + for i in range(p + 1, len(db)): + runmin = min(runmin, db[i]) + if db[i] > runmin + rerise_db: + stop = i + break + seg = db[p:stop] + t = np.arange(len(seg)) / rate + + for top, bottom in spans: + above_top = np.flatnonzero(seg <= top) + below_bot = np.flatnonzero(seg <= bottom) + if not len(above_top) or not len(below_bot): + continue + i1, i2 = above_top[0], below_bot[0] + if i2 <= i1 + 1: + continue + tt, ss = t[i1:i2 + 1], seg[i1:i2 + 1] + slope, _ = np.polyfit(tt, ss, 1) + if slope >= 0: + continue + return float(-60.0 / slope), (top, bottom) + return float("nan"), None + + +def attack_spectral_centroid(y, sr, attack=0.12, nperseg=2048, hop=512, + window=0.02): + """ + Attack spectral centroid: the energy-weighted mean STFT spectral + centroid over the first `attack` seconds after the RMS-envelope peak. + Discriminates, e.g., strike placement and damping on a cymbal. + + The constants (120 ms attack, 2048-sample Hann window, 512 hop) are + PROVISIONAL defaults reimplemented from the cymbal-comparison paper's + method description. + + Source: cymbal-comparison study (Jensenius). + + Args: + y (np.ndarray): Mono waveform. + sr (int): Sampling rate (Hz). + attack (float, optional): Analysis span after the envelope peak (s). + Defaults to 0.12. + nperseg (int, optional): STFT window length in samples. Defaults to 2048. + hop (int, optional): STFT hop in samples. Defaults to 512. + window (float, optional): RMS window length (s) used to locate the + envelope peak. Defaults to 0.02. + + Returns: + float: The attack spectral centroid in Hz (NaN if the attack segment + is too short to analyse). + """ + from scipy.signal import stft + y = np.asarray(y, float) + env, rate = rms_envelope(y, sr, window=window) + if len(env) == 0: + return float("nan") + i0 = int(np.argmax(env) / rate * sr) + seg = y[i0:i0 + int(attack * sr)] + if len(seg) < nperseg // 4: + return float("nan") + f, t, Z = stft(seg, sr, nperseg=min(nperseg, len(seg)), + noverlap=min(nperseg, len(seg)) - min(hop, len(seg) // 2)) + power = np.abs(Z) ** 2 + frame_energy = power.sum(axis=0) + cen = (power * f[:, None]).sum(axis=0) / (frame_energy + 1e-12) + return float((cen * frame_energy).sum() / (frame_energy.sum() + 1e-12)) diff --git a/tests/test_audiofeatures.py b/tests/test_audiofeatures.py new file mode 100644 index 0000000..739f34a --- /dev/null +++ b/tests/test_audiofeatures.py @@ -0,0 +1,99 @@ +"""Tests for musicalgestures._audiofeatures (onsets, T60, attack centroid). + +All ground truth is synthetic (see tests/_synth.py): click trains with known +onset times and tones with an exact exponential decay. +""" +import numpy as np +import pytest + +from musicalgestures import ( + rms_envelope, + spectral_flux, + spectral_flux_onsets, + energy_onsets, + t60_backward_decay, + attack_spectral_centroid, +) +from _synth import SR, click_train, decaying_tone + + +class TestRmsEnvelope: + def test_tracks_amplitude(self): + sr = SR + t = np.arange(sr) / sr + y = np.sin(2 * np.pi * 440 * t) + y[: sr // 2] *= 0.1 + env, rate = rms_envelope(y, sr, window=0.02) + assert rate == pytest.approx(50.0) + assert env[len(env) // 4] == pytest.approx(0.1 / np.sqrt(2), rel=0.05) + assert env[3 * len(env) // 4] == pytest.approx(1 / np.sqrt(2), rel=0.05) + + +class TestSpectralFluxOnsets: + def test_click_train(self): + truth = np.array([0.5, 1.0, 1.6, 2.3, 3.1]) + y = click_train(truth) + onsets = spectral_flux_onsets(y, SR) + assert len(onsets) == len(truth) + # STFT frame resolution: hop = 512 samples ~ 23 ms + assert np.max(np.abs(np.sort(onsets) - truth)) < 0.05 + + def test_flux_normalized(self): + y = click_train([0.2, 0.6]) + flux, times = spectral_flux(y, SR) + assert flux.max() == pytest.approx(1.0) + assert len(flux) == len(times) + + def test_silence_has_no_onsets(self): + y = np.zeros(SR) + assert len(spectral_flux_onsets(y, SR)) == 0 + + +class TestEnergyOnsets: + def test_click_train(self): + truth = np.array([0.4, 1.1, 1.9, 2.4]) + y = click_train(truth) + onsets = energy_onsets(y, SR) + assert len(onsets) == len(truth) + assert np.max(np.abs(np.sort(onsets) - truth)) < 0.05 + + def test_min_interval_merges_flams(self): + y = click_train([1.0, 1.03, 2.0]) + onsets = energy_onsets(y, SR, min_interval=0.06) + assert len(onsets) == 2 + + +class TestT60: + @pytest.mark.parametrize("true_t60", [0.8, 2.0, 6.0]) + def test_recovers_known_decay(self, true_t60): + y = decaying_tone(true_t60, dur=0.1 + 0.7 * true_t60) + t60, span = t60_backward_decay(y, SR) + assert span == (-5, -35) + assert t60 == pytest.approx(true_t60, rel=0.1) + + def test_t20_fallback_on_short_decay(self): + # Truncate before the -35 dB level is reached + true_t60 = 4.0 + y = decaying_tone(true_t60, dur=0.05 + true_t60 * 30 / 60 * 0.9) + t60, span = t60_backward_decay(y, SR) + assert span == (-5, -25) + assert t60 == pytest.approx(true_t60, rel=0.15) + + def test_no_decay_gives_nan(self): + y = np.sin(2 * np.pi * 440 * np.arange(SR) / SR) # steady tone + t60, span = t60_backward_decay(y, SR) + assert np.isnan(t60) and span is None + + +class TestAttackSpectralCentroid: + def test_orders_bright_vs_dark(self): + sr = SR + t = np.arange(int(0.5 * sr)) / sr + env = np.exp(-t * 8) + dark = np.sin(2 * np.pi * 300 * t) * env + bright = np.sin(2 * np.pi * 3000 * t) * env + c_dark = attack_spectral_centroid(dark, sr) + c_bright = attack_spectral_centroid(bright, sr) + assert c_dark == pytest.approx(300, rel=0.2) + assert c_bright == pytest.approx(3000, rel=0.2) + assert c_bright > c_dark From 7537f5b6bab3328c933ec91953f2e83e46dec396 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 08:53:03 +0200 Subject: [PATCH 5/6] feat: motiongram_data with vertical/horizontal orientation option Add a numpy-level motiongram core to _motionanalysis with a selectable orientation: the "vertical" (per-row mean, image row vs time) variant renders vertical trajectories -- e.g. a mallet's approach-and-rebound path -- directly, as used in the cymbal-comparison study; "horizontal" gives the per-column collapse. The existing image-producing pipelines (MgVideo.motiongrams and the ffmpeg videograms) already render both collapses as PNGs; what they do not offer is the motiongram as data on a chosen axis for further analysis, which is what the source study needed, so the orientation option is exposed here at the array level (with the mapping to the _mgh/_mgv outputs documented). A frame_diff=False mode yields the corresponding videogram. Tests: a synthetic falling-bar frame stack whose vertical motiongram ridge must descend monotonically, while its horizontal motiongram is uniform across columns. Source study: cymbal comparison (Jensenius). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/__init__.py | 1 + musicalgestures/_motionanalysis.py | 54 ++++++++++++++++++++++++++++++ tests/test_motiongram_data.py | 50 +++++++++++++++++++++++++++ 3 files changed, 105 insertions(+) create mode 100644 tests/test_motiongram_data.py diff --git a/musicalgestures/__init__.py b/musicalgestures/__init__.py index c7d5d10..bbf9919 100644 --- a/musicalgestures/__init__.py +++ b/musicalgestures/__init__.py @@ -108,3 +108,4 @@ def __init__(self): t60_backward_decay, attack_spectral_centroid, ) +from musicalgestures._motionanalysis import motiongram_data diff --git a/musicalgestures/_motionanalysis.py b/musicalgestures/_motionanalysis.py index 6a851a0..3ace1b6 100644 --- a/musicalgestures/_motionanalysis.py +++ b/musicalgestures/_motionanalysis.py @@ -2,6 +2,60 @@ import numpy as np +def motiongram_data(frames, orientation="vertical", frame_diff=True, normalize=True): + """ + Compute a motiongram as a plain numpy array from a stack of grayscale + frames, with a selectable orientation. + + With `orientation="vertical"` each (motion) frame is collapsed to its + per-row mean (the mean across image columns), and the resulting column + vectors are stacked over time into an (height, n) array -- image row vs + time. This "vertical approach" variant renders vertical trajectories + (e.g. a mallet's approach-and-rebound path toward an instrument) + directly. With `orientation="horizontal"` each frame is collapsed to + its per-column mean, giving a (width, n) array -- image column vs time + -- which renders horizontal (side-to-side) motion. + + This is the numpy-level counterpart of the image-producing motiongram + pipelines (`MgVideo.motiongrams`, whose `_mgh`/`_mgv` PNGs correspond to + the "vertical" and "horizontal" collapses here, up to transposition and + post-processing): use this function when you want the motiongram as + data for further analysis rather than as a rendered image. + + Source: cymbal-comparison study (Jensenius) -- vertical motiongram of + the mallet trajectory; building on the classic fourMs motiongram. + + Args: + frames (np.ndarray): Grayscale frames of shape (T, H, W). + orientation (str, optional): "vertical" (per-row mean; image row vs + time; shows vertical motion) or "horizontal" (per-column mean; + image column vs time; shows horizontal motion). Defaults to "vertical". + frame_diff (bool, optional): If True, collapse the absolute inter-frame + differences (a motiongram, T-1 time steps); if False, collapse the + frames themselves (a videogram, T time steps). Defaults to True. + normalize (bool, optional): If True, scale the result to [0, 1] by its + maximum. Defaults to True. + + Returns: + np.ndarray: The motiongram, of shape (H, T-1) for "vertical" or + (W, T-1) for "horizontal" (T instead of T-1 when `frame_diff` is + False). Time runs along the second axis. + """ + frames = np.asarray(frames, dtype=np.float32) + if frames.ndim != 3: + raise ValueError("motiongram_data expects frames of shape (T, H, W)") + data = np.abs(np.diff(frames, axis=0)) if frame_diff else frames + if orientation == "vertical": + gram = data.mean(axis=2).T # (H, T-1): image row vs time + elif orientation == "horizontal": + gram = data.mean(axis=1).T # (W, T-1): image column vs time + else: + raise ValueError("orientation must be 'vertical' or 'horizontal'") + if normalize: + gram = gram / (gram.max() + 1e-12) + return gram + + def centroid(image, width, height): """ Computes the centroid and quantity of motion in an image or frame. diff --git a/tests/test_motiongram_data.py b/tests/test_motiongram_data.py new file mode 100644 index 0000000..8df99a4 --- /dev/null +++ b/tests/test_motiongram_data.py @@ -0,0 +1,50 @@ +"""Tests for musicalgestures._motionanalysis.motiongram_data (orientation option).""" +import numpy as np +import pytest + +from musicalgestures import motiongram_data + + +def falling_bar_frames(T=30, H=40, W=60): + """A bright horizontal bar moving downward one row per frame.""" + frames = np.zeros((T, H, W), dtype=np.float32) + for i in range(T): + frames[i, 5 + i, :] = 255.0 + return frames + + +class TestMotiongramData: + def test_vertical_traces_falling_bar(self): + frames = falling_bar_frames() + gram = motiongram_data(frames, orientation="vertical") + assert gram.shape == (40, 29) # (H, T-1), time on axis 1 + # the motion energy at time t sits at rows 5+t / 6+t: the ridge descends + rows = gram.argmax(axis=0) + assert np.all(np.diff(rows) >= 0) + assert rows[0] in (5, 6) + assert rows[-1] in (33, 34) + + def test_horizontal_of_falling_bar_is_uniform(self): + frames = falling_bar_frames() + gram = motiongram_data(frames, orientation="horizontal") + assert gram.shape == (60, 29) # (W, T-1) + # a full-width bar moving vertically spreads evenly over columns + assert gram.std() / (gram.mean() + 1e-12) < 1e-6 + + def test_videogram_mode(self): + frames = falling_bar_frames() + gram = motiongram_data(frames, orientation="vertical", frame_diff=False) + assert gram.shape == (40, 30) # T columns, no differencing + + def test_normalization(self): + frames = falling_bar_frames() + gram = motiongram_data(frames) + assert gram.max() == pytest.approx(1.0) + raw = motiongram_data(frames, normalize=False) + assert raw.max() > 1.0 + + def test_bad_input_raises(self): + with pytest.raises(ValueError): + motiongram_data(np.zeros((5, 5))) + with pytest.raises(ValueError): + motiongram_data(falling_bar_frames(), orientation="diagonal") From 8c4a9ba43e153c4418550306252a4dd97edf8b81 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 09:17:16 +0200 Subject: [PATCH 6/6] fix: PR1 review findings (anchor one-to-one matching, guards, docstring accuracy) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/_alignment.py | 45 +++++++++++++++++++++++-------- musicalgestures/_audiofeatures.py | 16 ++++++----- musicalgestures/_pulse.py | 7 ++++- musicalgestures/_qom.py | 41 ++++++++++++++++++++++++---- tests/test_alignment.py | 31 +++++++++++++++++++-- tests/test_pulse.py | 4 +++ 6 files changed, 118 insertions(+), 26 deletions(-) diff --git a/musicalgestures/_alignment.py b/musicalgestures/_alignment.py index d8b8cc0..7d929ee 100644 --- a/musicalgestures/_alignment.py +++ b/musicalgestures/_alignment.py @@ -58,7 +58,7 @@ def xcorr_lag(x, y, fs, max_lag=1.5): # Near-exact ties are floating-point noise; prefer the smallest |lag|. cand = np.flatnonzero(cc >= cc.max() - 1e-9) best = cand[np.argmin(np.abs(lags[cand]))] - return lags[best] / fs, float(cc[best]) + return float(lags[best] / fs), float(cc[best]) def envelope_lag(x, y, rate, max_lag_s=1.5): @@ -142,13 +142,20 @@ def anchor_and_match(times_a, times_b, anchor_a=None, anchor_b=None, Per-take relative event alignment between two streams with independent clocks: both streams are shifted so that their anchor events -- by default the strongest event of each stream, physically the same moment - (e.g. the hardest strike) -- sit at t = 0; then every non-anchor event - of stream `a` is matched to the nearest event of stream `b` within + (e.g. the hardest strike) -- sit at t = 0; then non-anchor events of + stream `a` are matched to non-anchor events of stream `b` within +/- `window` seconds, and the signed offset (`b` minus `a`; positive = `b` later) is recorded. The anchor pair contributes an offset of 0 by - construction and is excluded. No absolute cross-stream synchronisation - is claimed: the result measures whether the two streams agree on the - RELATIVE timing of the remaining events. + construction and is excluded from BOTH streams -- `b`'s anchor is + dropped from the match pool symmetrically with `a`'s, so a non-anchor + `a` event near t = 0 cannot spuriously match `b`'s anchor. Matching is + one-to-one: every candidate pair within the window is considered in + order of increasing |offset| and greedily claimed, each `a` and `b` + event usable in at most one match (consistent with + `per_cycle_motion_delta`'s claim semantics), so a single `b` event + cannot be double-counted against multiple `a` events. No absolute + cross-stream synchronisation is claimed: the result measures whether + the two streams agree on the RELATIVE timing of the remaining events. The default matching window of 0.15 s is a PROVISIONAL default, reimplemented from the cymbal-comparison paper's method description. @@ -192,14 +199,30 @@ def pick_anchor(times, anchor, weights, name): b0 = pick_anchor(times_b, anchor_b, weights_b, "b") a = times_a - a0 b = times_b - b0 - a = a[np.abs(a) > 1e-12] # exclude the anchor itself + a = a[np.abs(a) > 1e-12] # exclude a's anchor + b = b[np.abs(b) > 1e-12] # exclude b's anchor, symmetrically if len(a) == 0 or len(b) == 0: return np.array([]) + + # All candidate pairs within the window, greedily claimed nearest-first + # so the match is one-to-one (an a event and a b event are each used + # in at most one pair), mirroring per_cycle_motion_delta's claiming. + ai, bi = np.meshgrid(np.arange(len(a)), np.arange(len(b)), indexing="ij") + diffs = b[bi] - a[ai] + within = np.abs(diffs) <= window + ai, bi, diffs = ai[within], bi[within], diffs[within] + order = np.argsort(np.abs(diffs)) + + claimed_a = np.zeros(len(a), dtype=bool) + claimed_b = np.zeros(len(b), dtype=bool) offsets = [] - for t in a: - j = int(np.argmin(np.abs(b - t))) - if abs(b[j] - t) <= window: - offsets.append(b[j] - t) + for k in order: + i, j = ai[k], bi[k] + if claimed_a[i] or claimed_b[j]: + continue + claimed_a[i] = True + claimed_b[j] = True + offsets.append(diffs[k]) return np.asarray(offsets, float) diff --git a/musicalgestures/_audiofeatures.py b/musicalgestures/_audiofeatures.py index 313478f..44644d2 100644 --- a/musicalgestures/_audiofeatures.py +++ b/musicalgestures/_audiofeatures.py @@ -139,13 +139,15 @@ def energy_onsets(y, sr, window=0.02, rel_threshold=0.15, min_interval=0.06): def t60_backward_decay(y, sr, window=0.02, spans=((-5, -35), (-5, -25)), rerise_db=6.0): """ - T60-style reverberation time by an ISO-3382-style backward-decay method: - the RMS envelope is converted to dB relative to its peak; from the peak - the decay is followed (stopping if the envelope re-rises by more than - `rerise_db` dB, marking a new onset), and T60 is estimated by linear - regression of the dB curve over the first available level span -- - by default -5 to -35 dB (a T30 measure, extrapolated x2), falling back - to -5 to -25 dB (T20, x3) when the deeper level is not reached. + T60-style reverberation time by an ISO-3382-inspired level-span + regression on the RMS envelope (without Schroeder backward + integration): the RMS envelope is converted to dB relative to its peak; + from the peak the decay is followed (stopping if the envelope re-rises + by more than `rerise_db` dB, marking a new onset), and T60 is estimated + by linear regression of the dB curve over the first available level + span -- by default -5 to -35 dB (a T30 measure, extrapolated x2), + falling back to -5 to -25 dB (T20, x3) when the deeper level is not + reached. The constants (20 ms window, -5/-35 with -5/-25 fallback, 6 dB re-rise stop) are PROVISIONAL defaults reimplemented from the cymbal-comparison diff --git a/musicalgestures/_pulse.py b/musicalgestures/_pulse.py index 3cb984f..ddf734c 100644 --- a/musicalgestures/_pulse.py +++ b/musicalgestures/_pulse.py @@ -249,11 +249,16 @@ def fit_accelerando(times, iois): tuple: `(ioi0, t_double, r2)` where `ioi0` is the fitted IOI at t=0 (s), `t_double` is the tempo-doubling time (s; `np.inf` if the sequence is not accelerating), and `r2` is the fit's coefficient of determination - on log2(IOI). + on log2(IOI). Degenerate input -- fewer than 3 valid (finite, + strictly positive) IOI points, too few to constrain the 2-parameter + fit -- returns `(nan, nan, nan)`. """ t = np.asarray(times, float) ioi = np.asarray(iois, float) m = np.isfinite(t) & np.isfinite(ioi) & (ioi > 0) + if m.sum() < 3: + nan = float("nan") + return nan, nan, nan t, ioi = t[m], np.log2(ioi[m]) A = np.vstack([t, np.ones_like(t)]).T (slope, intercept), *_ = np.linalg.lstsq(A, ioi, rcond=None) diff --git a/musicalgestures/_qom.py b/musicalgestures/_qom.py index 170da7a..a2f29ba 100644 --- a/musicalgestures/_qom.py +++ b/musicalgestures/_qom.py @@ -117,23 +117,39 @@ def band_limited_qom(pos, fs, lo=0.3, hi=15.0, order=4, auto_decimate=True): Returns: tuple: `(speed, fs_out)` where `speed` is the per-frame speed series - (length N-1, or shorter when decimated; empty if the input is too - short) and `fs_out` is its sampling rate (equal to `fs` unless - decimated). + (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 - hi_eff = min(hi, 0.9 * fs / 2) 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) @@ -303,12 +319,21 @@ def normalized_qom(landmarks, fs, scale=None, lo=0.3, hi=5.0, Returns: tuple: `(qom, speed, fs_out)` as in `group_qom`, with both `qom` and - `speed` divided by the body scale. + `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 @@ -335,8 +360,14 @@ def grid_qom(frames, grid=(6, 4), region=(0.0, 1.0, 0.0, 1.0), threshold=8.0): 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 diff --git a/tests/test_alignment.py b/tests/test_alignment.py index d0cfd6e..5064fed 100644 --- a/tests/test_alignment.py +++ b/tests/test_alignment.py @@ -56,10 +56,16 @@ def test_periodic_prefers_smallest_lag(self): assert r > 0.99 def test_envelope_lag_wrapper_agrees(self): + # Ground truth: y is built by rolling x by exactly 0.25 s, so the + # wrapper must recover that known lag (not merely equal whatever + # xcorr_lag happens to return). fs = 50.0 x = noisy_envelope(fs, seed=3) - y = np.roll(x, int(0.2 * fs)) - assert envelope_lag(x, y, fs) == xcorr_lag(x, y, fs) + shift = int(0.25 * fs) + y = np.roll(x, shift) + lag, r = envelope_lag(x[shift:-shift], y[shift:-shift], fs) + assert lag == pytest.approx(0.25, abs=1.5 / fs) + assert r > 0.9 class TestPerCycleMotionDelta: @@ -104,6 +110,27 @@ def test_requires_anchor_or_weights(self): with pytest.raises(ValueError): anchor_and_match([0.0, 1.0], [0.0, 1.0]) + def test_a_event_does_not_match_b_anchor(self): + # a has a non-anchor event within `window` of t = 0 (b's anchor, + # also at t = 0 after shifting); with no other b event nearby, this + # must NOT match -- b's anchor has to be excluded from the match + # pool symmetrically with a's, or this would spuriously match with + # a near-zero offset. + a = np.array([0.0, 0.05]) # anchor at 0.0, other event near it + b = np.array([0.0, 5.0]) # anchor at 0.0, no other nearby event + off = anchor_and_match(a, b, anchor_a=0.0, anchor_b=0.0, window=0.15) + assert len(off) == 0 + + def test_matching_is_one_to_one(self): + # Two a events both fall within `window` of the same single b + # event: only the nearer one may claim it, yielding exactly one + # match (not two, as many-to-one matching would give). + a = np.array([0.0, 1.0, 1.08]) + b = np.array([0.0, 1.02]) + off = anchor_and_match(a, b, anchor_a=0.0, anchor_b=0.0, window=0.15) + assert len(off) == 1 + assert off[0] == pytest.approx(0.02) # 1.02 - 1.0, the nearer pair + def test_offset_stats(self): off = np.array([-0.02, 0.0, 0.02, 0.04]) st = offset_stats(off) diff --git a/tests/test_pulse.py b/tests/test_pulse.py index 8c80f75..99f621b 100644 --- a/tests/test_pulse.py +++ b/tests/test_pulse.py @@ -108,6 +108,10 @@ def test_steady_pulse_gives_inf(self): _, t_double, _ = fit_accelerando(t, ioi) assert np.isinf(t_double) + def test_degenerate_input_gives_nan(self): + ioi0, t_double, r2 = fit_accelerando([0.0, 1.0], [2.0, 1.0]) + assert np.isnan(ioi0) and np.isnan(t_double) and np.isnan(r2) + class TestMotionOnsets: def test_recovers_ramp_onsets(self):