API reference¶
Every public function, grouped by the question being asked of a signal. The pages below are generated from the docstrings in the source, so a signature here is the signature the code has.
Quantity of motion¶
micromotion.qom ¶
Quantity of motion.
The signature measure of this research programme: the average speed of a body part, in millimetres per second, restricted to the micromotion band.
There is one definition and three named variants. The point of naming them is that they answer different questions and previously differed only by an undocumented line of filter code. Any figure or paper should be able to say which variant produced its numbers.
Definition¶
Band-limit each axis to 0.2-5 Hz, bring it to velocity, band-limit again, take the Euclidean norm across axes, and report the MEDIAN in mm/s.
The statistic is part of the definition rather than a presentation choice. The two are not
close: on accelerometer data the mean-to-median speed ratio is about 2, and one corpus record
carries a deposited mean of 12.79 mm/s beside a report quoting 11.12 for the same recordings
at the same band, the whole difference being this choice with neither document saying which it
had made. speed_from_position and speed_from_acceleration return the speed SERIES and
take no statistic, so the caller decides. Decide explicitly, and say which one a published
figure used.
Acceleration is brought to velocity by integration, position by differentiation. The second band-limiting is not cosmetic. Integrating a signal with any residual offset produces a ramp that dominates the result, and differentiating amplifies the high-frequency noise that the first pass was meant to exclude.
Variants¶
raw
The band as defined. Contains respiration, the ballistocardiac impulse and postural
sway together. This is what the deposited files report unless they say otherwise.
compensated
Respiration and cardiac activity removed: the lower edge is raised to 0.5 Hz and the
per-recording cardiac peak is notched out. What is left is postural micromotion.
tilt_corrected
For a single accelerometer only. A body-worn accelerometer cannot distinguish leaning
from translating: tilting into gravity produces an acceleration with no displacement.
Where a gyroscope is available the tilt component is estimated and removed. Measured
directly on the fNIRS session, tilt inflates raw QoM by 1.56x.
G
module-attribute
¶
G = 9.80665
Standard gravity, m/s^2.
Present because accelerometers export in g and the band-limited result must be in SI before integration. Getting this wrong is not hypothetical: every phone quantity of motion is then 9.80665x too large.
BANDS
module-attribute
¶
BANDS = {'micromotion': filters.BAND, 'wideband': filters.WIDEBAND, 'noresp': filters.NORESP_BAND, 'optical_legacy': filters.OPTICAL_LEGACY_BAND}
The four conventions in use, by name.
micromotion is 0.2-5 Hz and is the only one every device in the corpus can deliver, so
it is the one a cross-collection comparison must use.
wideband is 0.2-10 Hz, for jerk and other high-derivative measures that need the octave the
canonical band gives up. Only on collections sampled fast enough to reach it; check with
:func:effective_band, and do not infer the rate from a file's grid, which may be an upsample.
noresp is 0.45-5 Hz, the canonical band with its lower edge above respiration. On a
chest-worn sensor the 0.15-0.45 Hz stretch is dominated by respiratory chest tilt -- gravity
re-projected by the breathing ribcage, a rotation rather than a translation -- so the choice
between this and micromotion is a purpose decision: keep the respiratory term or exclude
it. See :data:~micromotion.filters.NORESP_BAND.
optical_legacy is the 10 Hz low-pass behind the published championship figures. It retains
sub-0.2 Hz postural drift and reads about 15 per cent higher. It is kept because those numbers
are in print, not because it is interchangeable with the others.
QomResult
dataclass
¶
Quantity of motion, with the series it was reduced from.
Source code in src/micromotion/qom.py
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 | |
binned ¶
binned(bin_s: float = 5.0)
Average speed in fixed-width bins, as a DataFrame.
The final bin is usually partial and is flagged rather than dropped, because a short bin is not comparable with a full one and silently including it inflated the deposited five-second series three- to fourteenfold.
Source code in src/micromotion/qom.py
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 | |
identify_acceleration_unit ¶
identify_acceleration_unit(acc, tol: float = 0.25) -> str
Which unit an accelerometer file is in, from its own values: "g", "mg" or "m/s^2".
A device that is mostly stationary is mostly measuring gravity, so the median vector norm of a resting recording sits near 1, 981 or 9.81 depending on the convention. Those are three orders of magnitude apart, and no plausible unit lies between them, so the identification is unambiguous whenever the recording is dominated by gravity.
THIS EXISTS BECAUSE NO FILE FORMAT IN THIS FIELD DECLARES ITS UNITS. Four accelerometers recorded simultaneously on one body in the Oslo corpus stored their values in three different conventions, none of them stated anywhere in the files. The method is Finn Upham's, from the analysis those recordings were made for: take the mean total acceleration while the participant lies still, and read off which constant it is near.
A UNIT ERROR IS THE HARDEST KIND TO NOTICE. It scales one recording by 9.8 or 981 and leaves every correlation, every rank statistic and every reliability estimate untouched, so nothing downstream complains. It shows up as a suspiciously ROUND ratio against a known value -- a real disagreement is ragged and a unit error is a constant.
tol is the fractional distance from a candidate at which the answer is still accepted.
Raises if the norm is near none of them, which means either the recording is not
gravity-dominated or the units are something this does not know about; in both cases guessing
would be worse than stopping.
identify_acceleration_unit(np.full((100, 3), [0.0, 0.0, 1.0])) 'g'
Source code in src/micromotion/qom.py
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 | |
velocity_from_acceleration ¶
velocity_from_acceleration(acc, fs: float, unit: str = 'm/s^2', lo: float = filters.BAND[0], hi: float = filters.BAND[1], notch_hz: float | None = None, integrate: str = 'rectangle') -> np.ndarray
Band-limited velocity, per axis, in mm/s, from acceleration.
acc is (n_samples, n_axes). unit is "m/s^2" or "g". Returns an array of
the same shape; :func:speed_from_acceleration is its Euclidean norm.
Use this when a descriptor needs the velocity vector rather than its magnitude -- jerk, spectral measures per axis, anything directional. Computing it from the speed alone is not equivalent, and integrating by hand invites a pipeline that differs from the rest of the corpus in the filter order or the quadrature rule.
integrate selects the quadrature rule, and the choice is not cosmetic. Both are in
common use and they differ by about 0.26 per cent on real phone data -- small, but a
systematic bias rather than noise, since the rectangle rule lags the signal by half a
sample.
Neither is universally right, so the default is "rectangle", which is what this
package's own reference numbers were computed with. Pass integrate="trapezoid" to
reproduce a pipeline that used the trapezoid rule; on one deposited value the two give
93.140 and 93.405 mm/s against a published 93.091.
Which rule the project should standardise on is an open question, deliberately not settled by this default.
Source code in src/micromotion/qom.py
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 | |
speed_from_acceleration ¶
speed_from_acceleration(acc, fs: float, unit: str = 'm/s^2', lo: float = filters.BAND[0], hi: float = filters.BAND[1], notch_hz: float | None = None, integrate: str = 'rectangle') -> np.ndarray
Band-limited speed, in mm/s, from acceleration.
acc is (n_samples, n_axes). unit is "m/s^2" or "g". This is the norm of
:func:velocity_from_acceleration; see that function for the integration options.
Source code in src/micromotion/qom.py
201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 | |
derivative ¶
derivative(x, fs: float) -> np.ndarray
Fourth-order central difference along the first axis.
The two-point central difference that numpy.gradient computes has the frequency
response sin(w*dt)/dt rather than w, so it increasingly under-reads towards Nyquist.
That does not matter at 200 Hz, where the top of the micromotion band is a twentieth of
Nyquist, but it matters at the 20 Hz common rate where the band edge is Nyquist itself:
on real optical data the two-point rule loses 4.9 per cent of the quantity of motion
across that resampling, and this rule loses 2.0 per cent.
A spectral derivative would be exact for a truly band-limited signal and is not used here, because these recordings do not begin and end at the same value and the implied wraparound step adds broadband energy that differentiation then amplifies. Measured on the same file, it inflated the result by 24 per cent.
Source code in src/micromotion/qom.py
219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 | |
velocity_from_position ¶
velocity_from_position(pos, fs: float, unit: str = 'mm', lo: float = filters.BAND[0], hi: float = filters.BAND[1]) -> np.ndarray
Band-limited velocity, per axis, in mm/s, from position.
pos is (n_samples, n_axes), normally the three coordinates of one optical marker.
unit is "mm" or "m". Returns an array of the same shape;
:func:speed_from_position is its Euclidean norm.
Pair with :func:velocity_from_acceleration when a descriptor must be computed the same
way across optical and accelerometer collections: take the velocity from whichever
function matches the recorded quantity, and everything downstream is identical.
Source code in src/micromotion/qom.py
243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 | |
speed_from_position ¶
speed_from_position(pos, fs: float, unit: str = 'mm', lo: float = filters.BAND[0], hi: float = filters.BAND[1]) -> np.ndarray
Band-limited speed, in mm/s, from position.
pos is (n_samples, n_axes), normally the three coordinates of one optical marker.
unit is "mm" or "m". This is the norm of :func:velocity_from_position.
Source code in src/micromotion/qom.py
269 270 271 272 273 274 275 276 277 278 279 280 281 | |
qom ¶
qom(data, fs: float, kind: str = 'acceleration', unit: str | None = None, variant: str = 'raw', band: str = 'micromotion', gyro=None, integrate: str = 'rectangle') -> QomResult
Quantity of motion for one recording.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
data
|
ndarray
|
(n_samples, n_axes) acceleration or position. |
required |
fs
|
float
|
Measured sampling rate. Use the rate measured from the timestamps,
not the nominal one; see :func: |
required |
kind
|
str
|
|
'acceleration'
|
unit
|
str
|
Defaults to |
None
|
variant
|
str
|
|
'raw'
|
band
|
str
|
|
'micromotion'
|
gyro
|
ndarray
|
(n_samples, 3) angular velocity in rad/s. Required by
|
None
|
integrate
|
str
|
|
'rectangle'
|
Returns:
| Name | Type | Description |
|---|---|---|
QomResult |
QomResult
|
The mean and median speed in mm/s, the full speed series, the rate it was
computed at, the variant, the notched cardiac frequency where one was found, and
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
Source code in src/micromotion/qom.py
311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 | |
remove_tilt ¶
remove_tilt(acc, gyro, fs: float, unit: str = 'm/s^2') -> np.ndarray
Subtract the gravity component that rotation moves between axes.
The sensor's orientation is tracked by integrating the gyroscope, the gravity vector is rotated into the sensor frame at each sample, and what remains is translation. The integration drifts, so the estimated gravity direction is high-passed back towards the measured one; this is a complementary filter, not an attitude estimator, and it is adequate only because the body barely moves.
Source code in src/micromotion/qom.py
403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 | |
tilt_fraction ¶
tilt_fraction(acc, gyro, fs: float, unit: str = 'm/s^2') -> dict
How much of a body-worn accelerometer's quantity of motion is tilt rather than travel.
An accelerometer cannot tell leaning from moving: rotating in the gravity field produces an acceleration with no displacement. Where a gyroscope is present the rotation is known, so the gravity component it accounts for can be removed and the two compared.
Measured on the fNIRS session this returns 1.56, meaning the raw figure is a little over half again the translational one. Do not read a single session as a population value; the point is that the inflation is measurable rather than assumed, and the assumption in circulation was 1.3 to 1.5.
Source code in src/micromotion/qom.py
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 | |
envelope ¶
envelope(x, fs, smooth=1.0, normalize=True)
Smooth, optionally z-scored envelope of a signal: Savitzky-Golay
smoothing (order 2, window smooth seconds) followed by
standardisation. Used to compare motion/audio envelopes across sources
on a common, amplitude-free scale.
Source: Westney-comparisons study (Jensenius).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
ndarray
|
Input 1-D signal. |
required |
fs
|
float
|
Sampling rate of the signal (Hz). |
required |
smooth
|
float
|
Smoothing window in seconds. None or 0 disables smoothing. Defaults to 1.0. |
1.0
|
normalize
|
bool
|
If True, z-score the result. Defaults to True. |
True
|
Returns:
| Type | Description |
|---|---|
|
np.ndarray: The smoothed (and optionally z-scored) envelope, same length as the input. |
Source code in src/micromotion/qom.py
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 | |
bin_series ¶
bin_series(x, fs, bin_s=1.0)
Mean of consecutive, non-overlapping bins of a signal (e.g. a per-second quantity-of-motion envelope from a per-frame speed series). Trailing samples that do not fill a whole bin are dropped.
Source: stillstanding study (Jensenius); also used in the Westney-comparisons study as a per-second envelope.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
ndarray
|
Input 1-D signal. |
required |
fs
|
float
|
Sampling rate of the signal (Hz). |
required |
bin_s
|
float
|
Bin length in seconds. Defaults to 1.0. |
1.0
|
Returns:
| Type | Description |
|---|---|
|
np.ndarray: One mean value per bin (empty if the signal is shorter than two bins). |
Source code in src/micromotion/qom.py
505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 | |
band_limited_qom ¶
band_limited_qom(pos, fs, lo=filters.BAND[0], hi=filters.BAND[1], 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.
.. note::
Defaults follow filters.BAND, like everything else in the package,
so a number cannot be quoted from a second band by accident.
:func:pose_qom pins its upper edge at 5 Hz because image-space
landmark jitter dominates above that. Since the band became 0.2-5 Hz
the two coincide, so that pin currently changes nothing; it is kept
because it is a statement about pose data rather than about the band.
Use :func:speed_from_position for new work, which band-limits
again after differentiating; this one does not, and reads high as a
result. See the interop guide for the comparison against
musicalgestures, which carries the same function.
.. warning::
The name overstates what this does, and it reads high as a result.
The position is band-limited; the speed derived from it is not.
Differentiation amplifies the high end, so the velocity carries energy
above hi that the stated band excludes, and none of it is removed.
Measured against :func:speed_from_position, which band-limits again
after differentiating, on a 200 Hz optical recording at a matched
band: this returns 3.1475 mm/s against 2.9754, i.e. 5.5 per
cent high. The decomposition is one-sided -- the differentiation rule
(first difference here, central difference there) accounts for 0.05 per
cent, and the missing second band-pass for the remaining 5.5. At the
respective defaults the gap is larger still, because this function
band-limits only the position, not the speed derived from it.
Prefer :func:qom or :func:speed_from_position for new work. This is
kept, unchanged, so that figures computed with musicalgestures
continue to reproduce -- not because it is the better measure. Whichever
you use, say which.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
pos
|
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. |
required |
fs
|
float
|
Sampling rate of the trajectory (Hz). |
required |
lo
|
float
|
Lower band edge (Hz). Defaults to |
BAND[0]
|
hi
|
float
|
Upper band edge (Hz), clipped to 0.9 x Nyquist.
Defaults to |
BAND[1]
|
order
|
int
|
Butterworth order of the direct band-pass. Defaults to 4. |
4
|
auto_decimate
|
bool
|
Enable the decimate+SOS low-band regime. Defaults to True. |
True
|
Returns:
| Name | Type | Description |
|---|---|---|
tuple |
|
lo=None gives a pure low-pass with no lower edge, which is the
:data:~micromotion.OPTICAL_LEGACY_BAND convention and what the pre-2020 optical
standstill studies used. It retains the slow postural drift the corpus band removes, so a
value computed that way is not comparable with one computed at BAND.
Raises:
| Type | Description |
|---|---|
ValueError
|
If the band is invalid, i.e. does not satisfy
|
Source code in src/micromotion/qom.py
531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 | |
accel_to_speed ¶
accel_to_speed(acc, fs, highpass=filters.BAND[0], 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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
acc
|
ndarray
|
Acceleration of shape (N, 3) in m/s^2 (or raw counts
with |
required |
fs
|
float
|
Sampling rate (Hz). |
required |
highpass
|
float
|
High-pass cutoff (Hz) used both before and
after integration. Defaults to |
BAND[0]
|
order
|
int
|
Butterworth order of the high-pass filters. Defaults to 2. |
2
|
normalize_gravity
|
bool
|
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. |
False
|
Returns:
| Type | Description |
|---|---|
|
np.ndarray: Speed series of length N (m/s). |
Source code in src/micromotion/qom.py
665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 | |
group_qom ¶
group_qom(points, fs, lo=filters.BAND[0], hi=filters.BAND[1], normalize='visible', **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.
.. warning::
normalize decides what the divisor is, and the number this returns
moves with it. Say which was used.
normalize="visible", the default, excludes each marker at the frames
where it was absent and averages over the rest. On twelve markers with a
realistic dropout pattern, a median of eight visible, it lands within
0.8 per cent of the unoccluded truth.
normalize="worn" averages over every marker at every frame instead.
Since band_limited_qom interpolates gaps, an occluded marker then
contributes a near-zero speed while still counting in the divisor, and
the result tracks how much the cameras saw: on the same data it reads
16 to 17 per cent low and its speed series correlates +0.25 to +0.70
with the per-frame count of visible markers. It is kept so that a figure
computed that way keeps reproducing; it is bit-for-bit identical on one
machine and agrees to about 1 part in 10^7 across platforms, since
filtfilt is not bit-reproducible between scipy builds.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
points
|
ndarray
|
Trajectories of shape (N, M, D): N frames, M markers/landmarks, D spatial dimensions. |
required |
fs
|
float
|
Sampling rate (Hz). |
required |
lo
|
float
|
Lower band edge (Hz). Defaults to |
BAND[0]
|
hi
|
float
|
Upper band edge (Hz). Defaults to |
BAND[1]
|
normalize
|
str
|
|
'visible'
|
**kwargs
|
Passed on to |
{}
|
Returns:
| Name | Type | Description |
|---|---|---|
tuple |
|
Source code in src/micromotion/qom.py
731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 | |
pose_qom ¶
pose_qom(landmarks, fs, lo=filters.BAND[0], 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.2-5 Hz), where higher bands are dominated by landmark
jitter rather than motion.
Source: Westney-comparisons study (Jensenius).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
landmarks
|
ndarray
|
Landmark trajectories of shape (N, L, 2) in pixels (a single landmark of shape (N, 2) is also accepted). |
required |
fs
|
float
|
Sampling rate (Hz, e.g. video frame rate). |
required |
lo
|
float
|
Lower band edge (Hz). Defaults to |
BAND[0]
|
hi
|
float
|
Upper band edge (Hz). Defaults to 5.0. |
5.0
|
**kwargs
|
Passed on to |
{}
|
Returns:
| Name | Type | Description |
|---|---|---|
tuple |
|
Source code in src/micromotion/qom.py
814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 | |
body_scale ¶
body_scale(landmarks, upper=(11, 12), lower=(23, 24))
Body-size scale (in the landmarks' own units, e.g. pixels) as the median
torso length: the distance from the midpoint of the upper landmarks
(shoulders) to the midpoint of the lower landmarks (hips). The torso
length is preferred over shoulder width because it stays robust in a
profile view, where the shoulder width collapses.
The default indices are MediaPipe Pose landmarks (11/12 shoulders, 23/24 hips).
Source: Westney-comparisons study (Jensenius).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
landmarks
|
ndarray
|
Landmark trajectories of shape (N, L, C) with C >= 2; only the first two coordinates are used. |
required |
upper
|
tuple
|
Indices of the two shoulder landmarks. Defaults to (11, 12). |
(11, 12)
|
lower
|
tuple
|
Indices of the two hip landmarks. Defaults to (23, 24). |
(23, 24)
|
Returns:
| Name | Type | Description |
|---|---|---|
float |
Median torso length (NaN if no finite frames). |
Source code in src/micromotion/qom.py
840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 | |
normalized_qom ¶
normalized_qom(landmarks, fs, scale=None, lo=filters.BAND[0], 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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
landmarks
|
ndarray
|
Landmark trajectories of shape (N, L, 2). |
required |
fs
|
float
|
Sampling rate (Hz). |
required |
scale
|
float
|
Precomputed body scale. Defaults to None (which
computes |
None
|
lo
|
float
|
Lower band edge (Hz). Defaults to |
BAND[0]
|
hi
|
float
|
Upper band edge (Hz). Defaults to 5.0. |
5.0
|
upper
|
tuple
|
Shoulder landmark indices for |
(11, 12)
|
lower
|
tuple
|
Hip landmark indices for |
(23, 24)
|
**kwargs
|
Passed on to |
{}
|
Returns:
| Name | Type | Description |
|---|---|---|
tuple |
|
Source code in src/micromotion/qom.py
872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 | |
grid_qom ¶
grid_qom(frames, grid=(6, 4), region=(0.0, 1.0, 0.0, 1.0), threshold=8.0)
Spatial grid quantity of motion from a stack of grayscale frames: the
absolute inter-frame difference is thresholded (small differences set to
zero to suppress sensor noise) and averaged within each cell of a
grid[0] x grid[1] grid laid over region, yielding one motion time
series per cell plus a per-cell mean-motion heatmap.
Source: Westney-comparisons study (Jensenius) -- audience-region motion mapping in a concert hall.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
frames
|
ndarray
|
Grayscale frames of shape (T, H, W). |
required |
grid
|
tuple
|
Grid size (columns, rows). Defaults to (6, 4). |
(6, 4)
|
region
|
tuple
|
Region of interest as fractions (x0, x1, y0, y1) of the frame. Defaults to the full frame. |
(0.0, 1.0, 0.0, 1.0)
|
threshold
|
float
|
Absolute-difference threshold below which pixel changes are zeroed (0-255 scale). Defaults to 8.0. |
8.0
|
Returns:
| Name | Type | Description |
|---|---|---|
tuple |
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
Source code in src/micromotion/qom.py
914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 | |
Filters and bands¶
micromotion.filters ¶
Band-limiting for the micromotion band.
One definition, used everywhere. Across the Still Standing repository 46 scripts defined their own filter and they did not all agree; the differences were invisible in the output and moved quantity of motion by up to 10 per cent.
The canonical band is 0.2-5 Hz, a zero-phase Butterworth of order 4 applied as second-order sections. The lower edge sits below the respiratory rate and above the postural drift that integration turns into a ramp; the upper edge is set by what the slowest device in the corpus can actually deliver.
BAND
module-attribute
¶
BAND = (0.2, 5.0)
The micromotion band, in Hz.
The lower edge was chosen by a sweep across seven optical datasets and 665 recordings: the between-dataset spread is 3.2 per cent at 0.15 Hz, 2.1 at 0.20, 2.7 at 0.25, 6.2 at 0.30 and 10.1 at 0.40. 0.2 Hz is a clear optimum, and it is also where the 20 Hz origin dataset stops being an outlier in either direction.
An edge this close to DC has to be checked against integration drift, since that is what it is there to control. It survives: on accelerometer data the ratio of mean to median speed, which rises if drift is leaking in, is 2.07 at 0.2 Hz against 2.00 at 0.3 -- flat.
The upper edge is set by deliverability, not by taste. A band above Nyquist is not a convention but a defect that returns a plausible number, and the ceiling must therefore be one that every device in the corpus can support:
============================ ================== ========== ========== collection sampling Nyquist 5 Hz? ============================ ================== ========== ========== phone, fused linear accel ~15 Hz 7.5 yes optical, 20 Hz subset 20 Hz 10 yes collaborators' audience data 10 Hz 5 at Nyquist optical and inertial, rest 100-256 Hz 50-128 yes ============================ ================== ========== ==========
A 10 Hz ceiling fails the first row on 354 of 355 days and sits exactly on Nyquist for the second. 5 Hz clears every one of them.
The first row is worth stating precisely, because the obvious reading is wrong. The phone's accelerometer runs at about 50 Hz. What runs at 15 Hz is the linear acceleration, which the logging app derives by fusing the accelerometer with the gyroscope and magnetometer to remove gravity -- and a fusion cannot output faster than its slowest input, which is the 15 Hz gyroscope. The deposited files carry that fused channel, so their usable Nyquist really is 7.5 Hz. The limit is the channel chosen, not the hardware.
What it costs. Band-limited speed mostly does not notice: measured over 466 person-recordings, 5 Hz against 10 differs by a median of 1.3 per cent, and 95 per cent of quiet-standing sway power lies below 1 Hz anyway. That is a distribution rather than a bound. Its ninetieth percentile is 3.2 per cent and its maximum 9.1, largest on optical collections sampled at 100 to 200 Hz, which genuinely resolve the octave between 5 and 10 Hz. Rankings survive: on one 199-recording optical collection the change moves the median 0.8 per cent and leaves the ranking at Spearman 0.996.
What it costs that matters. Jerk is two derivatives higher and lives in the discarded
region: at 5 Hz it is 37 to 66 per cent of its 10 Hz value, and the ranking shifts too. So jerk
must not be computed at this band on data that could support a wider one -- use
:data:WIDEBAND and say so. On the phone collection the wider jerk was never real: the 15 Hz
sensor cannot produce it, and computing it there inflated jerk 18 to 27 per cent with
interpolation.
WIDEBAND
module-attribute
¶
WIDEBAND = (0.2, 10.0)
0.2-10 Hz: the band for measures that need the octave :data:BAND gives up.
Jerk and other high-derivative quantities live between 5 and 10 Hz, where band-limited speed
does not. This band is for them, on collections sampled fast enough to deliver it -- check
with :func:effective_band first, and never assume it from a file's grid, which may be an
upsample of a much slower sensor.
It is deliberately not the default. A quantity computed here is not comparable with one
computed at :data:BAND, and it is not computable at all on the slower collections.
NORESP_BAND
module-attribute
¶
NORESP_BAND = (0.45, 5.0)
0.45-5 Hz: the micromotion band with its lower edge raised above respiration.
The 0.15-0.45 Hz stretch that :data:BAND keeps is, on a chest-worn sensor, dominated by
respiratory chest tilt. The ribcage turns the sensor as it expands, so what the accelerometer
reads there is gravity re-projected between axes -- a rotation, not the sensor travelling.
Measured across the year-long chest-phone record, 94 to 99 per cent of the power in that band
lies perpendicular to that day's gravity vector, implying a tilt of 0.05 to 0.14 degrees. This
band excludes it; a quantity computed here is micromotion above respiration.
The choice between :data:BAND and this one is a purpose decision, not a correctness one.
:data:BAND keeps the respiratory term, which makes it the cardiorespiratory torso measure and
the band every cross-collection comparison must use; this one drops the term, which makes it the
band for a postural question on a chest-worn sensor. The StillStanding365 record deposits
quantity of motion at both, as qom_mm_s and qom_045_5hz_mm_s, each with its band stated
-- which is the practice to copy, because a value at either band with the band unstated will
sooner or later be compared against one at the other.
It is not the compensated variant of :func:~micromotion.qom.qom, which raises the edge
further, to 0.5 Hz, and also notches the recording's own cardiac peak. This constant changes
the band and nothing else: the ballistocardiac impulse is still inside it.
OPTICAL_LEGACY_BAND
module-attribute
¶
OPTICAL_LEGACY_BAND = (None, 10.0)
A 10 Hz low-pass with no lower edge: the convention behind the published championship quantity of motion.
It is kept because those numbers are in print, not because it is interchangeable with
:data:BAND. Optical position is an absolute measurement, so sub-0.2 Hz postural drift in
it is real movement and there is no reason to discard it. A body-worn accelerometer cannot
offer the same choice: gravity is a DC term, and integrating any residual offset produces a
ramp that swamps the result. So the lower edge is optional for position and mandatory for
acceleration.
The two are not the same measure. On the 2015 championship the band-pass reads 15.5 per
cent below the low-pass. Any table that puts optical and accelerometer collections side by
side must therefore use :data:BAND throughout, and say so.
NARROW_BAND_RATIO
module-attribute
¶
NARROW_BAND_RATIO = 1000.0
Warn when the sampling rate exceeds the lower band edge by more than this.
The lower edge is what drives the conditioning, which is why the ratio is taken against it. The worked example below fails at 0.15 Hz, near the bottom of its band, not near the top; a test keyed on the upper edge measures the wrong thing and moves whenever the ceiling moves. It did: halving the canonical ceiling from 10 Hz to 5 doubled an upper-edge ratio while the conditioning was unchanged, and the warning began firing on every high-rate call.
A band-pass whose edges sit very close to zero in normalised frequency is numerically fragile, and how fragile depends on how it is realised. Second-order sections, which this module uses, stay accurate far longer than the transfer-function form -- but not forever, and a caller designing their own filter should be warned before they are bitten.
The failure is silent and it is not hypothetical. A 0.1-0.5 Hz third-order band-pass at
250 Hz, written in the usual butter(3, [lo/ny, hi/ny]) transfer-function form, has a
largest pole radius of 0.9979 and a measured passband gain of 0.84 at 0.15 Hz where it should
be 0.99. Nothing raises, nothing looks wrong, and every amplitude downstream is 16 per cent
low. The same design as second-order sections gives 0.9875.
The fix when this warns is to decimate first, so the band sits comfortably inside the new Nyquist range, then filter.
NYQUIST_MARGIN
module-attribute
¶
NYQUIST_MARGIN = 0.99
How close to Nyquist an upper band edge may sit, as a fraction.
A filter designed right at Nyquist has no transition band left, so the edge is pulled in.
The default is conservative. The Still Standing corpus uses 0.999 throughout, which matters
whenever the band edge is already near Nyquist -- a 10 Hz low-pass on 20 Hz data becomes
9.9 Hz here and 9.99 Hz there, and on real optical data that moved quantity of motion by
7e-5. Pass margin to match whatever convention the surrounding analysis uses.
bandpass ¶
bandpass(x, fs: float, lo: float | None = BAND[0], hi: float = BAND[1], order: int = ORDER, margin: float = NYQUIST_MARGIN)
Zero-phase band-limiting along the first axis.
lo=None gives a pure low-pass, which is the :data:OPTICAL_LEGACY_BAND convention.
Source code in src/micromotion/filters.py
220 221 222 223 224 225 226 227 228 229 230 | |
lowpass ¶
lowpass(x, fs: float, fc: float = BAND[1], order: int = ORDER, margin: float = NYQUIST_MARGIN)
Zero-phase low-pass along the first axis.
margin is how close to Nyquist fc may sit; see :data:NYQUIST_MARGIN.
Source code in src/micromotion/filters.py
233 234 235 236 237 238 239 240 241 242 | |
highpass ¶
highpass(x, fs: float, fc: float = BAND[0], order: int = ORDER)
Zero-phase high-pass along the first axis.
Used where the upper edge is meaningless because the rate is already near the band limit, and for gravity removal when no low-pass is wanted.
Source code in src/micromotion/filters.py
245 246 247 248 249 250 251 252 253 254 255 | |
notch ¶
notch(x, fs: float, f0: float, q: float = 6.0)
Zero-phase notch at f0 Hz.
Used to remove the cardiac peak when isolating postural micromotion. f0 is
normally found with :func:micromotion.spectral.cardiac_peak.
Source code in src/micromotion/filters.py
258 259 260 261 262 263 264 265 266 267 | |
edge_transient_samples ¶
edge_transient_samples(fs: float, lo: float | None = BAND[0], order: int = ORDER) -> int
Samples at each end that filtfilt contaminates.
A conservative estimate: the impulse response of the low edge, doubled for the
forward-backward pass. Callers should either trim this or flag it, as the deposited
five-second binning does with its edge column.
Source code in src/micromotion/filters.py
270 271 272 273 274 275 276 277 | |
effective_band ¶
effective_band(fs: float, lo: float | None = BAND[0], hi: float = BAND[1], margin: float = NYQUIST_MARGIN) -> tuple[float | None, float]
The band that will actually be applied at this sampling rate.
The requested upper edge is clamped to just below Nyquist, so a rate below twice hi
silently narrows the band. Ask before comparing two results computed at different rates:
at 8 Hz the canonical 0.2-5 Hz band becomes 0.2-3.96, which is a different measurement.
effective_band(100.0) (0.2, 5.0) effective_band(10.0) (0.2, 4.95)
Source code in src/micromotion/filters.py
280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 | |
Sampling rates and resampling¶
micromotion.resample ¶
Rate measurement and resampling.
Two rules, both learned the hard way.
Measure the rate, do not read it. Nominal rates in this corpus are wrong by up to 4.4 per cent, and one record's documented rate was out by a factor of 37.
Downsample, never upsample. Upsampling invents structure between samples, and every method that reads across scales treats the invention as real. An analysis that upsampled 20 Hz data to 25 Hz produced multifractal widths up to 6.6 where the plausible range is around 1. Nothing failed; the numbers were simply wrong.
COMMON_RATE
module-attribute
¶
COMMON_RATE = 20.0
Hz. The rate at which every collection can be compared, including the slowest.
It is the greatest common divisor of the corpus's optical rates -- 20, 100, 120 and 200 Hz -- so every recording reaches it by an integer decimation and none requires upsampling. That is its one virtue, and it is a real one: it is the only rate at which the natively-20 Hz origin study can be placed beside the rest at all.
It is a lossy rate, not a free one. The argument that the band stops at 10 Hz so 20 Hz discards nothing does not survive measurement, because a 10 Hz upper edge cannot be realised at 20 Hz: Nyquist sits exactly on it, the margin rule pulls the edge inside, and the anti-alias filter is already rolling off below it. Decimating 34 natively-200 Hz person-recordings and re-measuring moves quantity of motion by -2.09 per cent at the median, and between -0.30 and -10.57 per cent across recordings. A per-recording spread that wide is a distortion rather than a bias: it cannot be measured once and corrected away.
Prefer :data:HARMONISED_RATE where every series in the comparison can reach it. Use this one
when the comparison must include a 20 Hz collection, and say in the output that it is the lossy
view.
HARMONISED_RATE
module-attribute
¶
HARMONISED_RATE = 100.0
Hz. The preferred rate for a comparison whose series can all reach it.
Native for several collections, an exact halving from 200 Hz, and a 6-to-5 polyphase step from 120 Hz. Measured against native 200 Hz values it costs +0.02 per cent at the median and stays within +/-0.85 per cent, against -2.09 and up to -10.57 at 20 Hz. It also leaves the 0.2-5 Hz band comfortably inside Nyquist rather than sitting on it.
measured_rate ¶
measured_rate(t) -> float
Sampling rate in Hz, from a timestamp vector.
Sample count over elapsed span, deliberately not the reciprocal of the median interval. Where timestamps are rounded to whole milliseconds the intervals become a mixture of adjacent integers and their median is a quantisation artefact: on one dataset that route returns exactly 250 Hz for a recording that runs at 256, and on the pre-study phones it returned 636 Hz for a stream arriving at 106.
Source code in src/micromotion/resample.py
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 | |
rate_quality ¶
rate_quality(t) -> dict
How regular a timestamp vector actually is.
Returns the measured rate, the jitter, the largest gap, and counts of duplicated and backward timestamps. The balance-board files carry 4.7 per cent duplicates and 83 backward steps, and one pre-study file is 19 per cent covered because of a single 132-second gap; neither is visible from the rate alone.
Source code in src/micromotion/resample.py
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 | |
to_rate ¶
to_rate(x, fs_in: float, fs_out: float = COMMON_RATE)
Anti-alias resample to fs_out, refusing to upsample.
Raises rather than upsampling. If a series genuinely cannot reach the target it does not belong in a comparison built at that rate, and interpolating it in would corrupt the comparison silently.
Source code in src/micromotion/resample.py
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 | |
regularize ¶
regularize(t, x, fs_out: float | None = None, max_gap_s: float | None = None)
Put an irregularly sampled signal onto a uniform grid.
Sorts, drops duplicate and backward timestamps, then interpolates linearly. Samples
that fall inside a gap longer than max_gap_s are returned as NaN rather than
bridged, so that a 132-second hole cannot be mistaken for 132 seconds of stillness.
Written for the Wii balance board, whose timestamps are unsorted, 4.7 per cent duplicated, and arrive at a median 61.8 Hz that varies between files.
Source code in src/micromotion/resample.py
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 | |
interpolate_gaps ¶
interpolate_gaps(x, max_gap: int = 200)
Bridge short runs of NaN, leave long ones alone.
Filters cannot run across missing samples, so gaps have to be handled before anything else. Bridging a dropped frame is reconstruction; bridging a 469-second hole is invention, and the difference is only a matter of degree, which is why the threshold is explicit and the long gaps stay NaN for the caller to exclude.
Works column by column. Leading and trailing gaps are never filled, since there is nothing on one side to interpolate from.
Source code in src/micromotion/resample.py
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 | |
gap_report ¶
gap_report(x, fs: float) -> dict
Where the missing data is, and how it is distributed.
A single missing fraction hides the distinction that matters: one per cent scattered evenly is a usable recording, and one per cent in a single block in the middle is two recordings.
Source code in src/micromotion/resample.py
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 | |
Validation¶
micromotion.validate ¶
Checks that fail loudly on the errors this corpus makes silently.
Every check here exists because the failure it catches happened, went unnoticed, and produced a plausible number that someone used. None of them raised anything at the time. That is the common thread and the reason for the module: the corpus's characteristic failure is not a crash but a believable wrong answer, and the only defence is a check that runs on every build.
The intended use is a gate, not a report. Run :func:validate_series over everything a
harmonised table is about to be built from, and refuse to build if anything comes back at
"error". A finding at "warning" is something to record in the manifest beside the number
it affects.
Each check names the incident that motivated it, with its numbers, so that the tolerance is arguable rather than magic.
HUMAN_FLOOR_MM_S
module-attribute
¶
HUMAN_FLOOR_MM_S = 1.0
Below this band-limited speed, an optical trace is equipment rather than a person.
Calibrated on the Oslo Standstill Database rather than chosen: across 649 championship person-recordings the people span 2.326 to 17.705 mm/s and the 84 known tripod and reference-marker traces span 0.024 to 0.306, with an order of magnitude of empty space between them. 1.0 sits in that gap.
Finding
dataclass
¶
One thing wrong with one series or file.
Source code in src/micromotion/validate.py
29 30 31 32 33 34 35 36 37 38 39 40 | |
zero_triplets ¶
zero_triplets(x, where: str = '', max_fraction: float = 0.0) -> list[Finding]
Rows where every coordinate is exactly zero, which are gaps and not positions.
Qualisys writes a dropped frame as 0.000 0.000 0.000. That is a point on the laboratory
floor about a metre and a half below a standing head, so a reader that takes it literally
sees the head leave and return. A median-based measure barely notices; anything that sums or
integrates does. On one 2021 recording, 93 such frames out of 118698 gave the head a path
length of 119.8 m where the true figure is 11.0 m.
Exact zeros in all axes at once do not occur in real optical data, so the default tolerance
is zero. Raise max_fraction only for a sensor whose true output can sit at the origin.
Source code in src/micromotion/validate.py
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 | |
marker_average ¶
marker_average(markers, where: str = '', max_gap_fraction: float = 0.5) -> list[Finding]
Check a set of markers before averaging them into one position.
Averaging several markers into a single "head" or "trunk" position is routine and looks harmless. It is not, if the gaps have not been repaired first, because the usual repair happens at the end of a pipeline and an average destroys the evidence on the way in: the mean of two real coordinates and one zero triplet is a perfectly finite point that no later gap check will flag.
The damage is a clean multiplicative bias. Markers on one rigid segment move together, so
with n markers of which k are dead, the averaged position moves at about
(n - k) / n of the true amplitude -- and a speed derived from it is understated by the
same factor. One marker dead out of three is exactly two thirds, which reads as a third
less motion.
This happened. Four recordings in one 86-session collection carried a head marker that was never tracked, and their quantity of motion was reported 33.3 per cent low for as long as the collection existed, looking like unusually still standing rather than like a fault.
Pass a mapping of name to (n, 3) array. Repair each marker with
:func:micromotion.validate.zero_triplets and NaN before averaging, then use nanmean.
The same bias appears one level up in :func:micromotion.qom.group_qom, whose
normalize argument decides whether markers that were not visible count in the
divisor. Where positions are averaged here, check how speeds are averaged there.
Source code in src/micromotion/validate.py
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 | |
implausible_position ¶
implausible_position(x, where: str = '', axis: int = 2, min_fraction: float = 0.3, max_fraction: float = 2.5) -> list[Finding]
Samples that put a marker somewhere a body cannot be.
:func:zero_triplets catches a dropped frame written as three exact zeros. It cannot catch the
near miss: a reconstruction that lands close to the laboratory origin without being exactly
on it. Those samples pass every finiteness and sentinel test, because they are ordinary
numbers, and they are not rare enough to ignore -- a corpus of 1018 optical person-recordings
held two, one of which placed a head marker 139 mm below the floor.
The test is physical rather than statistical: a marker on a standing body stays within a band
around its own median height. Anything below min_fraction of that median, or above
max_fraction, is a tracking artefact rather than a posture.
The damage is uneven, which is why this is worth checking separately. A median-based measure barely notices 0.5 per cent of samples. Anything spatial is destroyed by them: on that Sverm recording the sway extent read 977 mm where the true figure is about 48.
Only meaningful for a marker whose median height is a real standing height, so recordings
whose median falls below 500 in the array's own units are skipped.
Source code in src/micromotion/validate.py
151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 | |
marker_noise ¶
marker_noise(x, fs: float, where: str = '', max_ratio: float = 5.0) -> list[Finding]
A marker that neither jumps nor drops out, but jitters.
:func:zero_triplets catches the dropped frame and :func:implausible_position catches the
reconstruction that lands near the laboratory origin. Neither can see the third failure, which
is a marker whose every sample is plausible and whose sample-to-sample noise is several times
what the body contributes. Nothing about such a trace looks wrong: it stays at head height, it
never leaps, and its median-based quantity of motion is perfectly ordinary.
It is destroyed only by measures that SUM. On one Sverm recording the band-limited quantity of motion is 4.95 mm/s -- the corpus median, an unremarkable standstill -- while the raw sample-to-sample path length runs at 79.18 mm/s, sixteen times higher. Plotted as cumulative distance beside 190 other recordings it was the obvious outlier, and it is not a person who moved.
The test compares the two. Raw path speed is the mean sample-to-sample displacement per second, which counts everything including the sensor's own jitter; band-limited speed keeps only the frequencies a standing body moves in. Their ratio is therefore how much of the measured path lies outside the band, and it is bounded below by 1 rather than by 0. Over 195 Sverm person-recordings it has a median of 1.39 and a 95th percentile of 2.37, then a gap to 4.5, 5.9, 10.8 and 16.2 -- so the default threshold sits in empty space rather than on a shoulder.
THE DENOMINATOR IS THE MEDIAN, and that was checked rather than assumed. Re-measured on
2026-08-15 both ways over the same recordings: against the median band-limited speed, which is
what this function uses, the ratio reproduces the figures above; against the MEAN it gives a
median of 1.22, a 95th percentile of 2.10 and a tail of 3.8, 5.2, 9.7 and 11.0, which is not
what is quoted here. So the calibration and the code agree about which statistic they are
dividing by. _curation/marker_jitter_sweep.py in the standstill corpus reports the MEAN-based
ratio in its own table, so its numbers are the second set and are not comparable with these.
Sampling rate matters and is already accounted for: a faster recording accumulates more raw path for identical behaviour, but it accumulates the same band-limited speed, so the ratio rises with rate. That is the point. It is asking how much of what a summing measure would count is not the body, and the answer legitimately depends on how often the sensor was asked.
Requires positions in millimetres. Returns nothing for a series too short to filter.
A GAP DOES NOT SILENCE IT, since 1.12.2. Until then a single non-finite sample anywhere in the
trace returned no finding at all, which is the failure this library exists to prevent: no
finding is what a clean recording returns, so one NaN made a jittering marker read as checked
and sound. The guard was there for a real reason -- a NaN propagates through diff into the
sum and through the filter into every sample of the band-limited series, so neither number
survives it -- but the answer is to measure what can be measured and say so, not to fall
silent. The check now runs on the longest contiguous finite run, reports the ratio for that
span, and says in the message that it is a span rather than the recording. Where no run is long
enough to filter it returns a WARNING saying the check could not run, so the gap is visible in
the same list as everything else.
This changes no verdict in the corpus it was written for: of 934 optical position recordings swept, 10 carry any gap at all and none of the 10 has a ratio anywhere near the threshold. It is a latent defect rather than one that has fired, and it is fixed because a check that cannot fail reads as coverage.
Source code in src/micromotion/validate.py
196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 | |
finite_fraction ¶
finite_fraction(x, where: str = '', min_finite: float = 0.8) -> list[Finding]
Whether enough of a series survived to measure, and whether any of it did.
A gap that runs off the start or end of a series cannot be interpolated — there is nothing on the far side — and a band-pass then spreads the surviving NaN across the whole recording. The result is an all-NaN series that is indistinguishable from an absent marker unless something looks. One Sverm 2012 recording lost a marker 548 s in and never regained it; the deposited value for it had been computed straight across the hole.
An emptied series is an error. A merely thin one is a warning, because the caller may legitimately be measuring the longest clean span instead.
Source code in src/micromotion/validate.py
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 | |
longest_finite_span ¶
longest_finite_span(x) -> tuple[int, int]
Start index and length of the longest run with no missing sample.
What to measure over when a gap cannot be bridged. Filtering across the gap is not an option and dropping the samples silently closes it, which is worse: the series then claims a duration it does not have.
Source code in src/micromotion/validate.py
313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 | |
timestamps ¶
timestamps(t, where: str = '') -> list[Finding]
Whether a timestamp column is usable as a clock.
Sorting a timestamp column into order is the tempting repair and the wrong one: it destroys the evidence that the clock misbehaved while leaving the samples in an order the sensor never produced. One balance-board collection carries 123 111 duplicate timestamps and 83 that step backwards, which is a device fault to be recorded, not a sort key to be fixed.
Source code in src/micromotion/validate.py
331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 | |
rate_agreement ¶
rate_agreement(t, documented_hz: float, where: str = '', tolerance: float = 0.02) -> list[Finding]
Whether the rate written down matches the rate the timestamps imply.
Measure the rate, do not read it. Documented rates in this corpus are wrong by up to 4.4 per cent, one record's was out by a factor of 37, and one championship's accelerometers turn out to run at 191.29–207.73 Hz against a nominal 200 — with the true rate a property of the individual device rather than the protocol, so three participants sharing a unit share a clock and everyone else does not.
A disagreement is an error rather than a warning because every frequency-domain measure downstream scales with it, and nothing further along can detect it.
Source code in src/micromotion/validate.py
362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 | |
held_samples ¶
held_samples(x, where: str = '', max_run: int = 50) -> list[Finding]
Long runs of an identical value, which mean a hold rather than a measurement.
A sensor sampled below the rate it is stored at is written out with each value repeated, and the file then claims a rate the data does not carry. The Delsys accelerometers in this corpus are stored at 2000 Hz and repeat every value, so their real rate is far lower and any spectrum computed at the stored rate is wrong above the true Nyquist.
Source code in src/micromotion/validate.py
388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 | |
frame_count ¶
frame_count(n: int, where: str = '') -> list[Finding]
A sample count sitting exactly on a 16-bit boundary, which means a silent truncation.
C3D stores its frame count in sixteen bits, so a conversion through it stops at 65535 and says nothing. Seven sessions in this corpus were nearly deposited that way: 327.7 s of a 360 s recording, complete-looking, with the last thirty-two seconds gone.
Source code in src/micromotion/validate.py
409 410 411 412 413 414 415 416 417 418 419 420 | |
edge_motion ¶
edge_motion(speed, fs: float, where: str = '', edge_s: float = 10.0, baseline_s: tuple[float, float] = (60.0, 300.0), factor: float = 2.0) -> list[Finding]
Whether a recording opens or closes with movement rather than standstill.
A deposited standstill recording should contain standstill and nothing else. In practice exports are trimmed by hand, or not trimmed at all, and what survives at the edges is people walking into position, settling, or being told the recording has ended. That inflates anything computed over a short window and is invisible in a whole-recording median.
speed is a band-limited speed series, whatever sensor it came from. The comparison is
against the recording's own settled interior rather than an absolute threshold, because the
quantity varies by two orders of magnitude across sensors in this corpus.
Returns one finding per affected end, at "warning": settling is a fact about the
recording to be recorded, not necessarily a fault to be fixed.
Source code in src/micromotion/validate.py
423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 | |
settling_time ¶
settling_time(speed, fs: float, baseline_s: tuple[float, float] = (60.0, 300.0), factor: float = 1.5, window_s: float = 5.0, max_s: float = 120.0) -> tuple[float, float]
How long each end of a recording takes to reach its settled level, in seconds.
Returns (head, tail): the time to trim from the start and from the end so that what
remains is within factor of the recording's own settled interior. Zero means that end is
already settled. Use it to choose a trim rather than guessing one: a fixed twelve seconds
was not enough for any recording in the collection it was chosen for.
The search stops at max_s and returns that value, which should be read as "still moving
when the search stopped" rather than as a measurement.
Source code in src/micromotion/validate.py
462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 | |
duplicate_files ¶
duplicate_files(paths, where: str = '') -> list[Finding]
Files in a set that are byte-identical.
A record shipped with the same recording under two names, one of them in no manifest, left over from a rename. It inflated the record and would have inflated any count taken by listing the directory.
Source code in src/micromotion/validate.py
497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 | |
too_still ¶
too_still(x, fs: float, where: str = '', floor_mm_s: float = HUMAN_FLOOR_MM_S, band: tuple[float, float] | None = None) -> list[Finding]
A trace too still to be a body: a tripod, a floor marker, a mount.
Every other check here asks whether a recording moved WRONGLY. This one asks whether it moved at all, and it exists because a name list cannot be trusted to be complete.
THE CASE THAT MOTIVATES IT, twice over. A corpus of optical recordings carried 84
tripod traces as participants for months, because the label pattern matching
equipment missed St01, ST1 and Tripod. Fixing the pattern was not enough: on
2026-08-19 four more entered the same corpus as static1 to static4, because the
pattern allowed a bare static and st1 but not a trailing digit on the word it
already contained. Both times the numbers looked ordinary -- a median is robust to a
tripod, so nothing downstream complained -- and both times what found them was a
check on the PHYSICS rather than on the name.
Use it that way. Run it over everything, and treat what it returns as a list of labels to add to whatever pattern you filter by, not as a filter itself: a check that quietly removes data reads as coverage.
Requires positions in millimetres. Returns nothing for a series too short to filter,
which is deliberate -- a short recording is a different complaint, and
:func:frame_count makes it.
Source code in src/micromotion/validate.py
533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 | |
validate_series ¶
validate_series(x, t=None, documented_hz: float | None = None, where: str = '', min_finite: float = 0.8, expect_positions: bool = True) -> list[Finding]
Run every applicable check on one series and return what is wrong with it.
x is (n_samples, n_axes) or one-dimensional. Pass t to check the clock, and
documented_hz to check it against what the record claims. Set expect_positions
False for a sensor whose output may legitimately be zero in every axis at once.
Source code in src/micromotion/validate.py
582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 | |
errors ¶
errors(findings) -> list[Finding]
Just the findings that should stop a build.
Source code in src/micromotion/validate.py
604 605 606 | |
raise_on_error ¶
raise_on_error(findings) -> None
Raise if anything in findings is an error, listing all of them.
A build that continues past a finding produces a plausible wrong number, which is the exact failure this module exists to prevent.
Source code in src/micromotion/validate.py
609 610 611 612 613 614 615 616 617 | |
Comparing over equal windows¶
micromotion.windows ¶
Comparing conditions over windows of the same length.
A paired contrast between two conditions is only about the conditions if the two were measured over comparable stretches of time. That is easy to assume and easy to get wrong, because segment duration usually follows the STIMULUS rather than the design: a running order plays a track for as long as the track lasts and leaves whatever gap it leaves, so "music" and "silence" end up different lengths without anyone deciding they should be.
It matters because a quantity averaged over a longer window settles further toward its middle. Two conditions measured over unequal windows are estimated with unequal smoothing, so the difference between them is changed by the schedule rather than by the participants. In the corpus this package was written for, the effect was to SUPPRESS: equalising the windows raised a music-versus-silence contrast at 0.5-1 Hz from +1.70 to +2.10 per cent and moved two further frequency bands from null to significant. Elsewhere in the same corpus it worked the other way and produced a false positive below 0.5 Hz. The direction is not predictable, which is the argument for checking rather than reasoning about it.
Two functions, in the order they should be used. :func:balance asks whether the windows differ by
condition at all, and is the one-line check that would have caught this years earlier than it was
caught. :func:equalise truncates every segment to a common length so that they do not.
>>> import numpy as np
>>> onset = np.array([0., 60., 120., 180.])
>>> offset = np.array([45., 105., 165., 240.]) # silences 45 s, music 60 s
>>> cond = np.array(["s", "s", "s", "m"])
>>> b = balance(onset, offset, cond)
>>> b.balanced
False
>>> new_offset = equalise(onset, offset)
>>> (new_offset - onset).tolist()
[45.0, 45.0, 45.0, 45.0]
CAP PER GROUP, NOT ACROSS THE STUDY. Where segments come from several recordings, editions or
sessions, pass by= so each group is equalised against its own shortest segment. A single cap
across the whole study equalises the conditions AND shortens every window, and those pull in
opposite directions: applied to six recording sessions whose segments ran from 20 to 180 s, a flat
cap made a real effect look like it had collapsed, purely by discarding the signal of the sessions
that had recorded longest.
Balance
dataclass
¶
What :func:balance found. Truthy when the windows are comparable.
by_condition maps each condition to (n, median, min, max) in seconds. ratio is the
longest condition median over the shortest, so 1.0 is perfect balance. balanced applies the
tolerance the caller asked for.
Source code in src/micromotion/windows.py
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 | |
balance ¶
balance(onset_s, offset_s, condition, by=None, tol: float = 0.05) -> Balance
Do the conditions being compared occupy windows of the same length?
Returns a :class:Balance, which is falsy when they do not. tol is how far the ratio of
condition medians may sit from 1.0 before the answer is no; the default of 0.05 allows the few
per cent that rounding a segment table to whole seconds produces.
PASS by= WHENEVER THE SEGMENTS COME FROM MORE THAN ONE RECORDING, SESSION OR EDITION, and
read the per-group result rather than the pooled one. Pooling averages the imbalances and can
hide a severe one almost completely. On the corpus this was written for, the pooled ratio over
six editions is 1.07 — a few per cent, easy to wave away — while one edition inside it sits at
9.00 and three others between 1.33 and 1.60. The pooled figure is the one that reassures, and
it is the wrong one. With by, ratio is the worst group's and groups holds each.
Run this before any paired contrast between conditions. It costs one line and answers a question that is otherwise invisible: nothing about a table of onsets and offsets announces that one condition is systematically longer than another, and no amount of checking the analysis code will reveal it, because the fault is in the design of the running order rather than in the arithmetic.
Source code in src/micromotion/windows.py
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 | |
equalise ¶
equalise(onset_s, offset_s, cap_s: float | None = None, by=None)
Truncate every segment to a common length, returning new offsets.
Each segment keeps its own onset and is cut from the start, so the equalised window is the beginning of the segment rather than a slice from its middle. That is deliberate: the start is the part every segment has.
cap_s sets the length; by default it is the shortest segment present, which keeps as much
of every segment as can be kept while making them equal. by groups the segments, so each
group is capped against its own shortest segment rather than against the study's — see the
module docstring for why a single cap across groups is the wrong tool.
Only offsets are returned. The onsets are unchanged, so the caller's other columns stay aligned and nothing needs reordering.
Source code in src/micromotion/windows.py
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 | |
Reading files¶
micromotion.io ¶
Readers for the layouts this corpus actually uses.
Each returns a :class:~micromotion.record.MotionRecord. Gap sentinels become NaN, units
are recorded rather than assumed, and the sampling rate is measured wherever the file
carries a timebase.
:func:read dispatches on content, not on the extension, because the extension lies: the
balance-board dumps are headerless and space-delimited whatever they are called, and the
Qualisys family puts three different header shapes behind one name.
Y_UP_COLLECTIONS
module-attribute
¶
Y_UP_COLLECTIONS: tuple[str, ...] = ()
Datasets whose vertical axis is Y rather than Z. Empty by default.
A system's frame is a property of how it was configured for a session, not of the system: the same OptiTrack rig can produce Y-up files for one study and Z-up for another. Where a dataset is known to be Y-up, name it here rather than compensating in downstream analysis.
Prefer rotating the data at source (X -> X, Y -> -Z_old, Z -> Y_old — a rotation, not an
axis swap). If you do, empty this in the same change: a reader that still claims Y for
rotated files hands every caller a horizontal axis as the vertical one, silently.
BOARD_MM
module-attribute
¶
BOARD_MM = (433.0, 238.0)
Wii Balance Board sensing area, width by depth, in millimetres.
Multiplying the normalised centre of pressure by these gives millimetres, which is what makes balance data comparable with the optical collections.
read_qualisys ¶
read_qualisys(path: str, drop_gaps: bool = True) -> MotionRecord
Qualisys or Qualisys-style TSV, in all three header shapes found in the corpus.
The shapes differ in what follows the ten KEY<TAB>value metadata lines:
- a
<marker> Xcolumn-name row, then data (2012, 2015, 2017, 2018, 2019, HpSp); - nothing, data begins immediately (2022, Bishop 2020);
- a
Frame/Timecolumn-name row, then data (MocapNoiseFloor, Solberg 2016).
The shape is detected by trying to parse the eleventh line as numbers, so a file that was exported differently still reads correctly.
Source code in src/micromotion/io.py
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 | |
read_sverm ¶
read_sverm(path: str) -> MotionRecord
Sverm plain-header export: Time then s<n>_head_{x,y,z}.
A separate reader because it shares nothing with the Qualisys export but the units. It silently returned zero series to a corpus analysis until one was written for it.
Source code in src/micromotion/io.py
151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 | |
read_ax3 ¶
read_ax3(path: str) -> MotionRecord
Axivity AX3 export with a ts,x,y,z header.
The rate is measured from the timestamps and differs meaningfully between files: it is a property of the physical logger, spanning 191.3-207.7 Hz across the 2024 units, and it can differ between recording sites in the same direction as the effect a study record's headline claim reports.
Source code in src/micromotion/io.py
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 | |
read_cwa ¶
read_cwa(path: str) -> MotionRecord
Axivity AX3 raw .cwa, the file the logger itself writes.
:func:read_ax3 reads the ts,x,y,z export made from one of these. Studies deposit the raw
file instead often enough that a reader is worth having: the Solberg 2015 dance record ships
twenty .cwa and no export, so without this every reuser writes their own decoder or installs
a second toolbox to get at deposited data.
THE TIMEBASE IS THE POINT, and it is why this returns wall-clock seconds rather than seconds
from the start. Each 512-byte block carries its own clock reading and a timestampOffset
saying which sample inside the block that reading applies to, so the time axis is rebuilt block
by block: a dropped block then leaves a gap in t instead of silently shifting every sample
after it forward. t is seconds from midnight on the recording's own day, which is the form a
session log tends to be written in; t0 carries that day.
TWO PACKINGS EXIST AND BOTH ARE IN THE WILD. The unpacked one stores three int16 per sample;
the packed one stores three 10-bit values and a shared 2-bit exponent in four bytes, and its
values must be scaled by 2**exponent to land on the same scale as the unpacked form. Getting
that exponent wrong is not obvious from the numbers --- it scales the whole recording by a
constant, which leaves every correlation and every rank statistic untouched. Check a decode with
:func:identify_acceleration_unit, which reads gravity off a still stretch and is the test this
reader was verified with: the twenty deposited Solberg files read 0.96 to 1.06 g.
Returns:
| Name | Type | Description |
|---|---|---|
MotionRecord |
MotionRecord
|
|
MotionRecord
|
measured rate in |
Source code in src/micromotion/io.py
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 | |
channel_rate ¶
channel_rate(t: ndarray, x: ndarray) -> float
How often a channel actually advances, not how often the file has a row for it.
A multi-sensor log is usually an interleaved union of streams that update at different rates, so the spacing of its rows belongs to no device. On one Physics Toolbox file the rows arrive at about 426 Hz while the accelerometer updates at 51 Hz and the fused channel at 15 Hz. Taking the row rate as the sensor rate is how a 100 Hz resampling grid and a 12.5 Hz decimation both came to be quoted as sampling rates in this project.
Counted as value changes over the elapsed span, which is robust to a channel repeating a
value and to occasional dropouts. The tempting alternative, 1 / median(diff(t)) at the
change points, is not: where updates arrive in bursts it returns the within-burst spacing,
which on that same file gives 680 Hz.
x may be one channel or an (n, k) block; the first column is used.
Source code in src/micromotion/io.py
372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 | |
channel_resolution ¶
channel_resolution(x, need: float | None = None) -> dict
The quantisation step of a channel, against the amplitude you mean to measure.
The neighbour of :func:channel_rate, and the same class of mistake. That one asks whether a
channel updates fast enough to carry your band; this one asks whether it resolves finely enough
to carry your amplitude. A channel can satisfy either and fail the other.
The case it was written for. Delsys EMG sensors carry a three-axis accelerometer alongside the muscle channel, and in a 2017 standstill recording those accelerometers step by 0.0395 m/s², identically on all twelve axes. Each axis therefore holds between eight and forty-nine distinct values across 720000 samples, where the muscle channel beside it holds about fifty thousand. The band-limited head acceleration being measured has a median of 0.033 m/s², so one quantisation step was larger than the entire signal. Correlating those accelerometers against anything returned about 0.03, which reads exactly like a real null and is not one.
An accelerometer bundled with another sensor is specified for that sensor's purpose. These are there to tell a standing muscle from a walking one and they do that well. Check the step against the amplitude you need before planning an analysis on a secondary channel.
Returns step, the modal spacing between adjacent distinct values, levels, how many
distinct values the channel holds, span, and ratio where need is given: the
amplitude you asked for divided by the step. A ratio below about 10 means the quantisation is a
visible part of your measurement; below 1 the signal is inside one step and cannot be recovered.
x may be one channel or an (n, k) block, in which case each column is reported.
Source code in src/micromotion/io.py
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 | |
read_phone ¶
read_phone(path: str, trim_clap_s: float = 0.0, *, channel: str = 'accel', trim_start_s: float | None = None, trim_end_s: float | None = None) -> MotionRecord
Physics Toolbox phone log, raw app export or cleaned tab-separated form.
The variant is detected from the first two lines, so both the app's own semicolon/
decimal-comma CSV and the cleaned TSV used by the deposited pipeline read correctly.
See _read_physics_toolbox_raw for what the raw format does that plain CSV readers
get wrong.
Which channel. A Physics Toolbox log carries two accelerations and they are not interchangeable.
channel="accel" (the default) reads gFx/gFy/gFz: the accelerometer itself,
total specific force including gravity, written in g and converted to m/s^2 here. Its
magnitude sits at about 1 g on a phone at rest. It updates at the accelerometer's own rate,
commonly 50 Hz but anywhere from 15 to 455 Hz depending on how the app was configured.
channel="fused" reads ax/ay/az: linear acceleration, already in m/s^2, with
gravity removed by fusing the accelerometer with the gyroscope and magnetometer. They are
not in g, whatever an older data dictionary said; reading them as g inflates every quantity
of motion by 9.80665.
Use accel for anything about how much something moved. The fusion cannot output faster
than its slowest input, so it advances at the gyroscope's rate — about 15 Hz — and in a
0.2-5 Hz band most of what it carries is its own noise floor rather than the body. At
standstill amplitudes the body sits below that floor, and the floor differs between
handsets, so two phones recording the same stillness disagree by a factor that looks like a
device calibration difference and is not. That is why accel is the default: on one
session the fused channel put a chest sensor at 1.65 mm/s against the accelerometer's 6.74,
and made a head sensor look 4.6 times more active than the chest when the true ratio is
1.12.
Use fused for the fusion's tilt correction, accepting the ceiling. Tilt is real, since
rotating a sensor swings gravity across its axes and no high-pass can remove that, but
measure it before assuming it dominates: reconstructed from a gyroscope on one chest
recording it accounted for 6 per cent of the accelerometer's band-limited content, and
removing it did not move that channel toward the fused one.
Whichever you choose, the other is in meta["extra"] and both are in the same unit, m/s^2.
meta["channel_rates"] gives each sensor's own rate, which is not the rate of the file's
rows: an interleaved log's row spacing belongs to no device. See channel_rate.
The rate is neither constant nor the nominal one. Physics Toolbox delivers whatever the
Android sensor stack gives it, so a log requested at 100 Hz arrives between roughly 100 and
170 Hz with millisecond-scale jitter, and differs between handsets recording the same event.
fs here is the measured mean rate over the span. Long dropouts are common: logging
stops when the app is backgrounded or the screen sleeps, and resumes silently, leaving gaps
of tens of seconds inside a file that otherwise looks continuous. Check meta["gaps"]
before treating a file as one recording; resample onto a uniform grid before filtering.
trim_clap_s drops that many seconds from each end. trim_start_s and
trim_end_s override it per end, so an opening clap can be removed without discarding
good data at the close: pass trim_start_s=35, trim_end_s=0.
Trim before plotting or transforming. A recording that opens with a synchronisation clap carries a transient that is a timing marker, not movement, and it can be two orders of magnitude above the standstill it precedes, reaching 10.29 m/s² in one recording against a body maximum of 0.155. Left in, it sets the y-axis of any plot, dominates any spectrum, and gives a peak detector a transient that is not a breath.
The settling that follows is slower and also worth removing. Measure it rather than guessing: both posture and heart rate can still be moving well beyond ten seconds.
Source code in src/micromotion/io.py
445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 | |
read_equivital ¶
read_equivital(path: str) -> MotionRecord
Equivital physiology CSV: accelerometer, ECG, respiration or RR.
Four of the five files per participant are delimited by comma-and-space and the fifth by a bare comma, so the separator is handled rather than assumed. The accelerometer is in raw counts, calibrated to g by the median vector magnitude, since a person standing still averages one g.
The rate is measured from the timestamps by span over count. Taking the median interval instead returns exactly 250 Hz for a recording that runs at 256, because the timestamps are rounded to whole milliseconds; that error was live in the deposited quantity of motion.
Source code in src/micromotion/io.py
603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 | |
read_balance_board ¶
read_balance_board(path: str) -> MotionRecord
Wii balance board dump: headerless, space-delimited, irregularly sampled.
Eight columns: timestamp in milliseconds, four corner load cells, total load, and centre of pressure in x and y, all normalised to 0-1.
Samples with no load carry a centre of pressure of exactly (0.5, 0.5), the midpoint of the board, which is a plausible-looking value for nobody standing on it. Those are returned as NaN. They are 13.8 per cent of the HpSp balance data.
Source code in src/micromotion/io.py
656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 | |
read_fnirs ¶
read_fnirs(path: str) -> MotionRecord
Artinis Brite wide export: 96 haemoglobin channels plus two IMUs.
Only the accelerometer and gyroscope are returned as data; the haemoglobin channels
are kept in meta because they are not motion. The accelerometer is in milli-g.
Column 0 is a sample index, not time, so the rate comes from the header. The absolute
start can be decoded from the TimeStampHi/Lo pair if it is needed.
Source code in src/micromotion/io.py
693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 | |
sniff ¶
sniff(path: str) -> str
Identify a file's layout from its first lines.
Source code in src/micromotion/io.py
744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 | |
read ¶
read(path: str, **kw) -> MotionRecord
Read any corpus motion file, dispatching on content.
Source code in src/micromotion/io.py
773 774 775 | |
micromotion.record ¶
The common structure every reader returns.
The corpus holds sixteen distinct file layouts across fourteen collections, and each one has already cost an analysis at least once: an axis convention that differs in a single edition, a rate that three documents disagree about, a gap code that looks like a valid measurement. A reader's job is to absorb all of that and hand back the same object regardless.
MotionRecord
dataclass
¶
One recording, in a form the rest of the package can use without special cases.
Attributes:
| Name | Type | Description |
|---|---|---|
data |
ndarray
|
(n_samples, n_channels) float array. Gaps are NaN, never a sentinel. |
fs |
float
|
Sampling rate in Hz, measured from the file where a timebase exists and taken from the header only where it does not. |
channels |
list[str]
|
Column names, one per column of |
kind |
str
|
|
unit |
str
|
|
vertical |
str
|
Which axis letter is up. Not always Z: the 2019 championship is Y-up, alone in the corpus. |
t |
ndarray | None
|
Timestamps in seconds where the file carries them. |
t0 |
object | None
|
Absolute start time where the file carries one. |
source |
str
|
Path the record was read from. |
meta |
dict
|
Anything else the header held, unmodified. |
Source code in src/micromotion/record.py
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 | |
markers
property
¶
markers: list[str]
Marker names, for records whose channels are <marker> X/Y/Z triplets.
marker ¶
marker(name: str) -> np.ndarray
The (n_samples, 3) block for one marker, by name.
Reading marker columns by name rather than by position is not a nicety. Six HpSp files break the documented 22-marker order, and the README's positional quick-start mis-assigned markers in every one of them until it was rewritten.
Source code in src/micromotion/record.py
77 78 79 80 81 82 83 84 85 86 87 | |
missing_fraction ¶
missing_fraction() -> float
Proportion of the array that is NaN.
Source code in src/micromotion/record.py
89 90 91 | |
qom ¶
qom(**kw)
Quantity of motion for this record, with kind and unit already filled in.
Passing a whole record computes across every channel at once, which is rarely what you want for a multi-marker file; select a marker first.
Source code in src/micromotion/record.py
93 94 95 96 97 98 99 100 101 102 103 | |
Motion capture¶
micromotion.mocap ¶
Motion-capture I/O and cross-modality utilities.
Pure numpy/scipy helpers ported from the "still standing" and Westney-comparisons studies:
- :func:
read_qtm_tsv-- a single robust reader for Qualisys Track Manager (QTM) TSV exports, consolidating the four/five near-duplicate loaders that were copy-pasted across the study scripts. - :func:
compare_modality_envelopes-- resample two motion envelopes onto a common per-second grid and correlate them (e.g. video-pose vs mocap validation). - :func:
dominant_frequency-- the dominant spectral peak of a signal within a band, via Welch.
.. note::
:func:compare_modality_envelopes deliberately takes precomputed 1-D
motion envelopes rather than computing quantity-of-motion internally, so
this module stays independent of the QoM machinery. The natural producer
of such envelopes is musicalgestures._qom.band_limited_qom followed by
a per-second binning (envelope / bin_series), arriving in a
sibling PR.
Source: still standing study and Westney-comparisons study (Jensenius).
read_qtm_tsv ¶
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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
path
|
str
|
Path to the |
required |
Returns:
| Name | Type | Description |
|---|---|---|
tuple |
|
Source code in src/micromotion/mocap.py
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 | |
dominant_frequency ¶
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.
A bare maximum, and it stays one, because published figures were computed
with it. It now warns when the band holds no peak at all, by
:func:~micromotion.spectral.is_band_floor, and returns the value
unchanged. Take that warning seriously on position and speed signals, whose
spectra fall steeply: on synthetic 1/f series with nothing in a 0.7-2.2 Hz
band this returns a median 1.07 times the lower edge and lands exactly on
that edge on none of forty of them, so an audit looking for values that sit
on a band boundary would not have found one. Use
:func:~micromotion.spectral.spectral_peak where a NaN is preferable to a
number, and :func:~micromotion.spectral.band_edge_sweep to test an
estimate that has already been computed.
Source: Westney-comparisons study (Jensenius), extended motion-feature analysis (motion dominant frequency of vertical trunk position).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
ndarray
|
1-D input signal. |
required |
fs
|
float
|
Sampling rate in Hz. |
required |
band
|
tuple
|
|
(0.3, 4.0)
|
Returns:
| Name | Type | Description |
|---|---|---|
float |
The dominant frequency in Hz, or |
Source code in src/micromotion/mocap.py
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 | |
compare_modality_envelopes ¶
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).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
env_a
|
ndarray
|
First 1-D motion envelope. |
required |
env_b
|
ndarray
|
Second 1-D motion envelope. |
required |
fs_a
|
float
|
Sampling rate of |
required |
fs_b
|
float
|
Sampling rate of |
required |
Returns:
| Name | Type | Description |
|---|---|---|
dict |
|
Source code in src/micromotion/mocap.py
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 | |
Spectral¶
micromotion.spectral ¶
Spectral helpers: finding the physiological peaks in a motion signal.
A body-worn accelerometer on someone standing still picks up two rhythms that are not movement in the intentional sense. Respiration sits at roughly 0.2-0.4 Hz and the ballistocardiac impulse, the recoil of the heart ejecting blood, at roughly 0.8-1.8 Hz. Both are inside or adjacent to the micromotion band, so isolating postural motion means locating them first.
CARDIAC_BAND
module-attribute
¶
CARDIAC_BAND = (0.7, 2.2)
Hz. 42-132 bpm, which covers rest through mild exertion.
DEFAULT_EDGE_FACTORS
module-attribute
¶
DEFAULT_EDGE_FACTORS = (0.7, 0.85, 1.0, 1.15, 1.3)
What :func:band_edge_sweep multiplies a band's lower edge by when given no edges.
SHARE_RULES
module-attribute
¶
SHARE_RULES = ('trapezoid', 'sum')
The quadrature rules a share can be taken under. See :func:band_share.
SHARE_INTERVALS
module-attribute
¶
SHARE_INTERVALS = ('closed', 'half_open')
The band-edge conventions a share can be taken under. See :func:band_share.
is_band_floor ¶
is_band_floor(f, p, band: tuple[float, float]) -> bool
Whether the largest value in band is where the band starts rather than a rhythm.
The test is :func:peak_from_spectrum's, so the package holds one peak rule and not two: is
there an interior local maximum of the spectrum divided by a log-log straight-line fit across
the band. True means there is not, and that whatever a bare argmax returned is a property
of the search band rather than of the body.
GIVE THIS AN UNFILTERED SPECTRUM. A band-pass applied before the transform builds its own
rising skirt inside the passband, and that skirt survives dividing out the log-log slope. On
twenty synthetic 1/f series with nothing in the band, the rule called the filtered spectrum a
peak 16 times and the raw spectrum once at a single-segment Welch, and 17 against 6 at the
averaging :func:_floor_spectrum uses. Every caller in this package therefore runs it on the
raw segment even where the estimate itself is taken from a filtered one. It is a difference of
sensitivity rather than of verdict: at the call sites here, where a warning follows if ANY
window is flagged, both choices still warn on all twenty series.
GIVE IT AN AVERAGED SPECTRUM TOO, which is the less obvious half. The rule asks whether any bin
stands a factor above a fitted slope, and in a lightly averaged Welch some bin always does, by
chance. On the same twenty 1/f series over a 0.7-2.2 Hz band this correctly returns True on
20 of 20 at about thirty Welch averages, on 14 of 20 at fifteen, and on 0 of 20 at the five
that mocap.dominant_frequency's own nperseg=2048 leaves — which is why every caller here
recomputes its own diagnostic spectrum instead of reusing the one the estimate came from. The
failure is one-sided: too few averages makes this MISS a band floor, never invent one. A short
record, or a band narrow enough that resolving it uses up the record, therefore gets a weaker
test rather than a wrong one.
Source code in src/micromotion/spectral.py
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | |
cardiac_peak ¶
cardiac_peak(x, fs: float, window_s: float = 60.0) -> float
Dominant frequency in the cardiac band, in Hz.
Pass the acceleration magnitude. Multiply by 60 for beats per minute. Returns NaN if the recording is too short for the band to be resolved.
This is a bare maximum inside the band, unlike :func:spectral_peak, and it stays one
because published beats-per-minute figures in this corpus were computed with it. What it
gained instead is a warning: when the band holds no peak it says so and returns the value
anyway. On synthetic 1/f series with nothing in the cardiac band it returns a median 1.05
times the 0.7 Hz lower edge, and lands exactly on that edge on only a quarter of them, so a
check for values sitting on the boundary would have passed most of them. Use
:func:spectral_peak where a NaN is preferable to a number, and :func:band_edge_sweep to
settle whether an estimate already computed is following its own band edge.
Source code in src/micromotion/spectral.py
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 | |
respiratory_peak ¶
respiratory_peak(x, fs: float, window_s: float = 120.0) -> float
Dominant respiratory frequency, in Hz. Multiply by 60 for breaths per minute.
This is measured in the time domain, from :func:detect_breaths, and NOT as a peak in a
periodogram. It used to be the latter and the change corrects a wrong number rather than
trading one convention for another.
Why, in short: a periodogram of belt or body motion is red, so the breathing bump sits on a much larger downward slope and never becomes the global maximum inside the band. Measured on the one dataset in this corpus with a ground truth -- Stillness2025's sixteen thoracic belts at 25.6 Hz -- the old version returned a median 7.5 breaths per minute where the belts' own breath timing gives 16.8 and where a resting adult breathes 12-20. It was not a calibration offset: it ranked those sixteen participants at Spearman -0.32 against their own breath timing, so it carried no usable information about who was breathing faster.
Four repairs were measured and all four rejected, which is why the spectral approach was abandoned rather than patched. Raising the band floor to 0.20 Hz leaves a 3.4 breaths-per-minute gap. Band-passing before the periodogram does not rescue it. The most prominent local maximum instead of the global one reaches Spearman +0.22. Dividing out a fitted power law before taking the maximum reaches +0.26 and biases the median high, to 21.5.
ONE REASON GIVEN FOR THE SECOND OF THOSE WAS WRONG, and it matters because it is the reason
people go on proposing it. This docstring used to say that band-passing "changes nothing
whatever, and cannot, because the maximum inside a band is unaffected by filtering inside that
same band". The maximum inside a band is NOT unaffected: a Butterworth is not flat inside its
own passband, its rising skirt reaches in, and multiplying a falling spectrum by that skirt
moves the maximum up. Measured on twenty synthetic 1/f series over a 0.12-0.40 Hz band, the
unfiltered maximum sits at 1.11 times the lower edge, a fourth-order zero-phase band-pass over
the same band moves it to 1.25 and a second-order one to 1.39, and the filtered and unfiltered
answers agree on 6 and 4 of the twenty. So band-passing moves the number by about a fifth and
still returns the edge; the repair fails because the answer is the edge either way, not
because filtering is a no-op. See :func:band_edge_sweep.
:func:cardiac_peak still uses the periodogram and is right to: its band sits above the slope
and the ballistocardiac impulse is a genuinely prominent peak, giving a median 75 bpm with an
interquartile range of 70-82 on a year of chest-phone data.
window_s is accepted for backward compatibility and is unused; breath detection does not
need a spectral window.
Source code in src/micromotion/spectral.py
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 | |
band_power ¶
band_power(x, fs: float, band: tuple[float, float], window_s: float = 60.0) -> float
Integrated power between two frequencies.
Trapezoid quadrature over a closed interval, [lo, hi]. A ratio of two of these is not
the same number as a ratio computed by summing bins over [lo, hi), which is what
:func:~micromotion.physio.spectral_band_fractions does; :func:band_share names the
convention it uses and can express either.
Source code in src/micromotion/spectral.py
162 163 164 165 166 167 168 169 170 171 172 173 174 | |
peak_from_spectrum ¶
peak_from_spectrum(f, p, band: tuple[float, float], require_peak: bool = True, min_excess: float = 2.0) -> dict
The peak rule, for a caller that already has a spectrum.
spectral_peak is this with a Welch in front of it. It exists separately because several
analyses compute one spectrum and read three bands off it, and recomputing the transform per
band to get at the rule would be both wasteful and an invitation to reimplement it locally.
A rule that is easier to copy than to import gets copied; that is how the reference markers
got into one of these analyses twice.
See spectral_peak for what the rule is and why.
Source code in src/micromotion/spectral.py
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 | |
spectral_peak ¶
spectral_peak(x, fs: float, band: tuple[float, float], window_s: float = 60.0, require_peak: bool = True, min_excess: float = 2.0) -> dict
Peak frequency in a band, or NaN when the band contains no peak.
Every spectrum has a maximum somewhere inside any band you choose, and on a falling
spectrum that maximum is the lowest bin. It is not a rhythm, it is the slope. With
require_peak the result is NaN unless the band holds something that stands above its own
background: an interior local maximum of the spectrum divided by a log-log straight-line fit
across the band, exceeding min_excess. is_peak and excess in the returned dict say
which happened and by how much.
WHY THE BASELINE IS A FITTED SLOPE AND NOT THE BAND MEDIAN. Measuring a peak against the median of its own band assumes the band is flat. Over a 1/f spectrum it is not: the median is dragged down by the high-frequency end, so the low bins clear any threshold without being peaks. On plain 1/f noise the band's largest value scores 3.7 against the median and 1.3 to 1.7 against the fitted slope, where a real rhythm scores 5 and up against either.
THE SIGNAL-TO-NOISE RATIO DOES NOT CATCH THIS, and an earlier version of this docstring advised using it that way. It is wrong in the one case that matters. The lowest bin of a 1/f spectrum has both the most power and the highest power-over-band-median of anything inside a band drawn above the knee, so raising an SNR threshold SELECTS the artefact rather than excluding it. In a year of daily standstill recordings, tightening the threshold from nothing to 5 took the share of days whose "respiration rate" sat exactly on the band floor from 21 per cent to 32, and moved the median from 10.5 breaths a minute to 9.0. Four analyses in the Oslo Standstill corpus reported a band edge as a measurement before this was found, one of them on 662 of 930 values, and one of those numbers had reached a book.
REJECTING THE EDGE BIN IS NOT ENOUGH EITHER. On a monotone slope, refusing the first bin moves the maximum to the second: the same 662 of 930 became 198 of 268, one bin along. Only asking whether the thing is a peak at all separates the two cases.
IT IS MORE ACCURATE AND NOT ONLY MORE CAUTIOUS. Dividing out the slope finds peaks the raw maximum misses: on 1/f noise plus a modest 0.25 Hz tone, the old answer is the band floor and this returns 0.25.
WHAT IT COSTS. It is conservative, and a rhythm weaker than about half that returns NaN even
though it is really there. That is the right direction to fail in -- a missing value rather
than a wrong one -- but it is a false negative, so a rate aggregated over many recordings will
be missing its weakest cases and the count of NaNs is part of the result rather than a
nuisance. Lower min_excess to trade the other way, knowingly.
snr is still the peak over the band median, unchanged, because callers report it.
require_peak=False restores the old behaviour exactly, for a caller who wants the largest
value in a band and knows that is what they are asking for.
Source code in src/micromotion/spectral.py
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 | |
band_edge_sweep ¶
band_edge_sweep(signals, fs: float, band: tuple[float, float], *, estimator=None, edges=None, reference=None) -> dict
Move the lower edge of a search band and see whether the answer follows it.
The question this answers is not "did the estimate land on the boundary" but "is the estimate the boundary". Those are different, and the second is the one that catches things. An estimator that reports the largest peak inside a band, run on a spectrum that falls steeply, returns a number that MOVES WITH the band edge at a near-constant multiple of it -- not the edge itself, a plausible-looking interior value a fifth or a third above it. Checking whether values sit ON a boundary passes that clean. Moving the boundary does not.
Pass one signal, or a sequence of signals from a collection. estimator is called as
estimator(item, fs, (lo, hi)) and must return a frequency in Hz; it defaults to the bare
in-band maximum, which is the thing usually under suspicion, and it can be any callable, so an
estimator from outside this package -- a whole pipeline, reading video -- can be audited by the
same rule. Items are passed through untouched, so they need only be whatever that callable
accepts. Warnings raised inside it are suppressed, since the sweep deliberately calls it in the
regime where it complains.
edges are the lower edges to try, defaulting to :data:DEFAULT_EDGE_FACTORS times
band[0]. The upper edge is held fixed throughout: this tests one boundary at a time.
Returns a dict:
edges, andfreq, the answer at each edge (the median across signals if there are several), withfreq_by_signalshaped(n_edges, n_signals)ratio,freq / edges, which is flat for an estimate that is the edgefactor, the least-squares multiple of the edge through the originrss_edgeandrss_constant, how well "the answer is c times the edge" and "the answer is a fixed frequency" each fit the sweepfollows, which isrss_edge < rss_constant: the edge explains the answers better than a rhythm doesragainstreferenceat each edge, andr_max, when a reference frequency per signal is given
WHY THE VERDICT COMPARES TWO FITS rather than thresholding a slope. A genuine rhythm returns
the same frequency at every edge below it, so a constant fits perfectly and the edge fits
badly. An estimate that is the edge fits c * edge and not a constant. The comparison needs
no threshold, and it covers both shapes of the failure at once: an answer sitting exactly on
the boundary is this with factor at 1.0.
WHAT IT CANNOT DO. Push the edge above the rhythm and every estimator follows it, correctly,
because the rhythm is no longer in the band. So keep the swept edges below the frequency you
expect; the default range spans 0.7 to 1.3 times a band edge that was presumably chosen to sit
below it. A follows verdict from edges that straddle the answer says nothing.
A single sweep is a strong test of one estimator on one collection. reference makes it
conclusive: if the estimate carries information about the body, it correlates with an
independent measurement of the same quantity at SOME edge. The case this was written from --
a heart rate read from a year of video -- reported 1.24 to 1.33 times its own 0.7 Hz lower
edge, moved from 40 to 116 beats a minute as the edge moved from 0.5 to 1.5 Hz, and never
correlated with a worn reference above 0.21 at any setting.
import numpy as np t = np.arange(0, 300, 0.1) rng = np.random.default_rng(0) tone = np.sin(2 * np.pi * 0.9 * t) + 0.5 * rng.normal(size=len(t)) band_edge_sweep(tone, 10.0, (0.5, 2.0))["follows"] False
Source code in src/micromotion/spectral.py
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 | |
band_rms ¶
band_rms(x, fs: float, band: tuple[float, float]) -> float
Root-mean-square amplitude within a band, in the input units.
Rate-independent, unlike an integrated measure, which makes it the right thing to quote when comparing instruments that sample at different rates.
Source code in src/micromotion/spectral.py
401 402 403 404 405 406 407 408 409 | |
band_power_fraction ¶
band_power_fraction(x, fs: float, bands: dict, window_s: float = 60.0) -> dict
Proportion of total power falling in each named band.
Pass something like {"respiratory": (0.1, 0.5), "cardiac": (0.7, 2.2)}. This is how
the chest-phone cardiac share that motivated the compensated quantity-of-motion variant
was established. The number that finding long circulated as -- 38 per cent -- did not
survive re-measurement under stated bands; the figure that reproduces is 25 per cent of
0.2-5 Hz acceleration power at 0.8-2.5 Hz, on the raw accelerometer channel. That history
is why :func:band_share exists and makes both bands mandatory; prefer it whenever the
result is going to be quoted.
WHICH CONVENTION THIS IS. Trapezoid quadrature over a closed interval, [lo, hi], with the
whole spectrum as the denominator rather than a named band -- the same arithmetic as
:func:band_power and as :func:band_share at its defaults, and NOT the same arithmetic as
:func:~micromotion.physio.spectral_band_fractions, which sums bins over [lo, hi) and
divides by a named band. The two are not interchangeable: on chest-accelerometer standstill
recordings the quadrature rule alone moves a respiratory share by up to 0.034 absolute on a
share of about 0.16, and the closure by up to 0.025. Do not compare a fraction from here with
one from there; :func:band_share can express either, and says which it used.
Source code in src/micromotion/spectral.py
412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 | |
band_share_from_spectrum ¶
band_share_from_spectrum(f, p, *, num_band: tuple[float, float], den_band: tuple[float, float], integrate: str = 'trapezoid', interval: str = 'closed') -> float
The share rule, for a caller that already has a spectrum.
:func:band_share is this with a Welch in front of it, and its docstring states the rule and
the incident behind it. The split mirrors :func:peak_from_spectrum: several analyses compute
one spectrum and read more than one share off it, and a rule that is easier to copy than to
import gets copied.
The deliverability check here is against the spectrum itself: a denominator edge above the highest frequency the spectrum reaches means the share is taken over a truncated denominator, and the warning says so. Everything else -- the mandatory bands, the containment rule, the non-finite warning, and both convention parameters -- is the same.
integrate is "trapezoid" or "sum" and interval is "closed" or
"half_open", exactly as in :func:band_share, which states what the four combinations
mean and what they cost. integrate="sum", interval="half_open" is the convention
:func:~micromotion.physio.spectral_band_fractions implements, and on the same spectrum the
two then agree to floating point; the defaults here are the convention
:func:band_power, :func:band_power_fraction and :func:band_share use.
Source code in src/micromotion/spectral.py
496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 | |
band_share ¶
band_share(x, fs: float, *, num_band: tuple[float, float], den_band: tuple[float, float], window_s: float = 60.0, integrate: str = 'trapezoid', interval: str = 'closed') -> float
Fraction of spectral power in one band over another. Both bands are mandatory.
Returns numerator-band power over denominator-band power, integrated on a Welch spectrum, as a number in 0 to 1. The numerator band must lie inside the denominator band, and neither has a default.
THE TWO BANDS ARE NOT THE WHOLE LABEL: THE ARITHMETIC IS PART OF IT. Two conventions for
turning a Welch spectrum into a band power are in use in this field, and this package used
both before it said so. integrate selects the quadrature rule -- "trapezoid", which
weights the two edge bins by a half, or "sum", which adds the bins in the mask (the
rectangle rule; the bin width cancels in a ratio, so a bin-summed fraction is exactly this).
interval selects what the band edges mean -- "closed", [lo, hi], or
"half_open", [lo, hi). The defaults, "trapezoid" and "closed", are what
:func:band_power, :func:band_power_fraction and this function have always computed, and
what the 25 per cent chest-phone cardiac share was measured under.
integrate="sum", interval="half_open" is the other convention in this package,
implemented by :func:~micromotion.physio.spectral_band_fractions, and it is what the
older cardiac and respiratory composition figures for chest-accelerometer standstill were
taken under; pass it to reproduce those rather than silently restating them.
BOTH CHOICES MOVE REAL SHARES BY MORE THAN THEY LOOK LIKE THEY SHOULD. Measured on chest-accelerometer standstill recordings with the mask held fixed so that only the quadrature rule changed, the respiratory share moved by up to 0.034 absolute on a share of about 0.16 -- over a fifth of the value. Closure costs up to 0.025 absolute on the same recordings, and it costs that much because analysts choose round band edges: at the conventional 60 s window the bin spacing is 1/60 Hz, so 0.40, 0.70, 2.20, 3.0, 5.0 and 8.0 Hz all land exactly on a bin, and closing the interval adds a whole bin at the numerator's upper edge and at the denominator's, which do not cancel. Four analysis scripts in one corpus were carried from the bin-sum convention onto the defaults here, so that every share in that corpus is taken the same way and each is comparable with the rest. The move republished one share from 58 to 59 per cent, another from 16.6 to 15.1, a fold from 3.1 to 3.2, and 18 of the 24 numbers in one table. That is a re-measurement rather than a refactor, and it is why these parameters exist: a convention has to be nameable before a corpus can decide to hold one.
THE TWO AGREE EXACTLY ON A FLAT BAND, which is why they look interchangeable. Where the
spectrum is flat, the trapezoid's half-weighted end bins remove precisely one bin's worth,
so "trapezoid", "closed" equals "sum", "half_open" to floating point. The
disagreement is driven by the slope inside the band, so it is largest at the low end of a
red spectrum -- the respiratory band, on every body-worn sensor here.
A share is comparable only with a share taken the same way. Record the rule and the closure beside the two bands whenever the number is going to be quoted.
WHY THERE ARE NO DEFAULTS. A share is a ratio of two integrals and it moves when either band moves. Four published-looking figures for the share of standstill motion -- 38, 43, 45 and 58 per cent -- circulated in one project and were quoted against one another as though they measured the same thing. Traced to their origins, each came from a hand-rolled fraction with a different, sometimes unstated, denominator; the 58 traces only to its own 0.10-3.0 Hz denominator, and the 45 is untraceable to any measurement at all. A share whose two bands are not stated beside it is not a reportable number. Report the domain (power of what quantity), the site (sensor and placement) and both bands, every time: acceleration power and position power weight the spectrum by a factor of frequency to the fourth relative to each other, so a share of one is not even approximately a share of the other from the same sensor on the same body.
WHAT IT REFUSES AND WHAT IT WARNS ABOUT. A numerator band wider than or outside the
denominator raises, because the result would not be a share. A denominator edge above what
fs can deliver warns, the way :func:~micromotion.filters.bandpass warns, because the
share silently becomes one over a narrower band -- and fs must be the channel's own rate,
not the file's row rate or a resampled grid's; measure it with
:func:~micromotion.io.channel_rate on an interleaved log. A non-finite input warns, as the
filters have since 1.7.0: one NaN makes the whole Welch spectrum NaN, so the share is NaN
rather than mostly right.
IF TWO CHANNELS MUST MEET AT ONE RATE FIRST, resample with
:func:~micromotion.resample.to_rate, which is an anti-aliased polyphase FIR resampler and
refuses to upsample. Plain interpolation onto a slower clock is not a resampler: measured in
the vest decomposition work, interpolating 256 Hz accelerometer axes to a belt's 25.6 Hz
folded high-frequency sensor noise into the cardiac band, which inflates exactly the kind of
share this function computes.
For a spectrum already in hand, :func:band_share_from_spectrum applies the same rule
without recomputing the transform. For several named bands over one common total,
:func:band_power_fraction remains the convenience -- it computes the default convention
here, over the whole spectrum rather than a named denominator. For the bin-sum convention
with a named denominator band, :func:~micromotion.physio.spectral_band_fractions is where
it already lives. This function is the one whose result is meant to be quoted, which is why
it makes the denominator explicit and lets the convention be named rather than assumed.
Source code in src/micromotion/spectral.py
541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 | |
mean_frequency ¶
mean_frequency(x, fs: float, band: tuple[float, float] = (0.1, 5.0), window_s: float = 60.0) -> float
Power-weighted mean frequency within a band, in Hz.
A different statistic from the peak: it moves when energy shifts between frequencies even if the dominant one does not change, so it is the more sensitive of the two to a gradual change in how a person is standing.
Source code in src/micromotion/spectral.py
651 652 653 654 655 656 657 658 659 660 661 662 663 | |
detect_breaths ¶
detect_breaths(x, fs: float, band: tuple[float, float] = RESPIRATORY_BAND, prominence_sd: float = 0.4, min_period_s: float = 2.0) -> dict
Breath timing from a respiration belt or a band-limited motion signal.
Returns peak and trough times and the cycle durations between them. Peaks and troughs are both returned because inhale and exhale are not symmetric, and a rate alone hides that.
Source code in src/micromotion/spectral.py
666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 | |
detect_breaths_adaptive ¶
detect_breaths_adaptive(x, fs: float, band: tuple[float, float] = RESPIRATORY_BAND, vel_frac: float = 0.55, baseline_hz: float = 0.2) -> dict
Breath detection that rejects chest movement which is not breathing.
Peak-and-trough detection, which :func:detect_breaths does, counts any sufficiently
prominent bump. On a belt worn by someone standing that includes postural sway, weight
shifts and swallows, all of which look like small breaths.
This asks a different question. A breath is a sustained rise in chest expansion whose velocity exceeds a threshold set from the signal's own positive-derivative mean, and which crosses an adaptive baseline -- a heavily low-passed copy of the signal rather than zero. A rise that never crosses that baseline did not start from an exhaled state and is discarded. That rejection step is what removes the sway.
Returns inspiration and expiration onsets, cycle durations and the rate. Expiration onset is taken as the end of the rise, which assumes passive expiration.
After Finn Upham's respiration work, reimplemented and used with permission.
Do not reach for this by default, on the evidence available here. It was added on the
expectation that rejecting non-breath rises would beat plain peak detection, and measured
against it that expectation did not hold. On twelve HpSp respiration-belt recordings the
two agree: median error against the spectral estimate 1.82 breaths per minute for both.
On eight chest accelerometers, which is the case the rejection step was
supposed to help, it is markedly worse -- median error against the same participant's belt
10.3 breaths per minute against 3.6 for :func:detect_breaths, over-counting throughout.
That may be the parameters rather than the idea; the velocity threshold is derived from a belt's amplitude distribution and an accelerometer's is not the same shape. It is kept because the approach is sound in its original setting and because someone should be able to tune it, not because it is currently the better detector.
Source code in src/micromotion/spectral.py
686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 | |
Physiology¶
micromotion.physio ¶
Physiology signal features for standstill / micromotion studies.
Two pure numpy/scipy surfaces ported from the "still standing" study:
- :func:
respiration_rate-- windowed breathing rate (breaths per minute) from a respiration waveform, via band-pass filtering and a Welch spectral peak per window. - :func:
spectral_band_fractions-- the fraction of a signal's Welch power falling in each of a set of caller-supplied named frequency bands. This is the generic "cardiorespiratory QoM" spectral-composition diagnostic with the heart-rate/respiration bands supplied by the caller, so the function carries no dependency on any particular physiological sensor. It sums bins over half-open bands, which is not the convention :mod:~micromotion.spectralcomputes shares under; its docstring says what the difference costs and how to ask for either.
Source: still standing study (Jensenius) -- Deichman / Equivital physiology analyses.
DEFAULT_TOTAL_BAND
module-attribute
¶
DEFAULT_TOTAL_BAND = (0.1, 8.0)
Hz. The denominator :func:spectral_band_fractions falls back on when none is given.
It is a fallback and not a standard. Published shares in this corpus rest on it, so it cannot move; a caller who does not name a denominator is warned rather than quietly given this one.
respiration_rate ¶
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.
THIS IS THE SHAPE THAT FAILS WORST ON A FALLING SPECTRUM, and the rate is unchanged but the function now says so. Band-passing to the same band the maximum is then searched in does not leave the maximum where it was: the filter's own rising skirt reaches into the passband and multiplies the falling spectrum by it, so the largest surviving value sits a fixed fraction above the lower edge. Measured on synthetic 1/f series with nothing in a 0.7-2.2 Hz band, the bare maximum returns a median 1.05 times the lower edge and this returns 1.41 times it, and neither lands on the edge itself. Sweeping the lower edge from 0.3 to 1.5 Hz drags this function's answer along at 1.17 to 1.61 times whatever it is set to.
So a window whose band holds no peak returns a plausible interior
frequency, not a boundary value, and an audit for values sitting on a
boundary passes it. Every window is therefore tested with
:func:~micromotion.spectral.is_band_floor on its own UNFILTERED
spectrum, which is the more sensitive of the two: on the filtered one the
skirt reads as a genuine peak on 17 of 20 such windows against 6 of 20 raw.
The function then warns once, naming how many windows failed. At this
function's own defaults, a 1/f series with no breathing in it returns 8 to
11 breaths a minute -- which is where this corpus was reading respiration
rates before any of this was found -- and warns on every one. Use
:func:~micromotion.spectral.spectral_peak where NaN is preferable to a
number, and :func:~micromotion.spectral.band_edge_sweep on a rate
already computed.
Source: still standing study (Jensenius), Deichman respiration analysis
(compute_qom_resp).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
waveform
|
ndarray
|
1-D respiration/breathing waveform. |
required |
fs
|
float
|
Sampling rate in Hz. |
required |
band
|
tuple
|
|
(0.1, 0.6)
|
window_s
|
float
|
Window length in seconds. Defaults to 30. |
30
|
step_s
|
float
|
Hop between windows in seconds. Defaults to 30. |
30
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
|
Source code in src/micromotion/physio.py
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 | |
spectral_band_fractions ¶
spectral_band_fractions(signal, fs, bands, *, total_band=None, 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.
WHICH CONVENTION THIS IS, AND WHAT ELSE IS IN THIS PACKAGE. Power is
bin-summed on the Welch grid over a half-open interval, [lo, hi), at
both the numerator and the denominator. That is one of the two spectral-share
conventions micromotion contains, and the other one is the default
everywhere else: :func:~micromotion.spectral.band_power,
:func:~micromotion.spectral.band_power_fraction and
:func:~micromotion.spectral.band_share integrate with the trapezoid rule
over a closed [lo, hi]. A share from here and a share from there are
not comparable, and until 1.11.0 neither docstring said so.
An earlier version of this docstring claimed the two "yield nearly identical results on the uniform frequency spacing of Welch". They do not. Measured on chest-accelerometer standstill recordings, holding the mask fixed so that only the quadrature rule changed, the respiratory share moved by up to 0.034 absolute on a share of about 0.16 -- over a fifth of the value. Interval closure costs a further 0.025 absolute on the same recordings, because analysts pick round band edges and at the conventional 60 s window (1/60 Hz bins) 0.40, 0.70, 2.20, 3.0, 5.0 and 8.0 Hz all land exactly on a bin, so closing the interval adds a whole bin at both the numerator's and the denominator's upper edge and the two do not cancel. The two conventions agree only where the band is flat, which the low end of a body-worn sensor's spectrum never is.
Nothing here is deprecated and no default has moved. This is a named
convention with a stated arithmetic, and the standstill composition figures
the corpus carried before it settled on the trapezoid-and-closed rule were
computed by exactly it, so this is where they reproduce. To compute a share
under the other convention, or to state
which convention a number was taken under, call
band_share(..., integrate="sum", interval="half_open"), which reproduces
this function on the same spectrum, or leave those parameters at their
defaults to get the trapezoid-and-closed convention. New work that will
quote a number should prefer :func:~micromotion.spectral.band_share: it
makes the denominator mandatory and the convention explicit.
total_band has no silent default. Left unset it falls back to
:data:DEFAULT_TOTAL_BAND, 0.1-8.0 Hz, and warns, because a share is a
ratio of two integrals and a denominator nobody wrote down is the failure
that put four incompatible standstill shares -- 38, 43, 45 and 58 per cent
-- into one project's writing. Passing the band explicitly, including
passing (0.1, 8.0), silences the warning and changes no number.
Source: still standing study (Jensenius), Deichman chest-QoM
cardiorespiratory spectral-composition analysis (deichman_full).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
signal
|
ndarray
|
1-D input signal. |
required |
fs
|
float
|
Sampling rate in Hz. |
required |
bands
|
dict
|
Mapping of band name to |
required |
total_band
|
tuple
|
|
None
|
nperseg_s
|
float
|
Welch segment length in seconds. Defaults to 20. |
20
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
Mapping of each band name to its power fraction in |
Source code in src/micromotion/physio.py
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 | |
respiration_onsets ¶
respiration_onsets(x, fs: float, *, lowpass_hz: float = 1.0, onset_frac: float = 0.55, baseline_hz: float = 0.2) -> dict
Inspiration and expiration onset times from a chest-expansion recording.
Inspiration onset is defined by chest-expansion velocity crossing a threshold rather than by a local minimum, which is what makes this robust on quiet standing: a belt on a standing body carries sway and weight shifts that produce local minima with no breath behind them. Expiration onset is the end of that rise, which assumes passive expiration.
The threshold is taken from the signal's own distribution -- onset_frac times the mean
positive velocity -- and candidate rises are then required to contain an upward crossing of
a heavily low-passed copy of the signal, so a rise that never returns to an exhaled baseline
is discarded.
Zero-crossing technique after Matsuda et al. and Upham (2018). Ported from Finn Upham's
respy (MIT, 2023) and reimplemented on numpy; see :func:respiratory_phases for why.
Returns inspiration_s, expiration_s, the normalised signal, and its velocity.
Source code in src/micromotion/physio.py
268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 | |
respiratory_phases ¶
respiratory_phases(x, fs: float, *, scale_high: float = 0.7, scale_low: float = 0.3, **kw) -> dict
Decompose a respiration recording into the phases of the breath cycle.
A breathing rate says how often; this says where in each cycle the body is. That matters here specifically, because the post-expiration pause is the moment in the cycle when the body is most nearly still, so relating breathing to micromotion wants phases rather than a rate.
Returns boolean masks over the input samples:
inspiration, expiration
the two half-cycles.
inspiration_high, expiration_high
high-flow moments judged against the whole recording -- the scale_high quantile of
velocity within each phase.
inspiration_v, expiration_v
high-flow moments judged within each breath, against that breath's own peak velocity.
Use these when breath size varies across the recording, which it does during settling.
post_expiration
the pause after expiration has slowed to scale_low of its own peak rate.
The defaults come from the coordination analysis in Upham (2018).
Ported from Finn Upham's respy (MIT, 2023) with permission, and reimplemented on numpy.
The port is not gratuitous: respy.Resp_phases assigns through df[col].loc[idx], which
under pandas copy-on-write silently does not write, so on pandas 2 and later it returns all
twelve of its phase columns empty without raising. Verified against respy 0.1.1 on
pandas 3.0.3, where every phase column came back 0.0 per cent populated.
References¶
Upham, F. (2018). Detecting the Adaptation of Listeners' Respiration to Heard Music. PhD thesis, New York University.
Source code in src/micromotion/physio.py
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 | |
Dynamics¶
micromotion.dynamics ¶
Scaling and nonlinearity measures for postural time series.
These describe how movement is structured in time rather than how much of it there is. They
are the part of this package with reuse value beyond the standstill corpus, and they are
also the part where a subtly wrong implementation produces a plausible number rather than an
error. Every function here is covered by a test against a process whose answer is known in
advance; see tests/test_dynamics.py.
One warning carried over from the corpus. These methods read across scales, so they are the
ones that upsampling corrupts. Resample downwards only, and use
:func:micromotion.resample.to_rate, which refuses to do otherwise.
iaaft ¶
iaaft(x, iters: int = 100, rng=None) -> np.ndarray
Iterative amplitude-adjusted Fourier transform surrogate.
Preserves both the amplitude distribution and the power spectrum while destroying any nonlinear structure, so it is the null hypothesis "this series is a linear Gaussian process observed through a monotonic transform".
Source code in src/micromotion/dynamics.py
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | |
phase_surrogate ¶
phase_surrogate(x, rng=None) -> np.ndarray
Phase-randomised surrogate: preserves the spectrum, not the distribution.
Cheaper than :func:iaaft and a weaker null, since a non-Gaussian amplitude
distribution alone can make a linear series look nonlinear against it.
Source code in src/micromotion/dynamics.py
45 46 47 48 49 50 51 52 53 54 55 56 57 58 | |
circular_shift_surrogate ¶
circular_shift_surrogate(x, rng=None, margin: float = 0.1) -> np.ndarray
Rotate a series in time.
Preserves everything about the series itself and destroys only its alignment with another, so it is the right null for a correlation between two recordings and the wrong one for a property of a single series.
Source code in src/micromotion/dynamics.py
61 62 63 64 65 66 67 68 69 70 71 72 | |
surrogate_test ¶
surrogate_test(x, statistic, n: int = 99, method=iaaft, rng=None) -> dict
Compare a statistic against its surrogate distribution.
Returns the observed value, the surrogate mean and standard deviation, a z score and a
two-sided p value. The p value uses the standard (count + 1) / (n + 1) form, which
cannot return zero.
Source code in src/micromotion/dynamics.py
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 | |
trev ¶
trev(x, tau: int = 1) -> float
Time-reversal asymmetry.
A linear Gaussian process looks the same run backwards; most nonlinear ones do not. The
statistic is the skew of the lag-tau increments, normalised by their variance.
Source code in src/micromotion/dynamics.py
94 95 96 97 98 99 100 101 102 | |
dfa ¶
dfa(x, smin: int | None = None, smax: int | None = None, nsc: int = 18, order: int = 1, *, fs: float | None = None, min_scale_s: float | None = None) -> dict
Detrended fluctuation analysis.
Returns the scaling exponent alpha. For reference, 0.5 is white noise, 1.0 is pink
noise, and 1.5 is Brownian motion. Values above 0.5 mean the series persists: a
deviation tends to be followed by more of the same.
THE SCALE FLOOR MATTERS MORE THAN IT LOOKS, and it is the argument most worth setting. The shortest scales of a real recording are dominated by measurement noise, which is white, so including them pulls the exponent toward 0.5. On standstill head-marker series the exponent moves by up to 0.15 between a floor of 8 samples and one of 0.3 s.
GIVE IT IN SECONDS WHERE YOU CAN. smin is a sample count, so the same call measures
different physical scales at different rates: 8 samples is 0.16 s at 50 Hz and 0.08 s at
100 Hz. Pass fs and min_scale_s instead and the floor is the same stretch of time
whatever the recording rate, which is what a comparison across a corpus needs.
A SINGLE EXPONENT MAY NOT DESCRIBE THE SERIES AT ALL. Where the scaling is multifractal the answer depends on which scales are included, and the sensitivity above is a symptom of that rather than of a bad estimator. Check with a multifractal width before quoting one number.
Source code in src/micromotion/dynamics.py
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 | |
mfdfa ¶
mfdfa(x, qs=None, smin: int = 16, smax: int | None = None, nsc: int = 20, order: int = 1) -> dict | None
Multifractal detrended fluctuation analysis.
Returns h(q), the singularity spectrum, the generalised Hurst exponent h2 and
the spectrum width. A width near zero means one scaling exponent describes the whole
series; a wide spectrum means different parts of it scale differently.
A width above about 2 for postural data is a signal to check the preprocessing rather than to celebrate: widths up to 6.6 were once produced entirely by upsampling.
Source code in src/micromotion/dynamics.py
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 | |
sda ¶
sda(x, y=None, fs: float = 25.0, maxlag: float = 10.0, nlags: int = 60) -> dict
Stabilogram diffusion analysis, after Collins and De Luca.
Postural sway behaves like two different processes at two timescales: over short intervals it drifts away from where it was, and over longer ones it is pulled back. The crossover between them is found by fitting two lines to the log-log mean-square displacement and taking their intersection.
Returns short- and long-term Hurst exponents and diffusion coefficients, and the
critical time and displacement. Hs above 0.5 with Hl below it is the normal
pattern: open-loop drift, then closed-loop correction.
Source code in src/micromotion/dynamics.py
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 | |
fit_two_region ¶
fit_two_region(dt, msd) -> dict
Fit two straight lines to a log-log curve and locate their crossover.
Generic: any measure with a scaling break can use it.
Source code in src/micromotion/dynamics.py
260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 | |
sampen ¶
sampen(x, m: int = 2, r: float | None = None) -> float
Sample entropy: how unpredictable the series is.
The negative log probability that two segments which match for m samples still match
for m + 1. Higher is less regular. r defaults to 0.2 standard deviations.
Source code in src/micromotion/dynamics.py
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 | |
ami ¶
ami(x, maxlag: int = 100, bins: int = 32) -> np.ndarray
Average mutual information against lag, for choosing an embedding delay.
Source code in src/micromotion/dynamics.py
308 309 310 311 312 313 314 315 316 317 318 | |
first_ami_minimum ¶
first_ami_minimum(x, maxlag: int = 100, smooth: int = 9) -> int
The conventional embedding delay: the first local minimum of mutual information.
The curve is smoothed first. Estimated from a histogram it is noisy enough that the first local minimum of the raw curve is usually spurious: on a clean sine, whose true answer is a quarter period, the unsmoothed rule returns a delay of 2.
When the curve has no local minimum within maxlag the answer is its smallest value,
which on a monotone curve is the end of the search and not a property of the signal. This
is a bounded search returning its own boundary, and it warns rather than passing a plausible
integer on silently: on synthetic 1/f² series it returns 96 with the default maxlag=100
every time, four short of the boundary only because the smoothing flattens the tail. Raise
maxlag until a minimum appears, or accept that this signal does not have one.
Source code in src/micromotion/dynamics.py
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 | |
embed ¶
embed(x, dim: int, tau: int) -> np.ndarray
Time-delay embedding. x may be one series or several columns.
Source code in src/micromotion/dynamics.py
355 356 357 358 359 360 361 362 363 | |
rqa ¶
rqa(x, dim: int = 3, tau: int | None = None, rr: float = 0.05, lmin: int = 2) -> dict
Recurrence quantification, at fixed recurrence rate.
The threshold is solved per plot so that exactly rr of pairs count as recurrent.
That matters: with a fixed absolute threshold, determinism partly measures how tightly
packed the trajectory is, so two recordings of different amplitude are not comparable.
Pass several columns to get multidimensional recurrence quantification, which is the form used for between-body coupling.
Leave tau unset unless you have a reason. A delay of 1 makes consecutive embedding
vectors share all but one coordinate, so recurrences chain into diagonals that reflect
the embedding rather than the dynamics: white noise embedded at tau=1 reports a
determinism of 0.57, against 0.08 at a properly chosen delay.
Source code in src/micromotion/dynamics.py
366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 | |
plv ¶
plv(a, b, fs: float, band: tuple[float, float] | None = None) -> dict
Phase-locking value between two signals.
Returns the locking strength from 0 to 1 and the preferred phase difference in radians. Band-pass first if the signals are broadband; locking is only meaningful within a band where each signal has a well-defined phase.
Source code in src/micromotion/dynamics.py
426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 | |
apen ¶
apen(x, m: int = 2, r: float | None = None) -> float
Approximate entropy.
The older sibling of :func:sampen, and biased: it counts each template as matching
itself, which pulls the estimate towards regularity, and the bias grows as the series
gets shorter. Provided because the balance literature reports it and comparisons need
it. For new work prefer :func:sampen, which drops the self-match.
Source code in src/micromotion/dynamics.py
445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 | |
dcca ¶
dcca(a, b, scales=None, order: int = 1) -> dict
Detrended cross-correlation between two non-stationary series.
An ordinary correlation between two signals that each wander is dominated by the wandering. This detrends both inside windows of many sizes and correlates what is left, giving a coefficient per timescale.
That per-scale answer is the point: two bodies can be uncorrelated second to second and
correlated over a minute, and a single number cannot express it. Returns rho against
scales, running from -1 to 1.
Source code in src/micromotion/dynamics.py
471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 | |
Postural geometry¶
micromotion.posture ¶
The shape and extent of standing still.
Quantity of motion says how fast a body part moved. These say where it went: how far it strayed, what area it covered, whether it swayed along one line or wandered in all directions.
The pair is more informative than either alone. Two people with the same quantity of motion can occupy regions differing several-fold, because speed and extent are close to independent in this data — how much someone moves predicts only about a quarter of how large a region they occupy.
sway_geometry ¶
sway_geometry(xy) -> dict
Principal axis, anisotropy and dispersion of a two-dimensional trace.
Give it the horizontal coordinates of a marker, or a centre-of-pressure track.
anisotropy runs from 0 to 1: zero when the sway is equally wide in every direction,
approaching one when it collapses onto a line. There are two other definitions of this
quantity in circulation in this project — one where 1 means isotropic and one reporting
only the angle — so the definition is stated rather than assumed. axis_deg is the
direction of greatest sway, and is axial: 10 degrees and 190 degrees describe the same
posture. Test it with :func:micromotion.circular.rayleigh_axial, never with the
ordinary Rayleigh test.
axis_deg IS IN THE RECORDING'S OWN FRAME AND IS NOT COMPARABLE ACROSS RECORDINGS
unless you know that the frames agree. Two things have to hold: the laboratory
coordinate convention, and which way the body was facing. A sway axis is anatomical --
it is front-to-back -- so turning the person turns this number without anything about
their posture changing.
Both failed in one corpus on one afternoon, which is why this paragraph exists. One championship edition is stored with its horizontal axes about 90 degrees from every other, at the same concentration, so a comparison of mean axes across editions measured the laboratories rather than the bodies. And a three-performer collection recorded many sessions with the performers standing in a circle facing each other, so within-person axial concentration came out at R = 0.23 to 0.48 with axes spread over nearly 180 degrees, against 0.75 to 0.85 in a group all facing one way -- and an analysis read that as a weak personal trait when varying facing predicts it exactly.
anisotropy and dispersion are rotation-invariant and carry neither problem.
So is the CONCENTRATION of a set of axes, which is why a clustering result can survive
where a mean angle cannot. If you need the angle itself across recordings, establish
the facing geometry first; if you cannot, use the invariant quantities and say so.
Source code in src/micromotion/posture.py
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 | |
ellipse_area_95 ¶
ellipse_area_95(xy) -> float
Area of the 95 per cent confidence ellipse, in the input units squared.
The standard summary of postural extent in the balance literature, which makes it the number to report when comparing against clinical work.
Source code in src/micromotion/posture.py
73 74 75 76 77 78 79 80 81 82 83 84 | |
heading_persistence ¶
heading_persistence(xy, speed_percentile: float = 20.0) -> dict
Does the trace reverse along a line, wander at random, or loop?
Takes the heading of each step, and averages the cosine of the change in heading from one step to the next. It measures smoothness, not shape: +1 is a trace whose direction barely changes between samples, 0 is a random walk, and -1 is a trace that reverses at every single step.
A back-and-forth sway along one line reads near +1, not -1, which is worth saying because the intuition runs the other way. Such a trace holds its heading for the whole of each excursion and reverses only at the turning points, so the reversals are a handful of steps among hundreds. Quiet standing in the championship corpus reads about 0.95. Reaching -1 takes a zigzag that flips direction at the sampling rate, which is a signature of noise rather than of movement.
straightness is the net displacement over the distance walked to achieve it, so it
is near 0 for someone who stays put however much they move, and near 1 for someone who
walks away in a straight line.
The slowest steps are dropped before averaging, speed_percentile of them by default.
A heading is the direction of a step, and the direction of a step that barely happened
is mostly noise; including them pulls the mean towards the 0 of a random walk.
.. warning::
This descriptor is unusually sensitive to how the series was brought to its sampling
rate, because a change of heading between consecutive samples is a different question
at every sampling interval. In the standstill corpus a bare polyphase resample read
0.15 here for the two 120 Hz editions --- which need a 5:12 conversion, where the
anti-alias filter has least room --- against about 0.95 for the others, and the split
was very nearly published as a difference between tracking systems. Through
:func:~micromotion.to_rate every edition reads 0.949 to 0.964. Bring series to a
common rate with that, and compare only series that share one.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
xy
|
Positions, shape (N, 2). Non-finite rows are dropped. |
required | |
speed_percentile
|
float
|
Percentile of step speeds below which steps are excluded from the heading average. Defaults to 20. |
20.0
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
dict
|
|
Source code in src/micromotion/posture.py
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 | |
path_length ¶
path_length(xy, fs: float | None = None) -> dict
Total distance travelled, and the rate at which it accumulated.
Unfiltered and undifferentiated, so it is not a quantity of motion and is not comparable with one: sensor noise adds to path length monotonically, which means a noisier device reports a longer path for an identical movement. Reported because the balance literature uses it, and because it correlates with head quantity of motion at 0.61 in this corpus, which is worth knowing but is not an equivalence.
Source code in src/micromotion/posture.py
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 | |
dispersion_radius ¶
dispersion_radius(xy, quantile: float = 0.95) -> float
Radius containing a given proportion of the samples, about the mean position.
A robust alternative to the ellipse area when the trace has excursions: a single lean that lasts two seconds changes the ellipse considerably and this hardly at all.
Source code in src/micromotion/posture.py
169 170 171 172 173 174 175 176 177 178 179 180 | |
shared_axis_projection ¶
shared_axis_projection(markers, reference, mask: str = 'reference') -> dict
Project several markers onto ONE axis, the reference marker's principal axis.
This is the difference between asking "does this marker sway along a line" and asking
"do these segments sway along the SAME line". :func:micromotion.principal_axis_projection
answers the first, per marker, giving each its own axis -- and two markers swaying at right
angles, each on its own axis, then correlate perfectly while sharing no direction of motion.
Correlating segments only means something once they are on a common axis, which is what
this builds.
The axis is anatomical -- front-to-back for a standing person -- so it is taken from one
named marker rather than from the pooled cloud, which would be dominated by whichever
segment moved most. axis_deg carries the same warning as
:func:sway_geometry's: it is in the recording's own frame and is not comparable across
recordings unless the laboratory convention and the direction the body faced are both known.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
markers
|
Mapping of name -> horizontal coordinates, each |
required | |
reference
|
Key naming the marker whose axis every other is projected onto. |
required | |
mask
|
str
|
Which samples count. |
'reference'
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
dict
|
|
dict
|
unit vector), |
|
dict
|
mask used when |
|
dict
|
finite samples in that marker's projection). |
Raises:
| Type | Description |
|---|---|
KeyError
|
If |
ValueError
|
If |
Source code in src/micromotion/posture.py
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 | |
segmental_coordination ¶
segmental_coordination(markers, reference, ratios=None, mask: str = 'reference', min_finite: float = 0.8, reduce=None) -> dict
Do a body's segments sway as one rigid link, or as several?
Projects every marker onto the reference marker's axis with
:func:shared_axis_projection, then reduces the result three ways: how strongly each pair
of segments agrees, how many independent axes the set spans, and which of a named pair
sways further. A body rocking at the ankles as a single inverted pendulum gives high
correlations, an effective degrees-of-freedom near one, and a head-to-hip amplitude ratio
above one -- the head is further from the ankle, so the same rotation carries it further.
The effective degrees of freedom is :func:micromotion.effective_dimensionality called on
the projections with rank=False, not a second implementation of the participation
ratio. Ranking is right for heavy-tailed descriptors and wrong here, where the columns are
sway signals and their actual covariance is the quantity of interest.
The correlations are pairwise-complete and each carries its own n. np.corrcoef
returns NaN if a single sample is missing, which drops whole sessions from a pair without
saying so -- one or two per pair in the corpus this comes from, under a heading that
reported one N for every row. Read n beside any correlation.
No body-part taxonomy is built in: markers is whatever the caller names, and
"inverted pendulum" is an interpretation of a ratio above one rather than something this
computes. Group membership stays with the study.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
markers
|
Mapping of name -> horizontal coordinates, each |
required | |
reference
|
Key naming the marker whose axis defines the shared direction. |
required | |
ratios
|
optional
|
Pairs |
None
|
mask
|
str
|
Passed to :func: |
'reference'
|
reduce
|
optional
|
Which markers enter the dimensionality reduction, by name.
Defaults to all of them. Pass an explicit list when |
None
|
min_finite
|
float
|
Markers finite on a smaller fraction of the usable
frames than this -- the reference-valid frames under |
0.8
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
dict
|
|
dict
|
given), |
|
dict
|
requested pairs), |
|
dict
|
entered the reduction) and |
Raises:
| Type | Description |
|---|---|
KeyError
|
If |
ValueError
|
If |
Source code in src/micromotion/posture.py
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 | |
Balance and centre of pressure¶
micromotion.balance ¶
Posturography and standstill-sway metrics for centre-of-pressure (CoP) and head/marker position signals.
This module ports the "still standing" study's posturography stack into pure numpy/scipy surfaces that operate on plain arrays -- no study-specific loaders, axis conventions, or marker loops. Three families of measures are provided:
- Sway amount / geometry -- :func:
cop_sway_metrics, :func:confidence_ellipse_area, :func:convex_hull_area. - Control dynamics / complexity -- :func:
stabilogram_diffusion(Collins-De Luca SDA), :func:dfa(detrended fluctuation analysis), :func:sample_entropy, :func:spectral_edges, :func:sway_texture, :func:principal_axis_projection. - Direction / extent -- :func:
sway_orientation, :func:axial_rayleigh, :func:spatial_extent.
The from-scratch SDA / DFA / sample-entropy implementations are validated in the test-suite against known-answer synthetic signals (white noise -> DFA alpha ~= 0.5 and SDA Hurst ~= 0.5; a sine -> low sample entropy relative to its shuffle).
Source: still standing study (Jensenius) -- posturography and micromotion analyses of the international "standstill" championships and related datasets.
confidence_ellipse_area ¶
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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
xy
|
ndarray
|
Point cloud of shape |
required |
conf
|
float
|
Confidence level in |
0.95
|
Returns:
| Name | Type | Description |
|---|---|---|
float |
Ellipse area in squared position units (e.g. mm^2), or
|
Source code in src/micromotion/balance.py
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 | |
convex_hull_area ¶
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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
xy
|
ndarray
|
Point cloud of shape |
required |
Returns:
| Name | Type | Description |
|---|---|---|
float |
Convex-hull area in squared position units, or |
Source code in src/micromotion/balance.py
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 | |
cop_sway_metrics ¶
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).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
xy
|
ndarray
|
CoP path of shape |
required |
t
|
ndarray
|
Per-sample timestamps in seconds. May be irregular. Defaults to None. |
None
|
fs
|
float
|
Constant sampling rate in Hz, used when |
None
|
freq_band
|
tuple
|
|
(0.1, 5.0)
|
resample_fs
|
float
|
Uniform rate in Hz onto which the path is interpolated before the spectral estimate. Defaults to 50.0. |
50.0
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
Metrics with keys |
Source code in src/micromotion/balance.py
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 | |
principal_axis_projection ¶
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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
xy
|
ndarray
|
Point cloud of shape |
required |
Returns:
| Type | Description |
|---|---|
|
np.ndarray: 1-D projection of shape |
Source code in src/micromotion/balance.py
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 | |
stabilogram_diffusion ¶
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).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
xy
|
ndarray
|
Sway path of shape |
required |
fs
|
float
|
Sampling rate in Hz. |
required |
short_max_s
|
float
|
Upper bound (s) of the short-term fitting window. Defaults to 0.6. |
0.6
|
long_min_s
|
float
|
Lower bound (s) of the long-term fitting window. Defaults to 1.5. |
1.5
|
n_lags
|
int
|
Number of log-spaced lags at which the MSD is evaluated. Defaults to 40. |
40
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
|
Source code in src/micromotion/balance.py
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 | |
dfa ¶
dfa(x, *, n_scales=18, min_scale=10, fs=None, min_scale_s=None)
Detrended fluctuation analysis (DFA) scaling exponent, as a float.
White noise gives alpha ~= 0.5, pink noise 1.0, a random walk 1.5.
This is a thin wrapper over :func:micromotion.dynamics.dfa, which is the single
implementation and additionally returns the scales and fluctuation curve. The float
return is kept because published code calls it that way.
min_scale is a SAMPLE COUNT, so the same call measures different physical scales at
different recording rates. Pass fs and min_scale_s to give the floor in seconds
instead; on a corpus spanning several rates that is what a comparison needs. See
:func:micromotion.dynamics.dfa for how much the floor is worth.
Shared with musicalgestures. The two implementations agree to 1.2 per cent on Brownian motion before being merged; the surviving one sits closer to the analytic answer of 1.5.
Source code in src/micromotion/balance.py
311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 | |
sample_entropy ¶
sample_entropy(x, m=2, r=0.2)
Sample entropy: the negative log probability that two segments matching for m
samples still match for m + 1. Higher means less regular.
r is a tolerance in units of the signal's standard deviation.
A thin wrapper over :func:micromotion.dynamics.sampen, the single implementation.
Shared with musicalgestures; the two agree to 0.2 per cent before
being merged.
Source code in src/micromotion/balance.py
336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 | |
spectral_edges ¶
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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
ndarray
|
1-D input signal. |
required |
fs
|
float
|
Sampling rate in Hz. |
required |
edges
|
tuple
|
Cumulative-power fractions in |
(0.5, 0.95)
|
nperseg
|
int
|
Welch segment length in samples. Defaults
to |
None
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
Mapping |
Source code in src/micromotion/balance.py
354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 | |
sway_texture ¶
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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
speed
|
ndarray
|
1-D speed signal (e.g. mm/s). |
required |
fs
|
float
|
Sampling rate in Hz. |
required |
frozen_threshold
|
float
|
Speed below which the signal is
considered "frozen", in the units of |
2.0
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
|
Source code in src/micromotion/balance.py
394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 | |
sway_orientation ¶
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.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
xy
|
ndarray
|
Point cloud of shape |
required |
Returns:
| Name | Type | Description |
|---|---|---|
dict |
|
Source code in src/micromotion/balance.py
430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 | |
axial_rayleigh ¶
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.
This takes DEGREES. :func:micromotion.circular.rayleigh_axial computes the same quantity
from RADIANS, and the two agree exactly on R and on the mean axis -- a report that they
disagree is a units error at the call site, not a difference of method. Feeding radians to this
function converts already-small numbers a second time, collapsing every angle towards zero and
inflating R towards 1, which is the usual way that mistake shows itself.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
angles_deg
|
ndarray
|
Axial angles in degrees. |
required |
Returns:
| Name | Type | Description |
|---|---|---|
dict |
|
Source code in src/micromotion/balance.py
464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 | |
spatial_extent ¶
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).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
pos
|
ndarray
|
Position trace of shape |
required |
fs
|
float
|
Sampling rate in Hz. |
required |
ellipse_conf
|
float
|
Confidence level for the ellipsoid volume. Defaults to 0.95. |
0.95
|
window_s
|
float
|
Window length in seconds for the within-window dispersion. Defaults to 20.0. |
20.0
|
vertical_axis
|
int
|
Index of the vertical axis ( |
None
|
Returns:
| Name | Type | Description |
|---|---|---|
dict |
|
Source code in src/micromotion/balance.py
509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 | |
Groups¶
micromotion.group ¶
Did these people move at the same moments?
Everything else in this package describes one recording, or relates a pair. This module asks the group question: given many people recorded simultaneously, did their movements coincide in time more than chance allows, and if so, when?
That question needs a null hypothesis with some care in it. Comparing each pair and averaging loses the timing. Shuffling samples destroys each person's own rhythm along with the alignment. The scheme used here shifts each person's event train by an independent random offset drawn from a bounded range: every individual keeps their own event rate and local structure, and only the alignment between people is destroyed. That is the null that "they moved together" should be tested against.
Two related things are provided. :func:event_train turns a continuous signal into a point
process of events. :func:coincidence_test asks whether those events line up across people,
returning both a single score for the whole recording and a p value at each moment, so that the
moments where the group actually converged can be located rather than inferred.
Method after Finn Upham's activity-analysis work, used with permission and reimplemented here;
see finn42/aa_test_package. The stilling-response statistic below is from Upham, Hoffding and
Rosas, "The Stilling Response" (Music & Science, 2024).
event_train ¶
event_train(x, fs: float, frame_s: float = 1.0, threshold: float | None = None, kind: str = 'increase') -> np.ndarray
Turn a continuous signal into a binary point process of events.
An event is a change across a window of frame_s seconds that exceeds threshold.
kind selects "increase", "decrease" or "change" for either direction.
threshold defaults to one standard deviation of the framed differences, which makes the
definition of "an event" relative to how much this person moves rather than absolute. That
matters when the group wears different sensors on different parts of the body: an absolute
threshold would count only the noisiest participants as ever doing anything.
Source code in src/micromotion/group.py
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 | |
coincidence_test ¶
coincidence_test(trains, fs: float, n_surrogates: int = 1000, shift_range_s: float = 30.0, frame_s: float = 1.0, rng=None) -> dict
Test whether events coincide across people more than chance.
trains is a sequence of equal-length binary event trains, one per person, already on a
common timebase.
Returns the observed coincidence count over time, the surrogate mean, a p value at every
sample, a surprise transform of it, and two summaries of the whole recording.
Read frac_significant before score. score is the mean of -log10(p) over every
sample, kept because it is the published statistic, but it is insensitive to exactly the
case these recordings usually present. Measured on twenty simulated people whose event
trains were independent apart from a few shared moments: three shared moments scored 0.32
and sixty scored 0.47, while wholly independent trains scored 0.34. The three-moment case
scored below the null case, because a mean over three thousand samples is dominated by the
99.9 per cent of them where nothing was happening. Over the same runs
frac_significant went 0.012, 0.034, 0.100, 0.200 against 0.004 for independence, which
is the separation you want.
The surrogates shift each person independently by up to shift_range_s seconds. Keep that
range comfortably longer than the timescale you are testing for and shorter than the
recording, or the null starts to resemble the data.
Source code in src/micromotion/group.py
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 | |
participation_ratio ¶
participation_ratio(series, event_times, fs: float, pre=(-3.0, -2.0), post=(0.0, 1.0)) -> np.ndarray
The fraction of people whose movement decreased after each event.
series is (n_people, n_samples) of quantity of motion on a common timebase;
event_times are moments of interest in seconds. For each event, every person's mean in
the pre window is compared with their mean in the post window, and the statistic is
the proportion who went down.
Counting signs rather than sizes is the point. Participants in these recordings wear different sensors in different places, and their quantity of motion differs by more than an order of magnitude, so any statistic that averages magnitudes is dominated by whoever wore the noisiest device. A proportion is immune to that, and to missing data.
Compare the result against :func:sliding_null, not against 0.5 — see that function for
why.
Source code in src/micromotion/group.py
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 | |
sliding_null ¶
sliding_null(series, fs: float, step_s: float = 0.1, **kw) -> np.ndarray
The same statistic computed at every moment of the recording.
This is the null distribution :func:participation_ratio should be judged against, and it
is better than a theoretical one for a reason worth stating: people standing still are not
a coin flip. Movement decays after any excursion, so at a randomly chosen moment rather more
than half of a group is usually already slowing down. Testing observed events against 0.5
would find an effect in any recording; testing them against this finds one only where the
events beat the recording's own baseline.
Compare the two with a one-sided Kolmogorov-Smirnov test.
Source code in src/micromotion/group.py
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 | |
sequential_stability ¶
sequential_stability(x) -> float
How steady a per-cycle measure is, ignoring slow drift.
The median absolute difference between consecutive values. Unlike a standard deviation it is unmoved by a gradual trend, so it answers "was this person's breathing steady" over a window in which their rate is also slowly changing — which is the usual situation.
Source code in src/micromotion/group.py
195 196 197 198 199 200 201 202 203 204 | |
Circular statistics¶
micromotion.circular ¶
Statistics for angles. This module is the family's reference for them.
Where circular statistics live across the four toolboxes. micromotion owns
them: the axial tests sway needs, the circular-linear correlation, the V-test.
ambiscape keeps a small circstats for the time-series end --- phase_stats,
relative_phase --- and the few primitives those need, rather than taking a
dependency on this package for six short functions. MGT re-exports this module
along with the rest of micromotion.
That division decayed once and would have again. On 2026-08-12 the two
Rayleigh implementations were found to disagree on about a fifth of random
cases: this one uses Wilkie's approximation, ambiscape used Zar's earlier
series expansion, both are published and neither was wrong, and nothing said
they were meant to match. ambiscape now matches this module, and
ambiscape/tests/test_circstats_agreement.py asserts it on every run rather
than trusting it. One further difference is deliberate and remains:
circ_corr returns a dict here and a float there, so a caller who swaps the
import gets a different shape --- changing either is an API break, and the
test pins the arithmetic so at least the numbers cannot drift on top of it.
Sway direction, the phase of a breath, the time of day a session happened: all are circular, and ordinary statistics quietly give wrong answers on them. The mean of 350 degrees and 10 degrees is 0, not 180.
Sway direction needs a further distinction. A body swaying forwards and backwards along one line has no preferred direction, only a preferred axis — 10 degrees and 190 degrees are the same posture. Statistics that treat those as opposite will report a person with a strong front-back sway as having no directional preference at all. The axial variants here double the angles before averaging, which is the standard fix, and they are the ones to use for sway.
circ_mean ¶
circ_mean(angles, weights=None) -> dict
Mean direction and concentration.
Returns the mean angle in radians and R, the resultant length, which runs from 0 for
angles spread uniformly around the circle to 1 for angles all identical. R is the
circular answer to "how consistent is this", and there is no separate standard deviation
worth reporting beside it.
Source code in src/micromotion/circular.py
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 | |
circular_sd ¶
circular_sd(R: float) -> float
Circular standard deviation, in radians, from a resultant length.
sqrt(-2 ln R): zero when every angle agrees, growing without bound as
the resultant vanishes. Added 2026-08-12 because this module owns circular
statistics for the family and musiscape needed it --- a consumer that
cannot get a primitive from the owner writes its own, which is how the
Rayleigh implementations drifted apart in the first place.
Source code in src/micromotion/circular.py
57 58 59 60 61 62 63 64 65 66 67 68 69 | |
rayleigh_from_R ¶
rayleigh_from_R(R: float, n: int) -> float
Rayleigh p from a resultant length and a count, without the angles.
:func:rayleigh is the entry point when the angles are to hand. This one
is for callers holding a summary --- a resultant length computed earlier,
or read from a table --- and uses the same Wilkie approximation, so the
two agree by construction rather than by intention.
Source code in src/micromotion/circular.py
72 73 74 75 76 77 78 79 80 81 82 83 | |
rayleigh ¶
rayleigh(angles) -> dict
Test whether angles are spread uniformly around the circle.
The null hypothesis is no preferred direction. A small p rejects it. Uses the standard approximation, which is accurate for n above about 10.
Note what this cannot detect: a perfectly bidirectional distribution, with half the
angles at 0 and half at 180, has a resultant length of zero and passes as uniform. Use
:func:rayleigh_axial when the data are axial, which sway is.
Source code in src/micromotion/circular.py
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 | |
rayleigh_axial ¶
rayleigh_axial(angles) -> dict
Rayleigh test for axial data, where an angle and its opposite are the same.
Doubles the angles, tests those, and halves the resulting mean direction back. This is the correct test for sway direction, marker orientation, and anything else defined on a line rather than a ray.
Source code in src/micromotion/circular.py
107 108 109 110 111 112 113 114 115 116 117 118 119 | |
axial_dispersion ¶
axial_dispersion(angles) -> float
Circular standard deviation of axial data, in radians.
Zero for a body swaying along one fixed line, rising towards the isotropic case.
Source code in src/micromotion/circular.py
122 123 124 125 126 127 128 129 130 131 132 | |
circ_corr ¶
circ_corr(a, b) -> dict
Correlation between two circular variables, after Jammalamadaka and Sengupta.
Source code in src/micromotion/circular.py
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 | |
circ_corr_linear ¶
circ_corr_linear(angles, x) -> dict
Correlation between a circular variable and a linear one.
The question behind "does quantity of motion depend on the time of day", or on the day of the year, where the predictor wraps and the outcome does not.
Source code in src/micromotion/circular.py
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 | |
vtest ¶
vtest(angles, mu: float) -> dict
Test for a preferred direction at a specified angle.
More powerful than :func:rayleigh when there is a prior expectation of where the
preference lies — a hypothesis that sway is front-to-back, for instance, rather than a
search for whichever direction happens to win.
Source code in src/micromotion/circular.py
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 | |
Alignment¶
micromotion.align ¶
Aligning recordings that share no clock.
The recurring problem in this corpus: two instruments recorded the same person at the same time and neither knows what time it was. A phone counts seconds from when its app started, an fNIRS headband from when its acquisition began, an audio recorder from when someone pressed the button. Putting them on one timeline means recovering the offset from the signals themselves.
Two methods, for two situations. When both signals carry the same physiological rhythm, track
that rhythm in each and cross-correlate the resulting curves: :func:instantaneous_rate then
:func:xcorr_lag. When the signals are of different kinds, different lengths, or sampled far
apart, search integer offsets directly with :func:search_lag, which tolerates all three.
Both report how confident they are, and neither should be used without looking. A cross-correlation always has a maximum; that the maximum means anything is a separate claim.
instantaneous_rate ¶
instantaneous_rate(x, fs: float, band: tuple[float, float] = CARDIAC_BAND, win_s: float = 20.0, step_s: float = 1.0, per_minute: bool = True)
Track a rhythm's frequency over time, by sliding a Welch estimate along it.
Returns (times, rate). With per_minute the rate is in beats or breaths per
minute; otherwise in Hz.
This turns a raw signal into something two instruments can be compared on even when they measure entirely different quantities. A haemoglobin trace and a chest accelerometer have no common units, but both carry a heartbeat, and the way that heart rate wanders over ten minutes is a signature specific enough to align them.
Each window is a bare maximum inside the band, and it stays one, because the values are what
published alignments were computed from. What is new is that the function counts the windows
whose band held no peak at all, by :func:~micromotion.spectral.is_band_floor, and warns once
if any did. On synthetic 1/f series with nothing in the cardiac band every window returns the
band's lowest bin, so the track is a flat line sitting on the edge — and two such flat lines
from two instruments will cross-correlate with each other perfectly while carrying nothing.
Where a NaN per window is preferable to a number, loop over
:func:~micromotion.spectral.spectral_peak instead; where a track already exists,
:func:~micromotion.spectral.band_edge_sweep settles whether it is the band edge.
Source code in src/micromotion/align.py
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 | |
xcorr_lag ¶
xcorr_lag(a, b, fs: float = 1.0, max_lag_s: float | None = None, min_r: float = 0.5, difference: bool = True) -> dict
Offset between two equally sampled signals, by cross-correlation.
Returns the lag in seconds and samples, the correlation at that lag, and confident,
which is whether that correlation reaches min_r. A positive lag means b starts
later than a.
That sentence was false until 1.13.0. The function returned the negative of what it
documented, and the test asserted the value the code produced rather than the value the
docstring promised, so the two agreed with each other and neither agreed with the
convention. A sign is exactly the kind of error a test written after the fact preserves.
musicalgestures.xcorr_lag documented and returned the opposite -- the correct --
sign throughout, which is how this was found.
Both series are differenced first, and that default is load-bearing rather than a
stylistic choice. Correlating two drifting series and taking the best lag is the classic
spurious-regression trap: over 200 pairs of independent random walks, the best-lag
correlation had a median of 0.47, exceeded 0.5 forty per cent of the time, and reached
0.98 at worst. No threshold can survive that. Differencing the same pairs gave a median
of 0.11 and never exceeded 0.16. Differencing is shift-equivariant, so the lag itself is
unaffected — it removes the drift, not the alignment. Pass difference=False only if
the inputs are already stationary and you know it.
The confidence test is a threshold on the correlation and deliberately not a measure of how sharply the peak stands out from the rest of the curve. Two such measures were tried and both get it backwards: a drifting series correlates highly with itself at every nearby lag, so its true peak looks unremarkable against the curve, while white noise gives a flat curve whose maximum looks strikingly sharp. A genuine alignment scored 0.98 and pure noise 0.07, but the noise had the higher peak-to-background ratio of the two.
Source code in src/micromotion/align.py
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 | |
search_lag ¶
search_lag(t_a, x_a, t_b, x_b, max_lag_s: float = 300.0, step_s: float = 1.0, min_overlap_s: float = 120.0, min_r: float = 0.5) -> dict
Offset between two irregular or unequal-length signals, by direct search.
Both series are put on a common grid, then every integer offset within max_lag_s is
scored by Pearson correlation over whatever the two share at that offset. Slower than
:func:xcorr_lag and far more tolerant: the series may differ in length, in sampling,
and in how much of the session they cover.
A positive lag means x_b starts later than x_a, the same convention as
:func:xcorr_lag. This function stated no convention at all until 1.13.0 and returned
the negative of that one; both were corrected together, since two alignment functions in
one module disagreeing about a sign is worse than either being wrong alone.
confident is True only when the best correlation reaches min_r and the overlap
reaches min_overlap_s. Recordings that fail the test are better left without an offset
than given one nobody can trust: a wrong alignment is harder to detect downstream than a
missing one.
Source code in src/micromotion/align.py
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 | |
find_transient ¶
find_transient(x, fs: float, threshold: float = 8.0, search_s: float | None = None, min_separation_s: float = 1.0)
Locate impulsive events: a hand clap, a tap on a sensor, a heel strike.
Returns their times in seconds. threshold is in robust standard deviations above the
median of the signal envelope, using the median absolute deviation so that the events
themselves do not inflate the scale they are measured against.
Where sessions open and close with a synchronisation clap, the usual handling is to trim a fixed window from each end. Detecting it instead gives the offset rather than discarding it, which is what turns a clap into an alignment rather than a nuisance.
Source code in src/micromotion/align.py
186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 | |
apply_lag ¶
apply_lag(t, lag_s: float)
Shift a timebase onto another recording's origin.
Source code in src/micromotion/align.py
224 225 226 | |
Feature vector¶
micromotion.features ¶
One canonical feature vector per recording.
Every attempt to compare recordings across this corpus -- clustering, identity classification,
condition classification, the dimensionality reductions in :mod:micromotion.descriptors -- needs
the same thing first: a fixed set of numbers describing one recording. Each attempt had been
inventing its own, which makes two results incomparable for reasons that have nothing to do with
the question either was asking.
This is that set, and it is the only thing of its kind the package offers. Models belong outside the package, in an analysis repository where a train/test split and its leakage are visible; what belongs here is the input they all start from.
Eleven descriptors, in three groups:
- amount and smoothness --
qom,jerk - frequency and texture --
centroid,f50,frozen,burst - sway geometry --
path,extent,area,anis,vert
The geometric five need true position and are nan for accelerometer collections. A chest-worn
accelerometer cannot give a sway ellipse, and doubly integrating one to fake it reports drift as
posture.
feature_vector ¶
feature_vector(x, fs: float, kind: str | None = None, unit: str | None = None, sensor_fs: float | None = None) -> dict | None
The eleven descriptors for one recording, or None if it is too short to describe.
kind and unit must be passed and are not guessed, because the collections do not
record the same quantity: the optical ones store position in mm and the accelerometer ones
store acceleration in g or m/s^2. Differentiating an acceleration series as though it were
position shifts every descriptor two derivatives up and still returns finite,
plausible-looking numbers, which is the failure mode this signature exists to prevent.
sensor_fs is the rate the device actually sampled at, which is not always the rate
the file is stored on: a uniform grid can be an upsample of a much slower sensor, and
checking the grid would admit a band the data cannot carry. Pass it whenever the two differ.
Everything derives from a band-limited velocity, so the filter order and the units are identical whichever way the recording arrived.
Source code in src/micromotion/features.py
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 | |
Equivalence testing¶
micromotion.equivalence ¶
Testing that an effect is absent, rather than failing to show that it is present.
Most of this corpus's interesting results are nulls: no environment effect, no seasonal rhythm, no coupling between performers, no association between sound level and quantity of motion. A non-significant test does not support any of those claims. It says the data are compatible with no effect, and equally compatible with an effect too small for this sample to resolve. Stated as "no effect was found", that is an overclaim, and reviewers say so.
Equivalence testing states the claim the reports actually want to make. Two one-sided tests (TOST) invert the usual logic: the null is that the effect is at least as large as some bound, and rejecting it supports the statement that the effect is smaller than that bound. The bound is a smallest effect size of interest, and choosing it is a scientific judgement, not a statistical one -- which is the point. "No effect" is not a testable claim; "smaller than half a millimetre per second" is.
Report both together. A result that is neither significant nor equivalent is genuinely inconclusive, and saying so is more honest than either alternative.
>>> tost_paired(before, after, bound=0.5) # doctest: +SKIP
{'equivalent': True, 'p': 0.004, ...}
tost_paired ¶
tost_paired(a, b, bound: float, alpha: float = 0.05) -> dict
Two one-sided tests on a paired difference, against a bound in the data's own units.
bound is the smallest difference worth caring about. Pass it in the units of a and
b -- mm/s for quantity of motion, breaths per minute for a rate -- rather than as a
standardised effect size, because a reader can argue about millimetres and cannot argue
about Cohen's d.
Source code in src/micromotion/equivalence.py
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 | |
tost_independent ¶
tost_independent(a, b, bound: float, alpha: float = 0.05) -> dict
The same for two independent samples, using Welch's standard error.
Source code in src/micromotion/equivalence.py
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 | |
equivalence_correlation ¶
equivalence_correlation(r: float, n: int, bound: float, alpha: float = 0.05) -> dict
Is a correlation smaller in magnitude than bound?
For the corpus's many "no association" results, which are reported as correlations rather than as mean differences. Works on Fisher's z, where the sampling distribution is normal with a standard error that does not depend on the correlation itself.
Takes r and n rather than the raw series, so it can be applied to a correlation that
is already published without recomputing it.
Source code in src/micromotion/equivalence.py
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 | |
interpret ¶
interpret(result: dict) -> str
One sentence a report can quote, naming the bound rather than hiding it.
Source code in src/micromotion/equivalence.py
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 | |
Dimensionality and reliability¶
micromotion.descriptors ¶
How many independent things is a descriptor set measuring, and how much of one is the person?
Two reductions that this corpus kept re-implementing per report, with the arithmetic drifting between copies. Both are here so there is one of each.
They answer questions that are easy to confuse. :func:effective_dimensionality asks how many
independent axes a set of measures spans -- whether eleven descriptors are eleven findings or
three wearing different names. :func:intraclass_correlation asks, of a single measure, how
much of its variance is the person rather than the occasion -- whether it is a trait or a state.
A name collision worth knowing about. group.participation_ratio is a different quantity
with a similar name: the fraction of a group whose movement decreased after an event. The
participation ratio of an eigenvalue spectrum, which is what effective dimensionality means
here, is deliberately not called that.
effective_dimensionality ¶
effective_dimensionality(x, rank: bool = True, by=None) -> dict
How many independent dimensions a set of descriptors spans.
x is (n_observations, n_descriptors). Returns the variance share of each component, the
number of components needed for 80 and 90 per cent, and the participation ratio
(sum lambda)^2 / sum(lambda^2) -- an effective count that needs no cutoff and is not
an integer.
Three choices are baked in because getting them wrong is what made earlier versions of this disagree with each other:
rank=True correlates the ranks rather than the values. Descriptors here are heavy-tailed
-- burstiness is a ratio of a 99th percentile to a median -- and on raw values a single
descriptor with a long tail dominates the first component and the answer becomes a statement
about that descriptor's outliers.
by standardises within groups before pooling, and should be the recording session, edition
or collection. Without it, a between-group difference in level appears as shared variance and
inflates the first component: descriptors do not become more correlated because two editions
were recorded at different rates, but they look it.
Columns that are constant or all-NaN are dropped, and the count is reported, because a degenerate column silently adds an eigenvalue of zero and deflates the participation ratio.
Source code in src/micromotion/descriptors.py
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 | |
intraclass_correlation ¶
intraclass_correlation(values, groups, log: bool | None = None) -> dict
The share of a measure's variance that is between groups rather than within them.
Fitted as a random-intercept mixed model, value ~ 1 with a random intercept per group,
which is the estimator the corpus uses for "is this a trait or a state": groups is the
person and the residual is the occasion.
log=None log-transforms when every value is strictly positive, because these are scale
quantities whose residuals are otherwise skewed; pass False to force the raw scale.
Returns the ICC, both variance components, the counts, and boundary. Check
boundary. A random-effect variance can be estimated at exactly zero, which is the
optimiser hitting the edge of the parameter space rather than a measurement of no group
effect -- the difference matters when the number of groups is small, and reporting 0.000
from such a fit implies a precision that is not there.
Source code in src/micromotion/descriptors.py
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 | |