diff --git a/musicalgestures/_mocap.py b/musicalgestures/_mocap.py
index 5085fe3..dcd8b36 100644
--- a/musicalgestures/_mocap.py
+++ b/musicalgestures/_mocap.py
@@ -1,214 +1,26 @@
-"""
-Motion-capture I/O and cross-modality utilities.
-
-Pure numpy/scipy helpers ported from the "still standing" and
-Westney-comparisons studies:
+"""Re-exported from the ``micromotion`` package.
 
-* :func:`read_qtm_tsv` -- a single robust reader for Qualisys Track Manager
-  (QTM) TSV exports, consolidating the four/five near-duplicate loaders that
-  were copy-pasted across the study scripts.
-* :func:`compare_modality_envelopes` -- resample two motion envelopes onto a
-  common per-second grid and correlate them (e.g. video-pose vs mocap
-  validation).
-* :func:`dominant_frequency` -- the dominant spectral peak of a signal within
-  a band, via Welch.
+These functions used to live here. They were moved to ``micromotion`` on 2026-07-29 so that
+one implementation of quantity of motion exists rather than two, and MGT now depends on that
+package instead of carrying its own copy. Behaviour is unchanged: this module's tests pass
+against ``micromotion`` unmodified.
 
-.. note::
-   :func:`compare_modality_envelopes` deliberately takes *precomputed* 1-D
-   motion envelopes rather than computing quantity-of-motion internally, so
-   this module stays independent of the QoM machinery. The natural producer
-   of such envelopes is ``musicalgestures._qom.band_limited_qom`` followed by
-   a per-second binning (``envelope`` / ``bin_series``), arriving in a
-   sibling PR.
+The dependency points this way round on purpose. ``micromotion`` needs only numpy, scipy and
+pandas, so someone analysing accelerometer data does not have to install a computer-vision
+stack; MGT already depends on ``ambiscape`` the same way, and neither of those packages
+imports MGT.
 
-Source: still standing study and Westney-comparisons study (Jensenius).
+Import from ``micromotion`` directly in new code.
 """
 
-import numpy as np
-
-
-def read_qtm_tsv(path):
-    """
-    Read a Qualisys Track Manager (QTM) TSV motion-capture export.
-
-    Consolidates the several near-duplicate loaders used across the studies
-    into one robust reader. It locates the ``MARKER_NAMES`` header row to
-    recover marker labels, autodetects where the numeric data block starts
-    (the first row whose first field parses as a float), drops a trailing
-    all-empty column produced by a trailing tab, converts exact-zero XYZ
-    triples (Qualisys gap fills) to ``NaN``, and falls back from UTF-8 to
-    latin-1 encoding. When a ``FREQUENCY`` header field is present the frame
-    rate is returned as well.
-
-    Source: still standing study and Westney-comparisons study (Jensenius) --
-    unifies the ``load_qtm`` variants in the balance/dynamics/circular/
-    spatial-range reports and the latin-1 variant in ``compare_mp_mocap``.
-
-    Args:
-        path (str): Path to the ``.tsv`` file.
-
-    Returns:
-        tuple: ``(marker_names, data, fs)`` where ``marker_names`` is a list
-            of ``M`` strings (empty if no header was found), ``data`` is a
-            float array of shape ``(T, M, 3)`` with gaps as ``NaN``, and
-            ``fs`` is the frame rate in Hz or ``None`` if not derivable from
-            the header.
-    """
-    import io
-
-    def _read_lines(enc):
-        with io.open(path, encoding=enc) as fh:
-            return fh.readlines()
-
-    try:
-        lines = _read_lines("utf-8")
-    except (UnicodeDecodeError, UnicodeError):
-        lines = _read_lines("latin-1")
-
-    marker_names = []
-    fs = None
-    data_start = None
-    for i, ln in enumerate(lines):
-        parts = ln.rstrip("\n").split("\t")
-        key = parts[0].strip().upper()
-        if key == "MARKER_NAMES":
-            marker_names = [p.strip() for p in parts[1:] if p.strip()]
-            continue
-        if key == "FREQUENCY" and len(parts) > 1:
-            try:
-                fs = float(parts[1])
-            except ValueError:
-                pass
-            continue
-        # first row whose leading field is numeric marks the data block
-        try:
-            float(parts[0])
-            data_start = i
-            break
-        except ValueError:
-            continue
-
-    if data_start is None:
-        raise ValueError(f"No numeric data block found in {path!r}")
-
-    rows = []
-    for ln in lines[data_start:]:
-        parts = ln.rstrip("\n").split("\t")
-        try:
-            rows.append([float(p) for p in parts if p != ""])
-        except ValueError:
-            continue
-    if not rows:
-        raise ValueError(f"No parseable data rows in {path!r}")
-
-    ncol = min(len(r) for r in rows)
-    arr = np.array([r[:ncol] for r in rows], dtype=float)
-
-    # A QTM data row may carry a leading time/frame column(s); the marker
-    # block is the trailing 3*M columns. Prefer the marker count from the
-    # header; otherwise infer the largest multiple of 3.
-    if marker_names:
-        M = len(marker_names)
-    else:
-        M = ncol // 3
-    offset = arr.shape[1] - 3 * M
-    if offset < 0:
-        # header lists more markers than columns present; fall back
-        M = arr.shape[1] // 3
-        offset = arr.shape[1] - 3 * M
-        marker_names = marker_names[:M]
-    block = arr[:, offset:offset + 3 * M]
-    data = block.reshape(block.shape[0], M, 3)
-
-    # exact-zero triples are Qualisys gap fills -> NaN
-    gap = np.all(data == 0, axis=2)
-    data[gap] = np.nan
-
-    return marker_names, data, fs
-
-
-def dominant_frequency(x, fs, band=(0.3, 4.0)):
-    """
-    Dominant frequency of a signal within a band, via a Welch spectrum.
-
-    Returns the frequency of the largest Welch power-spectral-density peak
-    inside ``band`` -- e.g. the dominant oscillation rate of a body-part
-    speed or vertical-position signal.
-
-    Source: Westney-comparisons study (Jensenius), extended motion-feature
-    analysis (motion dominant frequency of vertical trunk position).
-
-    Args:
-        x (np.ndarray): 1-D input signal.
-        fs (float): Sampling rate in Hz.
-        band (tuple, optional): ``(low, high)`` search band in Hz. Defaults
-            to ``(0.3, 4.0)``.
-
-    Returns:
-        float: The dominant frequency in Hz, or ``nan`` if the band is empty.
-    """
-    from scipy.signal import welch
-
-    x = np.asarray(x, dtype=float)
-    x = x[np.isfinite(x)]
-    if len(x) < 8:
-        return np.nan
-    lo, hi = band
-    f, P = welch(x - x.mean(), fs, nperseg=min(2048, len(x)))
-    mask = (f >= lo) & (f <= hi)
-    if not mask.any():
-        return np.nan
-    return float(f[mask][np.argmax(P[mask])])
-
-
-def compare_modality_envelopes(env_a, env_b, fs_a, fs_b):
-    """
-    Correlate two motion envelopes after resampling to a common grid.
-
-    Resamples both 1-D envelopes onto a shared one-sample-per-second grid
-    (by averaging within each second), truncates to the common length, and
-    returns the Pearson correlation -- the video-vs-mocap (or view-vs-view)
-    agreement measure. Both inputs are treated as already-computed motion
-    envelopes (e.g. per-frame band-limited quantity-of-motion), keeping this
-    function decoupled from the QoM computation itself. The per-second binning
-    uses an integer-rounded step, so non-integer frame rates (e.g. 29.97 fps)
-    drift slightly over long signals; this function is intended for validation
-    rather than precise alignment.
-
-    Source: still standing / Westney-comparisons study (Jensenius),
-    MediaPipe-vs-mocap validation (``compare_mp_mocap``).
-
-    Args:
-        env_a (np.ndarray): First 1-D motion envelope.
-        env_b (np.ndarray): Second 1-D motion envelope.
-        fs_a (float): Sampling rate of ``env_a`` in Hz.
-        fs_b (float): Sampling rate of ``env_b`` in Hz.
-
-    Returns:
-        dict: ``{"r", "n"}`` where ``r`` is the Pearson correlation of the
-            two per-second envelopes and ``n`` the number of common seconds
-            (``r`` is ``nan`` if fewer than three overlapping seconds or if
-            either resampled envelope is constant).
-    """
-    def _per_second(env, fs):
-        env = np.asarray(env, dtype=float)
-        nsec = int(len(env) / fs)
-        if nsec < 1:
-            return np.array([])
-        step = int(round(fs))
-        step = max(step, 1)
-        return np.array([np.nanmean(env[i * step:(i + 1) * step])
-                         for i in range(nsec)])
-
-    a = _per_second(env_a, fs_a)
-    b = _per_second(env_b, fs_b)
-    n = min(len(a), len(b))
-    if n < 3:
-        return dict(r=np.nan, n=n)
-    a = a[:n]
-    b = b[:n]
-    ok = np.isfinite(a) & np.isfinite(b)
-    if ok.sum() < 3 or a[ok].std() == 0 or b[ok].std() == 0:
-        return dict(r=np.nan, n=int(ok.sum()))
-    r = float(np.corrcoef(a[ok], b[ok])[0, 1])
-    return dict(r=r, n=int(ok.sum()))
+from micromotion.mocap import (  # noqa: F401
+    compare_modality_envelopes,
+    dominant_frequency,
+    read_qtm_tsv,
+)
+
+__all__ = [
+    "compare_modality_envelopes",
+    "dominant_frequency",
+    "read_qtm_tsv",
+]
diff --git a/musicalgestures/_physio.py b/musicalgestures/_physio.py
index 79115d4..c2eda23 100644
--- a/musicalgestures/_physio.py
+++ b/musicalgestures/_physio.py
@@ -1,141 +1,24 @@
-"""
-Physiology signal features for standstill / micromotion studies.
+"""Re-exported from the ``micromotion`` package.
 
-Two pure numpy/scipy surfaces ported from the "still standing" study:
+These functions used to live here. They were moved to ``micromotion`` on 2026-07-29 so that
+one implementation of quantity of motion exists rather than two, and MGT now depends on that
+package instead of carrying its own copy. Behaviour is unchanged: this module's tests pass
+against ``micromotion`` unmodified.
 
-* :func:`respiration_rate` -- windowed breathing rate (breaths per minute)
-  from a respiration waveform, via band-pass filtering and a Welch spectral
-  peak per window.
-* :func:`spectral_band_fractions` -- the fraction of a signal's Welch power
-  falling in each of a set of caller-supplied named frequency bands. This is
-  the generic "cardiorespiratory QoM" spectral-composition diagnostic with
-  the heart-rate/respiration bands supplied by the caller, so the function
-  carries no dependency on any particular physiological sensor.
+The dependency points this way round on purpose. ``micromotion`` needs only numpy, scipy and
+pandas, so someone analysing accelerometer data does not have to install a computer-vision
+stack; MGT already depends on ``ambiscape`` the same way, and neither of those packages
+imports MGT.
 
-Source: still standing study (Jensenius) -- Deichman / Equivital physiology
-analyses.
+Import from ``micromotion`` directly in new code.
 """
 
-import numpy as np
-
-
-def respiration_rate(waveform, fs, *, band=(0.1, 0.6), window_s=30,
-                     step_s=30):
-    """
-    Windowed respiration rate (breaths per minute) from a breathing waveform.
-
-    Each analysis window is band-pass filtered to the respiration band and
-    its dominant frequency is taken as the Welch spectral peak inside that
-    band; the rate is that frequency times 60. Windows advance by ``step_s``
-    seconds. The default band ``(0.1, 0.6)`` Hz corresponds to about
-    6-36 breaths/min. Each window must contain at least 15 seconds of valid
-    samples for spectral estimation.
-
-    Source: still standing study (Jensenius), Deichman respiration analysis
-    (``compute_qom_resp``).
-
-    Args:
-        waveform (np.ndarray): 1-D respiration/breathing waveform.
-        fs (float): Sampling rate in Hz.
-        band (tuple, optional): ``(low, high)`` respiration band in Hz.
-            Defaults to ``(0.1, 0.6)``.
-        window_s (float, optional): Window length in seconds. Defaults to 30.
-        step_s (float, optional): Hop between windows in seconds. Defaults to
-            30.
-
-    Returns:
-        dict: ``{"rate_bpm", "times_s", "median_bpm"}`` where ``rate_bpm`` is
-            the per-window rate (breaths/min, ``nan`` for windows without a
-            clear peak), ``times_s`` the window centre times in seconds, and
-            ``median_bpm`` the median across valid windows.
-    """
-    from scipy.signal import butter, filtfilt, welch
-
-    x = np.asarray(waveform, dtype=float)
-    x = x[np.isfinite(x)]
-    lo, hi = band
-    nyq = fs / 2.0
-    if len(x) < int(fs * window_s):
-        # single short window: still attempt one estimate
-        window_s = max(1.0, len(x) / fs)
-
-    b, a = butter(2, [lo / nyq, hi / nyq], btype="band")
-    xf = filtfilt(b, a, x - x.mean())
-
-    win = int(fs * window_s)
-    step = int(fs * step_s)
-    win = max(win, 1)
-    step = max(step, 1)
-
-    rates = []
-    times = []
-    for start in range(0, max(len(xf) - win + 1, 1), step):
-        seg = xf[start:start + win]
-        if len(seg) < int(fs * 15):  # need at least 15 s for spectral estimation
-            rates.append(np.nan)
-            times.append((start + win / 2) / fs)
-            continue
-        nperseg = min(len(seg), int(fs * window_s))
-        f, P = welch(seg, fs, nperseg=nperseg)
-        mask = (f >= lo) & (f <= hi)
-        if mask.any() and P[mask].sum() > 0:
-            fpk = f[mask][np.argmax(P[mask])]
-            rates.append(float(fpk * 60.0))
-        else:
-            rates.append(np.nan)
-        times.append((start + win / 2) / fs)
-
-    rate_bpm = np.array(rates, dtype=float)
-    return dict(rate_bpm=rate_bpm, times_s=np.array(times, dtype=float),
-                median_bpm=float(np.nanmedian(rate_bpm))
-                if np.isfinite(rate_bpm).any() else np.nan)
-
-
-def spectral_band_fractions(signal, fs, bands, *, total_band=(0.1, 8.0),
-                            nperseg_s=20):
-    """
-    Fraction of a signal's power in each of a set of named frequency bands.
-
-    Estimates the Welch power spectrum and, for each named band in ``bands``,
-    returns that band's summed power divided by the summed power in
-    ``total_band``. This is the generic spectral-composition diagnostic used
-    for the "cardiorespiratory QoM artifact" analysis (e.g. how much of a
-    chest-accelerometer QoM signal sits in a cardiac vs a respiration band),
-    with the bands supplied by the caller so there is no built-in dependence
-    on a heart-rate or respiration sensor. Power is bin-summed on the Welch
-    grid; the study source integrated with trapz, which yields nearly
-    identical results on the uniform frequency spacing of Welch.
-
-    Source: still standing study (Jensenius), Deichman chest-QoM
-    cardiorespiratory spectral-composition analysis (``deichman_full``).
-
-    Args:
-        signal (np.ndarray): 1-D input signal.
-        fs (float): Sampling rate in Hz.
-        bands (dict): Mapping of band name to ``(low, high)`` in Hz, e.g.
-            ``{"cardiac": (0.9, 1.3), "resp": (0.12, 0.5)}``.
-        total_band (tuple, optional): ``(low, high)`` reference band whose
-            power is the denominator. Defaults to ``(0.1, 8.0)``.
-        nperseg_s (float, optional): Welch segment length in seconds.
-            Defaults to 20.
-
-    Returns:
-        dict: Mapping of each band name to its power fraction in ``[0, 1]``
-            (``nan`` if the total band contains no power).
-    """
-    from scipy.signal import welch
+from micromotion.physio import (  # noqa: F401
+    respiration_rate,
+    spectral_band_fractions,
+)
 
-    x = np.asarray(signal, dtype=float)
-    x = x[np.isfinite(x)]
-    if len(x) < 8:
-        return {name: np.nan for name in bands}
-    nperseg = min(len(x), max(8, int(fs * nperseg_s)))
-    f, P = welch(x - x.mean(), fs, nperseg=nperseg)
-    tlo, thi = total_band
-    total = P[(f >= tlo) & (f < thi)].sum()
-    if total <= 0:
-        return {name: np.nan for name in bands}
-    out = {}
-    for name, (lo, hi) in bands.items():
-        out[name] = float(P[(f >= lo) & (f < hi)].sum() / total)
-    return out
+__all__ = [
+    "respiration_rate",
+    "spectral_band_fractions",
+]
diff --git a/musicalgestures/_posture.py b/musicalgestures/_posture.py
index 4751bc0..88ee353 100644
--- a/musicalgestures/_posture.py
+++ b/musicalgestures/_posture.py
@@ -1,633 +1,44 @@
-"""
-Posturography and standstill-sway metrics for centre-of-pressure (CoP) and
-head/marker position signals.
-
-This module ports the "still standing" study's posturography stack into
-pure numpy/scipy surfaces that operate on plain arrays -- no study-specific
-loaders, axis conventions, or marker loops. Three families of measures are
-provided:
+"""Re-exported from the ``micromotion`` package.
 
-* **Sway amount / geometry** -- :func:`cop_sway_metrics`,
-  :func:`confidence_ellipse_area`, :func:`convex_hull_area`.
-* **Control dynamics / complexity** -- :func:`stabilogram_diffusion`
-  (Collins-De Luca SDA), :func:`dfa` (detrended fluctuation analysis),
-  :func:`sample_entropy`, :func:`spectral_edges`, :func:`sway_texture`,
-  :func:`principal_axis_projection`.
-* **Direction / extent** -- :func:`sway_orientation`, :func:`axial_rayleigh`,
-  :func:`spatial_extent`.
+These functions used to live here. They were moved to ``micromotion`` on 2026-07-29 so that
+one implementation of quantity of motion exists rather than two, and MGT now depends on that
+package instead of carrying its own copy. Behaviour is unchanged: this module's tests pass
+against ``micromotion`` unmodified.
 
-The from-scratch SDA / DFA / sample-entropy implementations are validated in
-the test-suite against known-answer synthetic signals (white noise ->
-DFA alpha ~= 0.5 and SDA Hurst ~= 0.5; a sine -> low sample entropy relative
-to its shuffle).
+The dependency points this way round on purpose. ``micromotion`` needs only numpy, scipy and
+pandas, so someone analysing accelerometer data does not have to install a computer-vision
+stack; MGT already depends on ``ambiscape`` the same way, and neither of those packages
+imports MGT.
 
-Source: still standing study (Jensenius) -- posturography and micromotion
-analyses of the international "standstill" championships and related datasets.
+Import from ``micromotion`` directly in new code.
 """
 
-import numpy as np
-
-
-# ---------------------------------------------------------------------------
-# Sway amount / geometry
-# ---------------------------------------------------------------------------
-def confidence_ellipse_area(xy, conf=0.95):
-    """
-    Area of the confidence ellipse of a 2-D point cloud (e.g. a
-    centre-of-pressure trace).
-
-    The ellipse is the standard bivariate-Gaussian confidence region
-    ``area = pi * chi2_conf,2df * sqrt(det Cov)`` where ``Cov`` is the
-    2x2 covariance of the (mean-removed) points. For a CoP sway path this
-    is the classic 95% "sway-ellipse area".
-
-    Source: still standing study (Jensenius), HpSp balance analysis.
-
-    Args:
-        xy (np.ndarray): Point cloud of shape ``(T, 2)``.
-        conf (float, optional): Confidence level in ``(0, 1)``. Defaults to
-            0.95.
-
-    Returns:
-        float: Ellipse area in squared position units (e.g. mm^2), or
-            ``nan`` if fewer than three finite points are available.
-    """
-    from scipy.stats import chi2
-
-    xy = np.asarray(xy, dtype=float)
-    xy = xy[np.isfinite(xy).all(axis=1)]
-    if len(xy) < 3:
-        return np.nan
-    cov = np.cov(xy.T)
-    return float(np.pi * chi2.ppf(conf, 2) * np.sqrt(max(np.linalg.det(cov), 0.0)))
-
-
-def convex_hull_area(xy):
-    """
-    Area of the 2-D convex hull of a point cloud.
-
-    A non-parametric alternative to :func:`confidence_ellipse_area` for the
-    region occupied by a sway path: it makes no Gaussian assumption and is
-    driven by the outermost excursions.
-
-    Source: still standing study (Jensenius); complements the confidence
-    ellipse used in the balance reports.
-
-    Args:
-        xy (np.ndarray): Point cloud of shape ``(T, 2)``.
-
-    Returns:
-        float: Convex-hull area in squared position units, or ``nan`` if
-            fewer than three non-collinear finite points are available.
-    """
-    from scipy.spatial import ConvexHull, QhullError
-
-    xy = np.asarray(xy, dtype=float)
-    xy = xy[np.isfinite(xy).all(axis=1)]
-    if len(xy) < 3:
-        return np.nan
-    try:
-        return float(ConvexHull(xy).volume)  # 2-D "volume" == area
-    except QhullError:
-        return np.nan
-
-
-def cop_sway_metrics(xy, t=None, fs=None, *, freq_band=(0.1, 5.0),
-                     resample_fs=50.0):
-    """
-    Standard centre-of-pressure (CoP) sway metrics from a 2-D sway path.
-
-    Computes the classic posturographic descriptors: CoP path length and
-    path rate, the 95% confidence-ellipse area, medio-lateral (ML) and
-    antero-posterior (AP) ranges and standard deviations, the AP/ML range
-    and SD ratios, and the mean sway frequency of each axis (the
-    power-weighted mean frequency of a Welch spectrum inside ``freq_band``,
-    computed on a uniform grid at ``resample_fs``).
-
-    The first column of ``xy`` is treated as ML and the second as AP,
-    matching the study convention. Sampling time may be given either as an
-    explicit time vector ``t`` (seconds; may be irregular) or a constant
-    rate ``fs`` (Hz); if neither is supplied a rate of 1 Hz is assumed.
-
-    Source: still standing study (Jensenius), HpSp balance analysis
-    (``analyze_balance``).
-
-    Args:
-        xy (np.ndarray): CoP path of shape ``(T, 2)`` as ``[ML, AP]`` in
-            position units (e.g. mm).
-        t (np.ndarray, optional): Per-sample timestamps in seconds. May be
-            irregular. Defaults to None.
-        fs (float, optional): Constant sampling rate in Hz, used when ``t``
-            is not given. Defaults to None (interpreted as 1 Hz).
-        freq_band (tuple, optional): ``(low, high)`` band in Hz for the mean
-            sway frequency. Defaults to ``(0.1, 5.0)``.
-        resample_fs (float, optional): Uniform rate in Hz onto which the
-            path is interpolated before the spectral estimate. Defaults to
-            50.0.
-
-    Returns:
-        dict: Metrics with keys ``n``, ``dur``, ``fs_mean``, ``path_len``,
-            ``path_rate``, ``area95``, ``ml_range``, ``ap_range``,
-            ``ml_sd``, ``ap_sd``, ``ap_ml_range_ratio``,
-            ``ap_ml_sd_ratio``, ``mf_ml``, ``mf_ap`` and ``mf_mean``.
-    """
-    from scipy import signal as _sig
-
-    xy = np.asarray(xy, dtype=float)
-    if xy.ndim != 2 or xy.shape[1] != 2:
-        raise ValueError("xy must have shape (T, 2)")
-    ml = xy[:, 0]
-    ap = xy[:, 1]
-    n = len(xy)
-
-    if t is not None:
-        t = np.asarray(t, dtype=float)
-        t = t - t[0]
-    elif fs is not None:
-        t = np.arange(n) / float(fs)
-    else:
-        t = np.arange(n, dtype=float)
-    dur = float(t[-1]) if n > 1 else 0.0
-    fs_mean = (n - 1) / dur if dur > 0 else np.nan
-
-    ml_c = ml - np.nanmean(ml)
-    ap_c = ap - np.nanmean(ap)
-
-    dpath = np.sqrt(np.diff(ml) ** 2 + np.diff(ap) ** 2)
-    path_len = float(np.nansum(dpath))
-    path_rate = path_len / dur if dur > 0 else np.nan
-
-    ml_range = float(np.nanmax(ml) - np.nanmin(ml))
-    ap_range = float(np.nanmax(ap) - np.nanmin(ap))
-    ml_sd = float(np.nanstd(ml, ddof=1))
-    ap_sd = float(np.nanstd(ap, ddof=1))
-
-    area95 = confidence_ellipse_area(np.column_stack([ml_c, ap_c]), conf=0.95)
-
-    lo, hi = freq_band
-    tu = np.arange(0, dur, 1.0 / resample_fs) if dur > 0 else np.array([])
-
-    def _mean_freq(x):
-        if len(tu) < 4:
-            return np.nan
-        xu = np.interp(tu, t, x)
-        xu = _sig.detrend(xu)
-        nperseg = min(len(xu), int(resample_fs * 20))
-        if nperseg < 4:
-            return np.nan
-        f, P = _sig.welch(xu, fs=resample_fs, nperseg=nperseg)
-        band = (f >= lo) & (f <= hi)
-        if P[band].sum() <= 0:
-            return np.nan
-        return float(np.sum(f[band] * P[band]) / np.sum(P[band]))
-
-    mf_ml = _mean_freq(ml_c)
-    mf_ap = _mean_freq(ap_c)
-
-    return dict(
-        n=n, dur=dur, fs_mean=fs_mean,
-        path_len=path_len, path_rate=path_rate, area95=area95,
-        ml_range=ml_range, ap_range=ap_range, ml_sd=ml_sd, ap_sd=ap_sd,
-        ap_ml_range_ratio=ap_range / ml_range if ml_range else np.nan,
-        ap_ml_sd_ratio=ap_sd / ml_sd if ml_sd else np.nan,
-        mf_ml=mf_ml, mf_ap=mf_ap, mf_mean=float(np.nanmean([mf_ml, mf_ap])),
-    )
-
-
-# ---------------------------------------------------------------------------
-# Control dynamics / complexity
-# ---------------------------------------------------------------------------
-def principal_axis_projection(xy):
-    """
-    Project a 2-D (or N-D) point cloud onto its principal axis.
-
-    Runs a PCA on the mean-removed points and returns the 1-D coordinate
-    along the direction of greatest variance -- the natural 1-D reduction of
-    a sway path used by the dynamics/complexity measures. The PCA eigenvector
-    sign is arbitrary; the projection may be globally flipped across calls or
-    datasets.
-
-    Source: still standing study (Jensenius), sway-dynamics analysis.
-
-    Args:
-        xy (np.ndarray): Point cloud of shape ``(T, D)`` (typically
-            ``D == 2``).
-
-    Returns:
-        np.ndarray: 1-D projection of shape ``(T,)``.
-    """
-    xy = np.asarray(xy, dtype=float)
-    xy = xy - xy.mean(axis=0)
-    cov = np.cov(xy.T)
-    w, v = np.linalg.eigh(cov)
-    return xy @ v[:, -1]  # eigvecs ascending -> last is major axis
-
-
-def stabilogram_diffusion(xy, fs, *, short_max_s=0.6, long_min_s=1.5,
-                          n_lags=40):
-    """
-    Collins-De Luca stabilogram-diffusion analysis (SDA) of a sway path.
-
-    Fits the mean-square-displacement (MSD) curve
-    ``<[r(t+dt) - r(t)]^2>`` versus time-lag ``dt`` in log-log space and
-    reports a short-term and a long-term Hurst exponent (each ``slope / 2``)
-    plus the critical crossover time where the two regression lines
-    intersect. Persistent (open-loop) drift gives a short-term Hurst above
-    0.5; anti-persistent (closed-loop correction) gives a long-term Hurst
-    below 0.5.
-
-    Source: still standing study (Jensenius), sway-dynamics analysis;
-    method of Collins & De Luca (1993).
-
-    Args:
-        xy (np.ndarray): Sway path of shape ``(T, D)`` (``D >= 1``). A 1-D
-            input of shape ``(T,)`` is accepted and treated as a single
-            axis.
-        fs (float): Sampling rate in Hz.
-        short_max_s (float, optional): Upper bound (s) of the short-term
-            fitting window. Defaults to 0.6.
-        long_min_s (float, optional): Lower bound (s) of the long-term
-            fitting window. Defaults to 1.5.
-        n_lags (int, optional): Number of log-spaced lags at which the MSD
-            is evaluated. Defaults to 40.
-
-    Returns:
-        dict: ``{"H_short", "H_long", "crossover_s"}``. Entries are ``nan``
-            when a window contains fewer than three usable lags.
-    """
-    xy = np.asarray(xy, dtype=float)
-    if xy.ndim == 1:
-        xy = xy[:, None]
-    n = len(xy)
-    if n < 8:
-        return dict(H_short=np.nan, H_long=np.nan, crossover_s=np.nan)
-
-    lags = np.unique(
-        np.round(np.logspace(0, np.log10(max(n // 4, 2)), n_lags)).astype(int))
-    lags = lags[lags >= 1]
-    msd = np.array([np.mean(np.sum((xy[l:] - xy[:-l]) ** 2, axis=1))
-                    for l in lags])
-    t = lags / fs
-    ok = msd > 0
-    t, msd = t[ok], msd[ok]
-    short = t < short_max_s
-    long = t > long_min_s
-
-    def _slope_intercept(mask):
-        if mask.sum() < 3:
-            return np.nan, np.nan
-        a, b = np.polyfit(np.log(t[mask]), np.log(msd[mask]), 1)
-        return a, b
-
-    a_s, b_s = _slope_intercept(short)
-    a_l, b_l = _slope_intercept(long)
-    H_short = a_s / 2.0 if np.isfinite(a_s) else np.nan
-    H_long = a_l / 2.0 if np.isfinite(a_l) else np.nan
-    crossover = np.nan
-    if np.isfinite(a_s) and np.isfinite(a_l) and abs(a_s - a_l) > 1e-9:
-        crossover = float(np.exp((b_l - b_s) / (a_s - a_l)))
-    return dict(H_short=float(H_short), H_long=float(H_long),
-                crossover_s=crossover)
-
-
-def dfa(x, *, n_scales=18, min_scale=10):
-    """
-    Detrended fluctuation analysis (DFA) scaling exponent.
-
-    Integrates the mean-removed signal, then measures the RMS of the
-    linearly-detrended integrated profile within non-overlapping windows of
-    increasing size; the slope of ``log F(n)`` versus ``log n`` is the DFA
-    exponent ``alpha``. White noise gives ``alpha ~= 0.5``; a random walk
-    (Brownian) gives ``alpha ~= 1.5``; ``alpha == 1`` is 1/f noise.
-
-    Source: still standing study (Jensenius), sway-complexity analysis;
-    method of Peng et al. (1994).
-
-    Args:
-        x (np.ndarray): 1-D input signal.
-        n_scales (int, optional): Number of log-spaced window sizes.
-            Defaults to 18.
-        min_scale (int, optional): Smallest window size in samples. Defaults
-            to 10.
-
-    Returns:
-        float: The DFA exponent ``alpha`` (``nan`` if too short).
-    """
-    x = np.asarray(x, dtype=float)
-    x = x[np.isfinite(x)]
-    n = len(x)
-    if n < 4 * min_scale:
-        return np.nan
-    y = np.cumsum(x - x.mean())
-    scales = np.unique(
-        np.round(np.logspace(np.log10(min_scale),
-                             np.log10(n // 4), n_scales)).astype(int))
-    F = []
-    used = []
-    for m in scales:
-        k = n // m
-        if k < 1:
-            continue
-        seg = y[:k * m].reshape(k, m)
-        tt = np.arange(m)
-        rms = [np.sqrt(np.mean((s - np.polyval(np.polyfit(tt, s, 1), tt)) ** 2))
-               for s in seg]
-        fm = np.mean(rms)
-        if fm > 0:
-            F.append(fm)
-            used.append(m)
-    if len(F) < 2:
-        return np.nan
-    return float(np.polyfit(np.log(used), np.log(F), 1)[0])
-
-
-def sample_entropy(x, m=2, r=0.2):
-    """
-    Sample entropy (SampEn) of a 1-D signal.
-
-    Measures regularity/predictability: the negative log conditional
-    probability that sequences close (within tolerance ``r``) for ``m``
-    samples remain close for ``m + 1`` samples, self-matches excluded.
-    Lower values mean more repetitive/predictable signals. The signal is
-    z-scored internally so ``r`` is expressed as a fraction of its standard
-    deviation. Neighbour counts use a Chebyshev (max-norm) KD-tree.
-
-    Source: still standing study (Jensenius), sway-complexity analysis;
-    method of Richman & Moorman (2000).
-
-    Args:
-        x (np.ndarray): 1-D input signal.
-        m (int, optional): Embedding (template) length. Defaults to 2.
-        r (float, optional): Tolerance as a fraction of the signal SD.
-            Defaults to 0.2.
-
-    Returns:
-        float: Sample entropy (``nan`` if undefined, e.g. no matches).
-    """
-    from scipy.spatial import cKDTree
-
-    x = np.asarray(x, dtype=float)
-    x = x[np.isfinite(x)]
-    N = len(x)
-    if N < m + 2:
-        return np.nan
-    x = (x - x.mean()) / (x.std() + 1e-12)
-
-    def _count(mm):
-        emb = np.array([x[i:i + mm] for i in range(N - mm + 1)])
-        tree = cKDTree(emb)
-        pairs = tree.query_ball_point(emb, r, p=np.inf)
-        return sum(len(p) - 1 for p in pairs)  # exclude self-match
-
-    B = _count(m)
-    A = _count(m + 1)
-    if B == 0 or A == 0:
-        return np.nan
-    return float(-np.log(A / B))
-
-
-def spectral_edges(x, fs, *, edges=(0.5, 0.95), nperseg=None):
-    """
-    Spectral-edge frequencies of a signal.
-
-    Returns the frequencies below which a given cumulative fraction of the
-    Welch power spectrum lies. With the default ``edges`` the first value is
-    the median frequency (50% edge) and the second the 95% spectral-edge
-    frequency, two standard descriptors of sway spectral shape.
-
-    Source: still standing study (Jensenius), sway-dynamics analysis.
-
-    Args:
-        x (np.ndarray): 1-D input signal.
-        fs (float): Sampling rate in Hz.
-        edges (tuple, optional): Cumulative-power fractions in ``(0, 1)``.
-            Defaults to ``(0.5, 0.95)``.
-        nperseg (int, optional): Welch segment length in samples. Defaults
-            to ``min(2048, len(x))``.
-
-    Returns:
-        dict: Mapping ``"f<pct>"`` (e.g. ``"f50"``, ``"f95"``) to the edge
-            frequency in Hz.
-    """
-    from scipy import signal as _sig
-
-    x = np.asarray(x, dtype=float)
-    x = x[np.isfinite(x)]
-    if len(x) < 8:
-        return {f"f{int(round(e * 100))}": np.nan for e in edges}
-    if nperseg is None:
-        nperseg = min(2048, len(x))
-    f, P = _sig.welch(x - x.mean(), fs, nperseg=nperseg)
-    cP = np.cumsum(P)
-    if cP[-1] <= 0:
-        return {f"f{int(round(e * 100))}": np.nan for e in edges}
-    cP = cP / cP[-1]
-    return {f"f{int(round(e * 100))}": float(np.interp(e, cP, f))
-            for e in edges}
-
-
-def sway_texture(speed, fs, *, frozen_threshold=2.0):
-    """
-    Micro-texture of a sway speed signal: frozen fraction and burst rate.
-
-    Distinguishes a smooth wander from intermittent ballistic corrections.
-    The frozen fraction is the share of time the speed is below
-    ``frozen_threshold``; the burst rate is the number of upward threshold
-    crossings (onset of a velocity burst) per minute.
-
-    Source: still standing study (Jensenius), sway-texture analysis.
-
-    Args:
-        speed (np.ndarray): 1-D speed signal (e.g. mm/s).
-        fs (float): Sampling rate in Hz.
-        frozen_threshold (float, optional): Speed below which the signal is
-            considered "frozen", in the units of ``speed``. Defaults to 2.0.
-
-    Returns:
-        dict: ``{"frozen_fraction", "burst_rate"}`` where ``burst_rate`` is
-            in bursts per minute.
-    """
-    speed = np.asarray(speed, dtype=float)
-    speed = speed[np.isfinite(speed)]
-    if len(speed) < 2:
-        return dict(frozen_fraction=np.nan, burst_rate=np.nan)
-    frozen = float((speed < frozen_threshold).mean())
-    above = (speed >= frozen_threshold).astype(int)
-    bursts = int(np.sum(np.diff(above) == 1))
-    minutes = len(speed) / fs / 60.0
-    burst_rate = bursts / minutes if minutes > 0 else np.nan
-    return dict(frozen_fraction=frozen, burst_rate=burst_rate)
-
-
-# ---------------------------------------------------------------------------
-# Direction / extent
-# ---------------------------------------------------------------------------
-def sway_orientation(xy):
-    """
-    Principal sway-axis orientation and anisotropy of a 2-D point cloud.
-
-    A PCA of the mean-removed horizontal positions gives the orientation of
-    the major axis as an axial angle in ``[0, 180)`` degrees (undirected --
-    a line, not an arrow) and the anisotropy ``sqrt(lambda_max / lambda_min)``
-    of the sway ellipse. Anisotropy 1.0 is isotropic/circular; values above
-    ~1.3 indicate clearly directional sway.
-
-    Source: still standing study (Jensenius), sway-direction analysis.
-
-    Args:
-        xy (np.ndarray): Point cloud of shape ``(T, 2)``.
-
-    Returns:
-        dict: ``{"angle_deg", "anisotropy"}``. Both are ``nan`` if the
-            covariance is degenerate or there are too few finite points.
-    """
-    xy = np.asarray(xy, dtype=float)
-    xy = xy[np.isfinite(xy).all(axis=1)]
-    if len(xy) < 3:
-        return dict(angle_deg=np.nan, anisotropy=np.nan)
-    x = xy[:, 0] - xy[:, 0].mean()
-    y = xy[:, 1] - xy[:, 1].mean()
-    cov = np.cov(x, y)
-    w, v = np.linalg.eigh(cov)  # ascending eigenvalues
-    if w[0] <= 0:
-        return dict(angle_deg=np.nan, anisotropy=np.nan)
-    angle = float(np.degrees(np.arctan2(v[1, 1], v[0, 1])) % 180.0)
-    anisotropy = float(np.sqrt(w[1] / w[0]))
-    return dict(angle_deg=angle, anisotropy=anisotropy)
-
-
-def axial_rayleigh(angles_deg):
-    """
-    Axial Rayleigh test for a preferred orientation among axial angles.
-
-    Tests whether a sample of axial (undirected, ``[0, 180)`` deg) angles --
-    e.g. per-session principal sway axes -- clusters around a common
-    orientation. Angles are doubled to map the axial circle onto the full
-    circle before computing the mean resultant length ``R`` and the Rayleigh
-    p-value (small ``p`` with large ``R`` means a shared preferred axis).
-
-    Source: still standing study (Jensenius), sway-direction analysis.
-
-    Args:
-        angles_deg (np.ndarray): Axial angles in degrees.
-
-    Returns:
-        dict: ``{"R", "p", "mean_axis_deg", "n"}``.
-    """
-    a = 2 * np.radians(np.asarray(angles_deg, dtype=float))
-    a = a[np.isfinite(a)]
-    n = len(a)
-    if n < 2:
-        return dict(R=np.nan, p=np.nan, mean_axis_deg=np.nan, n=n)
-    C = np.mean(np.cos(a))
-    S = np.mean(np.sin(a))
-    R = float(np.hypot(C, S))
-    Z = n * R * R
-    p = float(np.exp(-Z) * (1 + (2 * Z - Z * Z) / (4 * n)))
-    mean_axis = float(np.degrees(0.5 * np.arctan2(S, C)) % 180.0)
-    return dict(R=R, p=p, mean_axis_deg=mean_axis, n=n)
-
-
-def spatial_extent(pos, fs, *, ellipse_conf=0.95, window_s=20.0,
-                   vertical_axis=None):
-    """
-    Spatial extent / occupied volume of a 3-D (or 2-D) position trace.
-
-    Complements sway magnitude (QoM) by describing *how large a region* a
-    marker occupies. Reports the RMS dispersion radius about the session
-    centroid, the Gaussian confidence-ellipsoid volume and its cube-root
-    radius, the mean within-window dispersion (which removes slow drift),
-    and a drift ratio ``full_dispersion / within_window_dispersion`` (``> 1``
-    when slow drift enlarges the occupied region over the session). When a
-    ``vertical_axis`` is given the drift is additionally split into
-    horizontal and vertical components.
-
-    Source: still standing study (Jensenius), spatial-range analysis
-    (``session_metrics``).
-
-    Args:
-        pos (np.ndarray): Position trace of shape ``(T, D)`` with ``D`` 2 or
-            3, in position units (e.g. mm).
-        fs (float): Sampling rate in Hz.
-        ellipse_conf (float, optional): Confidence level for the ellipsoid
-            volume. Defaults to 0.95.
-        window_s (float, optional): Window length in seconds for the
-            within-window dispersion. Defaults to 20.0.
-        vertical_axis (int, optional): Index of the vertical axis (``0``,
-            ``1`` or ``2``); when given, drift is decomposed into horizontal
-            and vertical parts. Defaults to None.
-
-    Returns:
-        dict: ``dispersion``, ``ellipsoid_volume``, ``ellipsoid_radius``,
-            ``within_window_dispersion``, ``drift_ratio`` and (when
-            ``vertical_axis`` is set) ``drift_horizontal`` and
-            ``drift_vertical``. Returns ``None`` if fewer than one window of
-            finite samples is available.
-    """
-    from scipy.stats import chi2
-
-    pos = np.asarray(pos, dtype=float)
-    if pos.ndim != 2:
-        raise ValueError("pos must have shape (T, D)")
-    D = pos.shape[1]
-    mask = np.isfinite(pos).all(axis=1)
-    P = pos[mask]
-    w = int(fs * window_s)
-    if len(P) < max(w, D + 2):
-        return None
-
-    centroid = P.mean(axis=0)
-    disp_vec = P - centroid
-    d = np.sqrt((disp_vec * disp_vec).sum(axis=1))
-    dispersion = float(np.sqrt((d * d).mean()))
-
-    cov = np.cov(P.T)
-    det = np.linalg.det(cov)
-    chi_d = chi2.ppf(ellipse_conf, D)
-    # volume of a D-dim Gaussian confidence ellipsoid
-    if D == 2:
-        vol = np.pi * chi_d * np.sqrt(max(det, 0.0))
-    else:  # D == 3
-        vol = (4.0 / 3.0) * np.pi * (chi_d ** 1.5) * np.sqrt(max(det, 0.0))
-    ellipsoid_volume = float(vol)
-    ellipsoid_radius = float(vol ** (1.0 / D))
-
-    def _within_window(Q):
-        cols = Q.shape[1] if Q.ndim == 2 else 1
-        Q = Q.reshape(len(Q), cols)
-        wr = []
-        for s in range(0, len(Q) - w + 1, w):
-            seg = Q[s:s + w]
-            cc = seg.mean(axis=0)
-            e = seg - cc
-            wr.append(np.sqrt((e * e).sum(axis=1).mean()))
-        return np.mean(wr) if wr else np.nan
-
-    def _full_dispersion(Q):
-        cols = Q.shape[1] if Q.ndim == 2 else 1
-        Q = Q.reshape(len(Q), cols)
-        cc = Q.mean(axis=0)
-        e = Q - cc
-        return np.sqrt((e * e).sum(axis=1).mean())
-
-    win_disp = float(_within_window(P))
-    drift_ratio = (dispersion / win_disp
-                   if win_disp and win_disp > 0 else np.nan)
-
-    out = dict(dispersion=dispersion, ellipsoid_volume=ellipsoid_volume,
-               ellipsoid_radius=ellipsoid_radius,
-               within_window_dispersion=win_disp, drift_ratio=drift_ratio)
-
-    if vertical_axis is not None:
-        hax = [k for k in range(D) if k != vertical_axis]
-
-        def _comp_drift(cols):
-            Q = P[:, cols]
-            wm = _within_window(Q)
-            full = _full_dispersion(Q)
-            return full / wm if wm and wm > 0 else np.nan
-
-        out["drift_horizontal"] = float(_comp_drift(hax))
-        out["drift_vertical"] = float(_comp_drift([vertical_axis]))
-
-    return out
+from micromotion.balance import (  # noqa: F401
+    axial_rayleigh,
+    confidence_ellipse_area,
+    convex_hull_area,
+    cop_sway_metrics,
+    dfa,
+    principal_axis_projection,
+    sample_entropy,
+    spatial_extent,
+    spectral_edges,
+    stabilogram_diffusion,
+    sway_orientation,
+    sway_texture,
+)
+
+__all__ = [
+    "axial_rayleigh",
+    "confidence_ellipse_area",
+    "convex_hull_area",
+    "cop_sway_metrics",
+    "dfa",
+    "principal_axis_projection",
+    "sample_entropy",
+    "spatial_extent",
+    "spectral_edges",
+    "stabilogram_diffusion",
+    "sway_orientation",
+    "sway_texture",
+]
diff --git a/musicalgestures/_qom.py b/musicalgestures/_qom.py
index a2f29ba..3e7f2d1 100644
--- a/musicalgestures/_qom.py
+++ b/musicalgestures/_qom.py
@@ -1,384 +1,38 @@
-"""
-Quantity-of-motion (QoM) signal cores for position, pose and accelerometer
-data.
+"""Re-exported from the ``micromotion`` package.
 
-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 used to live here. They were moved to ``micromotion`` on 2026-07-29 so that
+one implementation of quantity of motion exists rather than two, and MGT now depends on that
+package instead of carrying its own copy. Behaviour is unchanged: this module's tests pass
+against ``micromotion`` unmodified.
 
-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).
+The dependency points this way round on purpose. ``micromotion`` needs only numpy, scipy and
+pandas, so someone analysing accelerometer data does not have to install a computer-vision
+stack; MGT already depends on ``ambiscape`` the same way, and neither of those packages
+imports MGT.
 
-Sources: stillstanding study and Westney-comparisons study (Jensenius).
+Import from ``micromotion`` directly in new code.
 """
 
-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) and `fs_out` is its
-            sampling rate (equal to `fs` unless decimated). `speed` is empty
-            (and `fs_out` equals the input `fs`) when the input has fewer
-            than `int(fs) + 5` samples, or when it still contains non-finite
-            samples after per-dimension interpolation (i.e. a dimension had
-            fewer than 3 finite samples to interpolate from). In the
-            auto-decimate regime, `speed` is also empty (with `fs_out` the
-            decimated rate) when decimation leaves fewer than ~30 samples --
-            too few for a stable SOS band-pass.
-
-    Raises:
-        ValueError: If the band is invalid, i.e. does not satisfy
-            `0 < lo < hi <= 0.45*fs` (after `hi` is clipped to 0.9 x Nyquist).
-    """
-    from scipy import signal
-    pos = np.asarray(pos, float)
-    if pos.ndim == 1:
-        pos = pos[:, None]
-    pos = np.column_stack([_interp_nans(pos[:, i]) for i in range(pos.shape[1])])
-
-    hi_eff = min(hi, 0.9 * fs / 2)
-    if not (0 < lo < hi_eff <= 0.45 * fs + 1e-9):
-        raise ValueError("band must satisfy 0 < lo < hi <= 0.45*fs")
-
-    if len(pos) < int(fs) + 5 or not np.isfinite(pos).all():
-        return np.array([]), fs
-
-    if auto_decimate and fs / hi_eff >= 40:
-        q = int(min(13, fs // (20 * hi_eff)))
-        pos = signal.decimate(pos, q, axis=0, zero_phase=True)
-        fs_out = fs / q
-        if len(pos) < 30:
-            return np.array([]), fs_out
-        sos = signal.butter(2, [lo / (fs_out / 2), hi_eff / (fs_out / 2)],
-                            btype="band", output="sos")
-        filtered = signal.sosfiltfilt(sos, pos, axis=0)
-    else:
-        fs_out = fs
-        b, a = signal.butter(order, [lo / (fs / 2), hi_eff / (fs / 2)], btype="band")
-        filtered = signal.filtfilt(b, a, pos, axis=0)
-    speed = np.linalg.norm(np.diff(filtered, axis=0), axis=1) * fs_out
-    return speed, fs_out
-
-
-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. 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
-
-
-def grid_qom(frames, grid=(6, 4), region=(0.0, 1.0, 0.0, 1.0), threshold=8.0):
-    """
-    Spatial grid quantity of motion from a stack of grayscale frames: the
-    absolute inter-frame difference is thresholded (small differences set to
-    zero to suppress sensor noise) and averaged within each cell of a
-    `grid[0]` x `grid[1]` grid laid over `region`, yielding one motion time
-    series per cell plus a per-cell mean-motion heatmap.
-
-    Source: Westney-comparisons study (Jensenius) -- audience-region motion
-    mapping in a concert hall.
-
-    Args:
-        frames (np.ndarray): Grayscale frames of shape (T, H, W).
-        grid (tuple, optional): Grid size (columns, rows). Defaults to (6, 4).
-        region (tuple, optional): Region of interest as fractions
-            (x0, x1, y0, y1) of the frame. Defaults to the full frame.
-        threshold (float, optional): Absolute-difference threshold below which
-            pixel changes are zeroed (0-255 scale). Defaults to 8.0.
-
-    Returns:
-        tuple: `(series, heat)` where `series` has shape (T-1, rows*cols)
-            (cells in row-major order) and `heat` has shape (rows, cols) with
-            each cell's time-mean motion.
-
-    Raises:
-        ValueError: If `frames` is not 3-D (T, H, W), as in
-            `_motionanalysis.motiongram_data`.
-    """
-    frames = np.asarray(frames, dtype=np.float32)
-    if frames.ndim != 3:
-        raise ValueError("grid_qom expects frames of shape (T, H, W)")
-    T, H, W = frames.shape
-    gx, gy = grid
-    x0, x1, y0, y1 = region
-    xs = np.linspace(int(x0 * W), int(x1 * W), gx + 1).astype(int)
-    ys = np.linspace(int(y0 * H), int(y1 * H), gy + 1).astype(int)
-    d = np.abs(np.diff(frames, axis=0))
-    d[d < threshold] = 0.0
-    series = np.empty((T - 1, gy * gx), dtype=np.float32)
-    for r in range(gy):
-        for c in range(gx):
-            cell = d[:, ys[r]:ys[r + 1], xs[c]:xs[c + 1]]
-            series[:, r * gx + c] = cell.mean(axis=(1, 2))
-    heat = series.mean(axis=0).reshape(gy, gx)
-    return series, heat
+from micromotion.qom import (  # noqa: F401
+    accel_to_speed,
+    band_limited_qom,
+    bin_series,
+    body_scale,
+    envelope,
+    grid_qom,
+    group_qom,
+    normalized_qom,
+    pose_qom,
+)
+
+__all__ = [
+    "accel_to_speed",
+    "band_limited_qom",
+    "bin_series",
+    "body_scale",
+    "envelope",
+    "grid_qom",
+    "group_qom",
+    "normalized_qom",
+    "pose_qom",
+]
diff --git a/pyproject.toml b/pyproject.toml
index e70abbb..f9761f9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -24,6 +24,7 @@ classifiers = [
     "Programming Language :: Python :: 3.12",
 ]
 dependencies = [
+    "micromotion>=0.3",
     "numpy",
     "pandas",
     "matplotlib",