rm3100/characterize.py

#!/usr/bin/env python3
"""Characterise an RM3100 capture: noise floor, spectrum, stability.

    ./.venv/bin/python characterize.py capture_60s.csv

Produces a four-panel figure and a text summary:

  Amplitude spectral density  nT/sqrt(Hz) against the 1.2 nT/sqrt(Hz) the
                              manual quotes (Table 3-1), and against the
                              white-noise level implied by the sample sd.
  Allan deviation             where averaging stops helping and drift takes
                              over -- the honest measure of a noise floor.
  Residual distribution       after removing a slow trend, so a non-Gaussian
                              tail or quantisation shows up.
  Read latency                host-side diagnostics only -- the measurement
                              grid itself is uniform regardless.

Three further figures are written beside it, each answering something the four
panels collapse away:

  _spectrogram.png            whether a feature lasted the run or was a burst.
  _drift.png                  how far the chip's oscillator moved against the
                              host clock, and -- from the two-sample deviation
                              of the rate -- how much of that is real. These
                              parts drift 186-1241 ppm against a floor near
                              10 ppm, so it is a signal by a wide margin.
  _timebase.png               which clock each line is coherent on. A line at a
                              fixed frequency sharpens when the samples are
                              placed by the host clock; one locked to the
                              sampling is wrecked by it, which settles from a
                              single capture what compare.py needs two rates to
                              decide. Diagnostic only: every other figure here
                              uses the chip's grid, because neither time base is
                              right for both kinds of line.

Two noise figures are printed, and the difference between them matters. `sd` is
the plain standard deviation, which on a drifting capture measures the drift:
one run here reads 167.9 nT of sd on an axis whose actual noise is 12.1 nT.
`white sd` is sd(diff)/sqrt(2), a first difference that rejects anything slower
than the sample rate, and is what Table 3-1's 208/sqrt(cycle_count) should be
compared against.

`--trim` drops equal time from both ends. A capture usually starts while the
sensor is still settling and ends while it is being handled for the next run,
and neither belongs in a noise figure. Nothing is dropped unless asked.

The sample grid is exact, not assumed: the sensor samples on its own internal
schedule (manual sections 5.7.2, 5.8.2), logger.py writes a placeholder rather
than skip a grid point, and capture.py refuses a non-contiguous sample index.
Frequencies use the rate measured against the host clock, not the nominal table
value, which is ~6% out on this unit.

Rows lost during capture carry no data and are linearly interpolated on load, so
a capture with many of them will read as artificially smooth at high frequency.
capture.py's summary reports how many, and characterize.py prints it above.
"""

import argparse
import math
import sys
from collections import namedtuple

import matplotlib
matplotlib.use("Agg")
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np

import capture
import rm3100

SURFACE = "#fcfcfb"
TEXT_PRIMARY = "#0b0b0b"
TEXT_SECONDARY = "#52514e"
GRID = "#e3e2df"
REFERENCE = "#8a8880"

# Categorical slots 1-3; three peer axes, validated all-pairs in light mode.
AXES = [("x", "X", "#2a78d6"), ("y", "Y", "#eb6834"), ("z", "Z", "#1baf7a")]

# Total field is derived from the three axes rather than a fourth peer, so it
# takes neutral ink instead of a competing hue -- which also keeps the
# categorical set at the three slots that validate.
TOTAL = ("total", "Total", TEXT_PRIMARY)
SERIES = AXES + [TOTAL]

# Table 3-1: "Noise Density @ Max. Single-Axis Sample Rate".
SPEC_ASD_NT = 1.2

# Every panel is drawn on a fixed scale so two runs can be laid side by side and
# compared by eye. Autoscaling defeats that completely: a quieter capture simply
# redraws its own axis and looks identical to a noisy one. The ranges are wide
# enough for any configuration this driver offers -- cycle counts 50 to 800, so
# rates from ~2 to ~530 Hz and LSBs from 3 to 50 nT. Where an axis is
# logarithmic the range is whole decades, so the tick labels are powers of ten.
#
# The ASD frequency axis is the exception: linear, not logarithmic. A log axis
# gives the sub-hertz drift most of the width and crams every line worth naming
# -- mains, its harmonics, fs/4 -- into the last fifth. Linear spaces them
# evenly, which is how a spectrum is read for lines. It costs the low-frequency
# decades, which the Allan panel beside it covers better anyway, and it cuts off
# above 100 Hz: at cycle count 100 that leaves out 100-150 Hz, where nothing but
# the noise floor has ever appeared.
ASD_HZ = (0.0, 100.0)
ASD_NT = (5e-1, 5e1)            # floor ~1 nT/rtHz, lines to ~15
ALLAN_TAU_S = (1e-3, 1e3)       # one sample period to a 1000 s average
ALLAN_NT = (5e-1, 5e1)
RESIDUAL_NT = (-150.0, 150.0)   # +/-6 LSB at cycle count 100
RESIDUAL_DENSITY = (0.0, 0.045)
LATENCY_MS = (0.0, 20.0)        # covers 1.9 ms at cc=50 to 13.7 ms at cc=400
LATENCY_COUNT = (0.8, 1e5)

# compare.py plots the same two quantities normalised by |B|, so they get their
# own fixed ranges in ppm. Same reasoning: two invocations should overlay. Its
# frequency axis stays logarithmic: that figure overlays captures whose rates
# differ several-fold, and log is what puts them on comparable footing.
FRACTIONAL_ASD_PPM = (1e0, 1e3)
FRACTIONAL_ALLAN_PPM = (1e0, 1e3)
FRACTIONAL_ASD_HZ = (1e-2, 1e3)

# Spectrogram segment length. Longer than SEGMENT_SECONDS because the job here
# is to resolve one line from another -- 40 s gives 0.025 Hz, enough to separate
# a mains line from a sample-locked one even when they land close together --
# and a ~300 s capture still yields ~100 columns at this overlap.
SPECTROGRAM_SECONDS = 40.0
# 31/32. Overlapping segments are not independent, so this buys no new
# information -- what it buys is columns: the hop is 1.25 s instead of 5 s, so a
# feature that lasts a few seconds is drawn as a few columns rather than falling
# between two. It costs only render time, and the axes are wide enough to show it.
SPECTROGRAM_OVERLAP = 0.96875
# Spectrograms are drawn in ppm of |B| per sqrt(Hz), not nT/sqrt(Hz), for the
# same reason compare.py's figures are: a capture whose gain is 7% lower records
# 7% fewer nanotesla of the same noise, and a colour scale shared across captures
# would show that as a quieter sensor.
#
# Sharing a colour scale also requires the *rendering* to be identical, which is
# easier to get wrong. Two things are pinned below to make it so:
#
#   SPECTROGRAM_BIN_HZ     one displayed frequency step, reached by averaging
#                          native bins. Since the window is a fixed *duration*,
#                          native resolution is the same for every capture, so a
#                          fixed step means the same number of bins averaged and
#                          therefore the same degrees of freedom per cell. This
#                          is the one that matters: dof is what sets how widely
#                          cells scatter, and so where the colour distribution
#                          sits. Before it, a 150 Hz capture had 6.6x more bins
#                          averaged into each pixel than a 38 Hz one, its speckle
#                          was suppressed, and that read as a lower noise level
#                          when nothing differed.
#
# A single fixed frequency axis would remove the last of it -- the ~1.8x
# difference in how many binned rows still fall in one pixel -- but only by
# drawing a 38 Hz capture into the bottom quarter of a 160 Hz axis and leaving
# three quarters blank. `--max-freq` does that for anyone who wants it; the
# default is each capture's own Nyquist, which is worth the residual.
#
# The range spans the floor (~30 ppm/rtHz here) to a line an order of magnitude
# above it.
SPECTROGRAM_PPM = (15.0, 200.0)
SPECTROGRAM_BIN_HZ = 0.1        # 4 native bins at a 40 s window: 8 dof per cell

# Turbo by default: a spectrogram is read for narrow lines against a broad floor,
# and hue steps make those jump out in a way a single hue cannot. Turbo is the
# right rainbow to reach for -- it is jet without the sharp lightness kinks that
# fake banding.
#
# It is worth being accurate about what it is not. Measured over its own 256
# steps, turbo's L* runs 12 to 91 but is *not* monotonic: it climbs to a bright
# yellow mid-scale then darkens into red, and 129 of 255 steps go down. Jet is
# 108 of 255 with a worse worst-case step, so turbo is smoother, not ordered.
# Neither survives colour vision deficiency well. Two alternatives are a flag
# away: `viridis`, which is genuinely monotonic (0 of 255 steps decrease), and
# `sequential`, the single-hue ramp below (L* 97 to 17, hue 213 +/- 1 deg).
SPECTROGRAM_COLORMAP = "turbo"
SEQUENTIAL_STEPS = ["#f4f8fd", "#cfe0f5", "#93bbe9", "#4e8ddb",
                    "#2a78d6", "#1b5091", "#0e2a4d"]

# Lower edge of the band every broadband figure is quoted over. Below a few Hz
# the spectrum is drift, not noise floor, and it varies far more between runs
# than the floor does -- so including it would compare environments rather than
# sensors.
BAND_LO_HZ = 3.0
# Upper edge, as a fraction of Nyquist. This has to stay clear of the corner of
# any filter a compared series has been through, or the comparison scores one
# path partway down a rolloff and reads it as a quieter sensor. The binding one
# is decimation: DECIMATE_CUTOFF_FRACTION is where that filter turns over, so
# BAND_NYQUIST_FRACTION / 2 is held below it with margin to spare.
BAND_NYQUIST_FRACTION = 0.8

# Welch segment duration. Frequency resolution is 1/this, so 20 s of segment
# resolves 0.05 Hz -- fine against a band starting at 3 Hz -- and a 5-minute
# capture still yields ~30 overlapping segments to average.
SEGMENT_SECONDS = 20.0

# A sample-locked tone is reported at this many standard errors. Under the null
# the statistic is Rayleigh(1), so 5 is a one-in-270,000 bin and stays clear of
# the ~20 fractions tested. Real artefacts here reach 30.
LINE_SIGMA = 5.0

# Drift figure. The left panel fits the rate on windows of DRIFT_WINDOW_S; the
# right sweeps that window length to show what it costs. 10 s is a compromise
# and not the optimum -- sigma_y bottoms out nearer 2 s on these parts -- but 2 s
# windows put 145 points on a 290 s capture, which draws as a band rather than a
# trend. At 10 s there are ~29, the floor is still only 13-34 ppm, and the drift
# being plotted is 186-1241 ppm, so the error bars stay a minor part of the span.
DRIFT_WINDOW_S = 10.0
# Log-spaced, from a window that holds a few hundred samples to a third of a
# 300 s capture -- past that there are fewer than three windows and `rate_allan`
# drops the point rather than estimate a deviation from one difference.
RATE_TAUS = (1.0, 1.5, 2.0, 3.0, 5.0, 7.0, 10.0, 15.0,
             20.0, 30.0, 50.0, 70.0, 100.0)
DRIFT_PPM = (-800.0, 800.0)     # the worst run here spans 1241 ppm peak to peak
RATE_TAU_S = (1.0, 1e2)
RATE_ALLAN_PPM = (1e0, 1e3)

# Time-base figure. The host clock is fitted against sample index by a
# polynomial of this degree: high enough to follow a warm-up transient, low
# enough that it cannot chase the per-read scheduling jitter, which carries no
# information about the chip. Degree 1 is exactly the uniform grid already in
# use, and nothing above 3 changed any measured amplitude on these captures --
# consistent with `calibrate.rate_allan`, which finds nothing above the noise
# below a ~2 s timescale.
TIMEBASE_DEGREE = 3
# Half-width of the scan around each line, and how finely it is sampled. A
# 290 s capture resolves 1/T = 0.0034 Hz, so +/-0.05 Hz is ~30 resolution widths
# -- wide enough to show the shape of a peak and whether it moved, narrow enough
# that the peak is not a single pixel.
TIMEBASE_SPAN_HZ = 0.05
TIMEBASE_POINTS = 401
TIMEBASE_PANELS = 3
# A verdict needs a peak that is actually a peak: this is how far the scan's
# maximum must stand above the scan's own median before the ratio of the two
# time bases means anything.
#
# Local contrast rather than significance against the broadband noise, because
# that is the question. A line can sit at a respectable sigma across the whole
# band and still have no resolvable peak in a 0.1 Hz window -- the 60 Hz alias
# in |B| at cycle count 400 does exactly that, ranking third at 9 sigma while
# its scan is visibly flat -- and comparing where two flat curves happen to peak
# is comparing noise. 3x is a little above what a smooth curve reaches by
# chance and well below the 4-6x a real line here shows.
TIMEBASE_MIN_CONTRAST = 3.0


def band_for(fs):
    """The (lo, hi) analysis band for a capture sampled at fs."""
    return BAND_LO_HZ, BAND_NYQUIST_FRACTION * fs / 2


def welch_asd(v, fs, nperseg=4096):
    """Amplitude spectral density in units/sqrt(Hz) via Welch's method."""
    nperseg = min(nperseg, len(v) // 4 * 2 or len(v))
    step = nperseg // 2
    window = np.hanning(nperseg)
    # Normalisation for a one-sided PSD with this window.
    scale = 1.0 / (fs * (window ** 2).sum())

    segments = []
    for start in range(0, len(v) - nperseg + 1, step):
        seg = v[start:start + nperseg]
        # Linear detrend: removes DC and any slow ramp that would smear
        # energy across the low-frequency bins.
        seg = seg - np.polyval(np.polyfit(np.arange(nperseg), seg, 1),
                               np.arange(nperseg))
        spectrum = np.abs(np.fft.rfft(seg * window)) ** 2 * scale
        spectrum[1:-1] *= 2.0          # fold negative frequencies
        segments.append(spectrum)

    psd = np.mean(segments, axis=0)
    freqs = np.fft.rfftfreq(nperseg, 1.0 / fs)
    return freqs[1:], np.sqrt(psd[1:])   # drop DC bin


def allan_deviation(v, fs, points=40):
    """Overlapping Allan deviation of the signal against averaging time tau."""
    n = len(v)
    max_m = n // 4
    ms = np.unique(np.geomspace(1, max(max_m, 2), points).astype(int))
    taus, devs = [], []
    cumulative = np.concatenate([[0.0], np.cumsum(v)])

    for m in ms:
        # Bin means of length m, taken at every offset (overlapping).
        means = (cumulative[m:] - cumulative[:-m]) / m
        diffs = means[m:] - means[:-m]
        if diffs.size < 2:
            continue
        taus.append(m / fs)
        devs.append(np.sqrt(0.5 * np.mean(diffs ** 2)))
    return np.array(taus), np.array(devs)


def white_sd(v):
    """Per-sample noise with drift removed, in the units of v.

    The first difference of white noise has twice its variance, so dividing by
    sqrt(2) recovers the original sd -- while any component slower than the
    sample rate is differenced away. On a capture that drifts, this and the
    plain sd measure entirely different things, and only this one is the noise.

    It is not free: differencing has a +6 dB/octave response, so a narrowband
    line near Nyquist is weighted about twice as heavily as the broadband floor.
    Use band_stats for a figure that is flat across the band.
    """
    if len(v) < 2:
        return 0.0
    return float(np.diff(v).std() / np.sqrt(2))


def segment_length(fs, seconds=SEGMENT_SECONDS):
    """Welch segment length in samples for a fixed segment *duration*.

    Holding the duration fixed rather than the sample count is what makes two
    captures at different rates comparable. A fixed nperseg gives them different
    frequency resolutions and, worse, different numbers of segments to average:
    the median of a Welch estimate is biased low by an amount that depends on
    that count, so a slower capture would read as a quieter one by a few percent
    for no physical reason at all.
    """
    return max(64, int(round(fs * seconds)))


def band_stats(v, fs, band=None, nperseg=None):
    """Broadband figures over a stated band: (median ASD, RMS, peak, peak Hz).

    The median is the floor -- robust to the handful of bins a line occupies --
    while the RMS integrates everything in the band including those lines, so
    the two together say how much of the band is line rather than floor.
    """
    lo, hi = band if band is not None else band_for(fs)
    freqs, asd = welch_asd(v, fs, nperseg=nperseg or segment_length(fs))
    keep = (freqs >= lo) & (freqs <= hi)
    if not keep.any():
        nan = float("nan")
        return nan, nan, nan, nan
    f, a = freqs[keep], asd[keep]
    return (float(np.median(a)), float(np.sqrt(np.trapezoid(a ** 2, f))),
            float(a.max()), float(f[np.argmax(a)]))


def fir_lowpass(cutoff, ntaps):
    """Linear-phase low-pass FIR: a sinc truncated by a Blackman window.

    `cutoff` is in cycles per sample (0 to 0.5). Written out rather than taken
    from scipy, which is not a dependency here and would be a large one for one
    filter. Blackman buys a -74 dB stopband for a transition band of about
    5.5/ntaps, which is the trade that matters: the whole point of the filter is
    that what it stops does not fold back.

    Normalised to unit gain at DC, so decimating does not rescale the field.
    """
    if not 0 < cutoff < 0.5:
        raise ValueError(f"cutoff {cutoff} must be in (0, 0.5) cycles/sample")
    ntaps = int(ntaps) | 1          # odd, so the delay is a whole sample
    if ntaps < 3:
        raise ValueError(f"ntaps {ntaps} is too short to filter anything")
    n = np.arange(ntaps) - (ntaps - 1) / 2
    h = 2 * cutoff * np.sinc(2 * cutoff * n) * np.blackman(ntaps)
    return h / h.sum()


# Where the anti-alias filter turns over, as a fraction of the decimated rate.
# It must sit above BAND_NYQUIST_FRACTION / 2 so that nothing compared is
# scored on the rolloff, and below 0.5 so that nothing folds back. The
# transition band left over is (0.5 - this)/k, and Blackman gives about
# 5.5/ntaps, so the tap count has to grow with k: 128*k+1 leaves roughly 2x
# margin on that and is odd by construction.
DECIMATE_CUTOFF_FRACTION = 0.45
assert BAND_NYQUIST_FRACTION / 2 < DECIMATE_CUTOFF_FRACTION < 0.5


def decimate(v, k, method="fir"):
    """Downsample by an integer k, low-passing first so nothing folds in.

    Two filters, because they answer different questions:

      boxcar  average k consecutive samples. This is what the sensor itself
              does over a cycle count, so it is the right comparison for "would
              a slower cycle count have given me this?" -- and it shares the
              boxcar's poor stopband, which is why the sensor aliases.
      fir     a proper anti-alias filter. This is what you would actually do in
              post, and it is the one that removes out-of-band interference
              instead of folding it.

    Returns the decimated series; the new rate is fs/k.
    """
    k = int(k)
    if k < 1:
        raise ValueError(f"decimation factor {k} must be at least 1")
    if k == 1:
        return np.asarray(v, dtype=float)
    v = np.asarray(v, dtype=float)
    if method == "boxcar":
        n = len(v) // k * k
        if n < k:
            raise ValueError(f"{len(v)} samples cannot be decimated by {k}")
        return v[:n].reshape(-1, k).mean(axis=1)
    if method != "fir":
        raise ValueError(f"unknown decimation method {method!r}")
    h = fir_lowpass(DECIMATE_CUTOFF_FRACTION / k, 128 * k + 1)
    if len(v) <= 2 * len(h):
        raise ValueError(
            f"{len(v)} samples is too few for a {len(h)}-tap filter; "
            f"decimating by {k} needs at least {2 * len(h) + 1}")
    # Discard a full filter length at each end rather than half: 'same' pads
    # with zeros, so the taper reaches ntaps//2 in and a half-length trim would
    # leave the tail of it in the data.
    return np.convolve(v, h, mode="same")[len(h):-len(h)][::k]


SampleLine = namedtuple("SampleLine", "numerator period amplitude sigma")


def sample_locked_lines(v, max_period=8, threshold=LINE_SIGMA):
    """Tones sitting at a simple rational fraction of the sample rate.

    Returns SampleLine(numerator, period, amplitude, sigma) for each surviving
    fraction numerator/period of fs, strongest first. `amplitude` is the peak
    amplitude of the tone in the units of v.

    What this finds that a spectrum cannot distinguish is *what a line is locked
    to*. A line at a fixed frequency lands on a different fraction of fs when
    the rate changes; one that stays at fs/4 across captures whose rates differ
    is locked to the sampling, so it is an artefact of the measurement rather
    than a field. Comparing two captures settles it -- see compare.py.

    Evaluating the transform at exactly j/period rather than reading a spectrum
    puts the whole capture behind one number, which is what makes a 1 nT tone
    detectable at 30 sigma under 20 nT of broadband noise.

    Only fractions in lowest terms are tested, so each frequency is reported
    once. Without that, a period-4 pattern would also be reported at 6 and 8,
    which are multiples of it and carry no additional information.

    A cubic trend is removed first, so drift cannot leak into a low fraction.
    """
    v = np.asarray(v, dtype=float)
    if len(v) < 4 * max_period:
        return []
    index = np.arange(len(v))
    v = v - np.polyval(np.polyfit(index, v, 3), index)
    sd = v.std()
    if sd == 0:
        return []
    found = []
    for period in range(2, max_period + 1):
        for numerator in range(1, period // 2 + 1):
            if math.gcd(numerator, period) != 1:
                continue
            if 2 * numerator == period:
                # Nyquist has no phase to fit -- it is a real alternation, so
                # one degree of freedom rather than two, and a tighter null.
                amplitude = abs(float((v * (-1.0) ** index).mean()))
                sigma = amplitude / (sd / np.sqrt(len(v)))
            else:
                phase = np.exp(-2j * np.pi * numerator / period * index)
                amplitude = 2 * abs(complex((v * phase).mean()))
                # Real and imaginary parts each have variance 2 sd^2 / len, so
                # the magnitude is Rayleigh(1) in these units under the null:
                # 5 sigma is a one-in-270,000 fluctuation.
                sigma = amplitude / (sd * np.sqrt(2 / len(v)))
            if sigma >= threshold:
                found.append(SampleLine(numerator, period, float(amplitude),
                                        float(sigma)))
    return sorted(found, key=lambda line: -line.sigma)


def spectrogram(v, fs, seconds=SPECTROGRAM_SECONDS, overlap=SPECTROGRAM_OVERLAP):
    """(times, freqs, ASD) -- how the spectrum of a series evolves.

    The Welch average that `band_stats` reports is this, collapsed along time.
    Keeping the time axis is what distinguishes a line that was there all along
    from a burst that a single average would smear into a raised floor, and it
    is the only view that shows the difference at a glance.

    Segments are long and heavily overlapped: the point here is frequency
    resolution fine enough to separate a mains line from a sample-locked one, and
    at these capture lengths there are plenty of samples to spend on it.
    """
    nperseg = max(64, int(round(fs * seconds)))
    if len(v) < nperseg:
        raise ValueError(f"{len(v)} samples is shorter than one "
                         f"{seconds:g} s window ({nperseg} samples)")
    step = max(1, int(round(nperseg * (1.0 - overlap))))
    window = np.hanning(nperseg)
    scale = 1.0 / (fs * (window ** 2).sum())
    index = np.arange(nperseg)
    starts = range(0, len(v) - nperseg + 1, step)
    columns, times = [], []
    for start in starts:
        seg = v[start:start + nperseg]
        seg = seg - np.polyval(np.polyfit(index, seg, 1), index)
        spectrum = np.abs(np.fft.rfft(seg * window)) ** 2 * scale
        spectrum[1:-1] *= 2.0
        columns.append(np.sqrt(spectrum[1:]))
        times.append((start + nperseg / 2) / fs)
    freqs = np.fft.rfftfreq(nperseg, 1.0 / fs)[1:]
    return np.array(times), freqs, np.array(columns).T


def bin_frequency(freqs, asd, step=SPECTROGRAM_BIN_HZ):
    """Average a spectrogram onto a fixed frequency step. Returns (freqs, asd).

    Averaging in *power* -- the array is already amplitude, so it is squared and
    rooted around the mean -- because that is what adds degrees of freedom. The
    point is not to save pixels but to give every capture the same dof per cell:
    the window is a fixed duration, so native resolution is identical across
    captures, and a fixed step therefore averages an identical number of bins.
    """
    width = max(1, int(round(step / (freqs[1] - freqs[0]))))
    if width == 1:
        return freqs, asd
    usable = len(freqs) // width * width
    grouped = (asd[:usable] ** 2).reshape(-1, width, asd.shape[1])
    return freqs[:usable].reshape(-1, width).mean(axis=1), np.sqrt(grouped.mean(axis=1))


def alias_of(frequency, fs):
    """Where `frequency` lands after sampling at fs, and whether it folded."""
    folded = abs(frequency - round(frequency / fs) * fs)
    return folded, frequency > fs / 2


def dither_check(v, lsb, longest=1024):
    """Does averaging still recover resolution below one quantiser step?

    Returns [(n, sd, sd_in_lsb, ratio_to_ideal), ...] for block averages of n
    samples. Ideal is sd(1)/sqrt(n), so a ratio near 1 means averaging is
    buying everything it should.

    This is the check the whole cycle-count choice rests on. Section 5.1 warns
    that quantisation "generally dictates working above a cycle count of ~30",
    and the failure it warns about is specific: if the sensor's own noise is
    small against the LSB, samples stop straddling the boundary, the quantiser
    stops being dithered, and averaging stalls at a fixed fraction of an LSB no
    matter how long you average. What that looks like here is a ratio that
    climbs while the LSB column stops falling.

    A slow trend is removed first. Otherwise drift -- which averaging cannot
    remove and is not supposed to -- is what stops the average shrinking, and
    the answer would be about the site rather than the quantiser.
    """
    v = np.asarray(v, dtype=float)
    index = np.arange(len(v))
    v = v - np.polyval(np.polyfit(index, v, 5), index)
    base = white_sd(v)
    if base == 0:
        return []
    out = []
    n = 1
    while n <= longest and len(v) // n >= 16:
        usable = len(v) // n * n
        means = v[:usable].reshape(-1, n).mean(axis=1)
        sd = float(means.std())
        out.append((n, sd, sd / lsb, sd / (base / np.sqrt(n))))
        n *= 4
    return out


def decade_ticks(ax, which="both", minimum_decades=1.0):
    """Label a log axis at powers of ten, with unlabelled minors between.

    Matplotlib falls back to labelling minor ticks when a log axis spans less
    than a decade, which produces a scale reading 2x10^0, 3x10^0, 4x10^0 --
    dense, and not what a log axis is read as. Decades are the convention, so
    the range is widened to hold at least one rather than the labelling being
    changed to suit a narrow range.
    """
    axes = {"x": [ax.xaxis], "y": [ax.yaxis],
            "both": [ax.xaxis, ax.yaxis]}[which]
    for axis in axes:
        setter = ax.set_xlim if axis is ax.xaxis else ax.set_ylim
        lo, hi = ax.get_xlim() if axis is ax.xaxis else ax.get_ylim()
        if lo > 0 and hi > lo:
            short = minimum_decades - math.log10(hi / lo)
            if short > 0:
                pad = 10 ** (short / 2)
                setter(lo / pad, hi * pad)
        axis.set_major_locator(mticker.LogLocator(base=10.0))
        axis.set_minor_locator(
            mticker.LogLocator(base=10.0, subs=tuple(range(2, 10))))
        axis.set_minor_formatter(mticker.NullFormatter())


def reference_lines(fs):
    """Frequencies worth marking on a spectrum: [(hz, label, is_alias), ...].

    Mains is drawn where it actually lands, not where it is generated. Above
    Nyquist it has folded, and the folded position is the only place it can be
    seen -- so that is what gets the line, labelled as an alias so it is never
    mistaken for a real signal at that frequency.
    """
    marks = [(fs / 4, "fs/4", False), (fs / 2, "fs/2 = Nyquist", False)]
    for mains in (60.0, 120.0):
        landed, folded = alias_of(mains, fs)
        if landed < fs / 2 * 0.995:
            marks.append((landed,
                          f"{mains:g} Hz" + (" alias" if folded else ""),
                          folded))
    return sorted(marks)


def coherent_amplitude(v, t, freqs):
    """Amplitude of a tone at each of `freqs`, for samples taken at times `t`.

    Evaluates 2*|mean(v * exp(-2*pi*i*f*t))| directly. That is what a spectrum
    computes, but without requiring the samples to be evenly spaced, which is
    the entire reason this exists rather than a call to `welch_asd`: the
    question here is which *time base* a line is coherent on, and one of the two
    candidates is not uniform.

    A tone present at f throughout the record sums in phase and returns its own
    amplitude; anything else averages towards zero as 1/sqrt(n). A sample placed
    at the wrong time by a fraction d of a cycle contributes cos(2*pi*d) instead
    of 1, so coherence is lost at a rate set by the tone's frequency -- which
    makes this the sharpest test available of whether a time base is right, and
    also the reason a low-frequency line can say nothing about one.

    The mean is removed first, or a DC offset leaks into the lowest frequencies.
    """
    v = np.asarray(v, dtype=float)
    v = v - v.mean()
    t = np.asarray(t, dtype=float)
    return np.array([2.0 * abs(complex((v * np.exp(-2j * np.pi * f * t)).mean()))
                     for f in np.atleast_1d(freqs)])


def host_time_base(cap, degree=TIMEBASE_DEGREE):
    """Sample times taken from the host clock, smoothed against sample index.

    The alternative to `cap.elapsed`, which is `sample_index * dt_true` and so
    asserts the rate never moved. It did: `calibrate.rate_stability` measures
    186 to 1241 ppm of movement across the captures here.

    Raw host timestamps carry that drift, but also the scheduling jitter of each
    individual read, which says nothing about when the chip actually sampled.
    A low-order polynomial keeps the first and discards the second. Degree 1
    reduces exactly to `cap.elapsed`, which makes the uniform grid the degree-1
    member of this family rather than a separate idea, and nothing above degree
    3 changed any measured amplitude on these captures.

    Returned relative to the first sample, like `cap.elapsed - cap.elapsed[0]`.
    Shifting to that origin is what makes the degree-1 case *exactly* the
    uniform grid rather than merely parallel to it: the fit carries a non-zero
    intercept, since a least-squares line through a curve does not pass through
    its first point.
    """
    index = cap.sample_index - cap.sample_index[0]
    host = cap.system_time - cap.system_time[0]
    if len(index) <= degree + 1:
        return index * cap.dt_true
    fitted = np.polyval(np.polyfit(index, host, degree), index)
    return fitted - fitted[0]


def sibling_path(path, suffix):
    """Name a companion figure from the main one: `x_noise.png` -> `x_drift.png`.

    The trailing `_noise` is stripped first, so the extra figures do not come
    out named `_noise_drift` and `_noise_spectrogram`.
    """
    stem = path.rsplit(".", 1)[0]
    if stem.endswith("_noise"):
        stem = stem[:-len("_noise")]
    return f"{stem}_{suffix}.png"


def make_spectrogram(cap, path, seconds=SPECTROGRAM_SECONDS,
                     overlap=SPECTROGRAM_OVERLAP,
                     colormap=SPECTROGRAM_COLORMAP, max_freq=None):
    """Write a time-frequency plot of |B|, marking the lines worth watching."""
    fs = cap.true_rate_hz
    times, freqs, asd = spectrogram(cap.total, fs, seconds, overlap)
    native = freqs[1] - freqs[0]
    import calibrate
    rate_drift = calibrate.rate_stability(cap)
    freqs, asd = bin_frequency(freqs, asd)
    # Fractional, so a gain difference between captures is not read as a noise
    # difference when they share a colour scale.
    field = float(np.linalg.norm([cap.x.mean(), cap.y.mean(), cap.z.mean()]))
    asd = asd / field * 1e6
    dof = 2 * max(1, int(round(SPECTROGRAM_BIN_HZ / native)))
    top = max_freq or fs / 2

    fig, ax = plt.subplots(figsize=(13.5, 7), dpi=150)
    fig.patch.set_facecolor(SURFACE)
    ax.set_facecolor(SURFACE)
    ramp = (mcolors.LinearSegmentedColormap.from_list("sequential",
                                                      SEQUENTIAL_STEPS)
            if colormap == "sequential" else colormap)
    # imshow rather than pcolormesh: the grid is regular in both axes, and at
    # this resolution that is thousands of times fewer objects to draw.
    # 'antialiased' resamples when there are more rows than pixels. 'nearest'
    # would drop a one-bin line entirely depending on where it fell, which is
    # exactly the feature this plot exists to show.
    image = ax.imshow(asd, origin="lower", aspect="auto", cmap=ramp,
                      norm=mcolors.LogNorm(*SPECTROGRAM_PPM),
                      extent=[times[0], times[-1], freqs[0], freqs[-1]],
                      interpolation="antialiased")
    bar = fig.colorbar(image, ax=ax, pad=0.02, extend="both", fraction=0.04)
    bar.set_label("ppm of |B| per √Hz", color=TEXT_SECONDARY, fontsize=10)
    bar.ax.tick_params(colors=TEXT_SECONDARY, labelsize=9, length=0)
    bar.outline.set_visible(False)

    # Labels sit inside the axes on a surface-coloured chip: outside they get
    # clipped by the colorbar, and a spectrogram has no margin to give away.
    # Anything at the very top is the axis boundary itself -- fs/2 when the plot
    # runs to Nyquist -- where a line adds nothing and its label lands on the
    # title.
    marks = [m for m in reference_lines(fs) if freqs[0] <= m[0] <= top * 0.97]
    for slot, (hz, label, folded) in enumerate(marks):
        ax.axhline(hz, color=TEXT_PRIMARY, linewidth=1.0,
                   linestyle=":" if folded else "--", alpha=0.85)
        # Lines can land within a label's height of each other -- fs/4 and a
        # mains alias were 2.4 Hz apart here. Staggering across the width keeps
        # both readable without moving either line.
        ax.annotate(f"{label}  {hz:.2f} Hz",
                    xy=(0.008 + 0.17 * (slot % 3), hz),
                    xycoords=("axes fraction", "data"),
                    xytext=(0, 7), textcoords="offset points",
                    va="bottom", ha="left", color=TEXT_PRIMARY, fontsize=9,
                    bbox=dict(boxstyle="round,pad=0.25", facecolor=SURFACE,
                              edgecolor="none", alpha=0.85))

    ax.set_title("Spectral density of |B| over time", loc="left",
                 color=TEXT_PRIMARY, fontsize=13, fontweight="bold", pad=8)
    ax.set_xlabel("elapsed (s)", color=TEXT_SECONDARY, fontsize=10)
    ax.set_ylabel("frequency (Hz)", color=TEXT_SECONDARY, fontsize=10)
    ax.set_ylim(0.0, top)
    if fs / 2 < top:
        # Above Nyquist there is no data, and the blank says so -- that missing
        # bandwidth is a real difference between configurations, not a gap.
        ax.annotate(f"no data above Nyquist, {fs / 2:.1f} Hz",
                    xy=(0.5, (fs / 2 + top) / 2), xycoords=("axes fraction", "data"),
                    ha="center", va="center", color=TEXT_SECONDARY, fontsize=10)
    ax.tick_params(colors=TEXT_SECONDARY, labelsize=9, length=0)
    for side in ("top", "right"):
        ax.spines[side].set_visible(False)
    for side in ("left", "bottom"):
        ax.spines[side].set_color(GRID)

    # Independent windows in the whole record. The columns far outnumber these,
    # because they overlap -- so cell-to-cell scatter is the two degrees of
    # freedom of a single periodogram, not a changing spectrum. What is real is
    # what stays put across many columns.
    independent = max(1, int(cap.duration / seconds))
    fig.text(0.5, 0.965,
             f"{cap.path.rsplit('/', 1)[-1]} — {len(cap.sample_index):,} samples "
             f"at {fs:.2f} Hz, cycle count {cap.cycle_count}",
             color=TEXT_PRIMARY, fontsize=10, ha="center")
    fig.text(0.5, 0.935,
             f"oscillator drifted {rate_drift * 1e6:.0f} ppm across the run — "
             f"a line at fs/4 smears by {rate_drift * fs / 4:.4f} Hz.  "
             f"{seconds:g} s windows, {native:.3f} Hz native binned to "
             f"{SPECTROGRAM_BIN_HZ:g} Hz ({dof} dof/cell), {overlap:.2%} overlap, "
             f"{len(times)} columns from {independent} independent windows. "
             f"Dashed = real, dotted = alias."
             + "  Colour scale and dof per cell are fixed across captures.",
             color=TEXT_SECONDARY, fontsize=9, ha="center")
    fig.tight_layout(rect=[0, 0, 1, 0.925])
    fig.savefig(path, facecolor=SURFACE)
    plt.close(fig)
    return path


def make_drift(cap, path, window=DRIFT_WINDOW_S):
    """Write the oscillator's drift against the host clock, with its error bar.

    Two views of one measurement. The left panel is the chip's rate through the
    run; the right is how well that rate can be known as a function of how long
    it is measured for -- which is what decides whether the left panel shows a
    signal or its own noise. On these parts it is emphatically a signal: 186 to
    1241 ppm of drift against a floor near 10 ppm.
    """
    import calibrate
    centres, rates = calibrate.window_rates(cap, window)
    if len(rates) < 3:
        raise ValueError(f"{cap.duration:.0f} s gives {len(rates)} window(s) of "
                         f"{window:g} s, too few to show a trend")
    taus, devs = calibrate.rate_allan(cap, RATE_TAUS)
    if not len(taus):
        raise ValueError(f"{cap.duration:.0f} s is too short for any averaging "
                         f"time in {RATE_TAUS[0]:g}-{RATE_TAUS[-1]:g} s")
    # The error bar is the two-sample deviation at the window length actually
    # plotted. Not the scatter of the points themselves, which is the drift this
    # figure exists to show, and using it would beg the question.
    floor = float(np.interp(window, taus, devs))
    best = int(np.argmin(devs))
    ppm = (rates / rates.mean() - 1.0) * 1e6
    span = calibrate.rate_stability(cap) * 1e6

    fig, axs = plt.subplots(1, 2, figsize=(13.5, 5.5), dpi=150)
    fig.patch.set_facecolor(SURFACE)
    for ax in axs:
        ax.set_facecolor(SURFACE)

    a = axs[0]
    a.axhline(0.0, color=REFERENCE, linestyle="--", linewidth=1.0)
    a.errorbar(centres, ppm, yerr=floor * 1e6, fmt="o", markersize=3.5,
               color=AXES[0][2], ecolor=REFERENCE, elinewidth=1.0, capsize=2.5,
               label=f"{window:g} s windows, ± σ_y = {floor * 1e6:.0f} ppm")
    trend = np.polyfit(centres, ppm, 1)
    a.plot(centres, np.polyval(trend, centres), color=TEXT_PRIMARY, linewidth=1.5,
           label=f"trend {trend[0] * cap.duration:+,.0f} ppm across the run")
    a.set_title("Chip oscillator against the host clock", loc="left",
                color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8)
    a.set_xlabel("elapsed (s)")
    a.set_ylabel("rate error (ppm of the run mean)")
    # The y axis is fixed so two runs can be compared; the x axis is the
    # capture's own length, which pinning would simply clip on a longer one.
    a.set_ylim(*DRIFT_PPM)
    a.annotate("the sample grid stays uniform in index —\nthis is its spacing "
               "in seconds moving.\nGain does not follow: it is set by the "
               "ratio\nof two on-die clocks, which drift together.",
               xy=(0.02, 0.05), xycoords="axes fraction",
               color=TEXT_SECONDARY, fontsize=9)

    a = axs[1]
    a.loglog(taus, devs * 1e6, "o-", color=AXES[0][2], linewidth=1.5,
             markersize=4.5, label="σ_y(τ), measured")
    a.plot([taus[best]], [devs[best] * 1e6], "o", markersize=11,
           markerfacecolor="none", markeredgecolor=TEXT_PRIMARY,
           markeredgewidth=1.4, zorder=5)
    a.annotate(f"{devs[best] * 1e6:.1f} ppm at τ = {taus[best]:g} s",
               xy=(taus[best], devs[best] * 1e6), xytext=(0, -26),
               textcoords="offset points", ha="center",
               color=TEXT_PRIMARY, fontsize=9)
    # Anchored at the longest tau and drawn back to the minimum: white noise
    # cannot rise, so a branch parallel to this is a frequency ramp, and one
    # steeper than it is not.
    rising = taus[best:]
    a.plot(rising, devs[-1] * 1e6 * rising / taus[-1], linestyle=":",
           color=REFERENCE, linewidth=1.5,
           label="∝ τ — a deterministic frequency ramp")
    a.set_title("How well the rate can be measured", loc="left",
                color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8)
    a.set_xlabel("averaging time τ (s)"); a.set_ylabel("σ_y (ppm)")
    a.set_xlim(*RATE_TAU_S); a.set_ylim(*RATE_ALLAN_PPM)
    decade_ticks(a)
    a.annotate("falling = measurement noise averaging down;\nrising = drift the "
               "measurement has resolved",
               xy=(0.02, 0.05), xycoords="axes fraction",
               color=TEXT_SECONDARY, fontsize=9)

    for ax in axs:
        ax.grid(True, which="both", color=GRID, linewidth=0.7)
        ax.set_axisbelow(True)
        for side in ("top", "right"):
            ax.spines[side].set_visible(False)
        for side in ("left", "bottom"):
            ax.spines[side].set_color(GRID)
        ax.tick_params(colors=TEXT_SECONDARY, labelsize=9, length=0)
        ax.xaxis.label.set_color(TEXT_SECONDARY)
        ax.yaxis.label.set_color(TEXT_SECONDARY)
        ax.legend(frameon=False, fontsize=9, labelcolor=TEXT_SECONDARY,
                  loc="upper right")

    fig.text(0.5, 0.955,
             f"{cap.path.rsplit('/', 1)[-1]} — {len(cap.sample_index):,} samples "
             f"at {cap.true_rate_hz:.2f} Hz, cycle count {cap.cycle_count}",
             color=TEXT_PRIMARY, fontsize=10, ha="center")
    # Two figures for the drift, deliberately: the peak-to-peak spread is what
    # `rate_stability` feeds to the uncertainty propagation and what the rest of
    # the tooling quotes, while the fitted trend uses every window and is the
    # better estimate of the warm-up itself. They differ by how much of the
    # transient a coarse window averages over, and printing only one invites the
    # other to be read off the legend and mistaken for a discrepancy.
    fig.text(0.5, 0.915,
             f"drifted {span:.0f} ppm peak to peak, "
             f"{abs(trend[0] * cap.duration):.0f} ppm by the fitted trend, "
             f"against a floor of {devs[best] * 1e6:.1f} ppm at τ = "
             f"{taus[best]:g} s — resolved {span / (devs[best] * 1e6):.0f}× over.",
             color=TEXT_SECONDARY, fontsize=9, ha="center")
    fig.tight_layout(rect=[0, 0, 1, 0.90])
    fig.savefig(path, facecolor=SURFACE)
    plt.close(fig)
    return path


def make_time_base(cap, path):
    """Write, line by line, which clock each one is coherent on.

    A line at a fixed frequency lives in wall time, so it sharpens when the
    samples are placed by the host clock. A line locked to the sampling lives in
    sample index, and the same substitution wrecks it. Measuring both ways
    therefore separates a real signal from an artefact of measuring *within a
    single capture*, where compare.py needs two at different rates to do it.

    The catch is that the test only has force on a strong line: coherence is
    lost in proportion to the tone's frequency, so a low-frequency or aliased
    line says little, and a weak one returns noise. Both are marked as such
    rather than being given a verdict they cannot support.
    """
    fs = cap.true_rate_hz
    chip_t = (cap.sample_index - cap.sample_index[0]) * cap.dt_true
    host_t = host_time_base(cap)
    # The same cubic detrend `sample_locked_lines` applies, and for the same
    # reason: drift is not a line, and left in place it leaks into everything
    # below a few tenths of a hertz.
    order = np.arange(len(cap.total))
    v = cap.total - np.polyval(np.polyfit(order, cap.total, 3), order)
    # Against the quantiser rather than against zero. A capture whose |B| never
    # moves by even a thousandth of an LSB has no spectrum -- every sample is
    # the same integer -- but the cubic detrend still leaves ~1e-11 nT of
    # floating-point residue, so an exact test would let it through and then
    # find "lines" in rounding error.
    if v.std() < cap.lsb_nt * 1e-3:
        raise ValueError("|B| is constant to within a thousandth of an LSB; "
                         "no lines to test")
    null = v.std() * np.sqrt(2 / len(v))

    candidates = [(line.numerator / line.period * fs,
                   f"{line.numerator}/{line.period} of fs", False)
                  for line in sample_locked_lines(v)]
    candidates += [(hz, label, folded)
                   for hz, label, folded in reference_lines(fs)
                   if not label.startswith("fs/")]

    scored = []
    for hz, label, folded in candidates:
        if not 0 < hz < fs / 2:
            continue
        sigma = float(coherent_amplitude(v, chip_t, hz)[0]) / null
        scored.append((sigma, hz, label, folded))
    scored.sort(reverse=True)
    scored = scored[:TIMEBASE_PANELS]
    if not scored:
        raise ValueError("no lines to compare time bases on")

    # Floored rather than purely per-panel: a capture with one or two lines
    # worth testing would otherwise be narrower than its own caption.
    fig, axs = plt.subplots(1, len(scored),
                            figsize=(max(9.8, 4.7 * len(scored)), 5.2),
                            dpi=150, squeeze=False)
    fig.patch.set_facecolor(SURFACE)
    for ax, (sigma, hz, label, folded) in zip(axs[0], scored):
        ax.set_facecolor(SURFACE)
        freqs = np.linspace(hz - TIMEBASE_SPAN_HZ, hz + TIMEBASE_SPAN_HZ,
                            TIMEBASE_POINTS)
        chip = coherent_amplitude(v, chip_t, freqs)
        host = coherent_amplitude(v, host_t, freqs)
        ratio = host.max() / chip.max()
        contrast = chip.max() / np.median(chip)
        if folded:
            verdict = "aliased — folding assumes uniformity, not a clean test"
        elif contrast < TIMEBASE_MIN_CONTRAST:
            verdict = f"no resolved peak ({contrast:.1f}× its floor) — no verdict"
        elif ratio > 1.05:
            verdict = "coherent in wall time — external"
        elif ratio < 0.95:
            verdict = "coherent in sample index — an artefact of measuring"
        else:
            verdict = "no preference between the two"

        ax.plot(freqs - hz, chip, color=AXES[0][2], linewidth=1.3,
                label=f"chip index grid — {chip.max():.3f} nT")
        ax.plot(freqs - hz, host, color=AXES[1][2], linewidth=1.3,
                label=f"host clock, degree {TIMEBASE_DEGREE} — {host.max():.3f} nT")
        # Padded so the verdict can sit between title and axes without either
        # landing on the traces, which reach the top of a tightly scaled panel.
        ax.set_title(f"{label} at {hz:.3f} Hz   {ratio:.2f}×", loc="left",
                     color=TEXT_PRIMARY, fontsize=11, fontweight="bold", pad=22)
        ax.annotate(verdict, xy=(0.0, 1.01), xycoords="axes fraction",
                    color=TEXT_SECONDARY, fontsize=8.5, va="bottom")
        ax.set_xlabel(f"offset from {hz:.3f} Hz (Hz)")
        ax.set_ylabel("coherent amplitude (nT)")
        ax.set_xlim(-TIMEBASE_SPAN_HZ, TIMEBASE_SPAN_HZ)
        # Headroom for the legend, which would otherwise sit on the peak the
        # panel exists to show.
        ax.set_ylim(0.0, max(chip.max(), host.max()) * 1.45)
        ax.grid(True, which="both", color=GRID, linewidth=0.7)
        ax.set_axisbelow(True)
        for side in ("top", "right"):
            ax.spines[side].set_visible(False)
        for side in ("left", "bottom"):
            ax.spines[side].set_color(GRID)
        ax.tick_params(colors=TEXT_SECONDARY, labelsize=9, length=0)
        ax.xaxis.label.set_color(TEXT_SECONDARY)
        ax.yaxis.label.set_color(TEXT_SECONDARY)
        ax.legend(frameon=False, fontsize=8.5, labelcolor=TEXT_SECONDARY,
                  loc="upper right")

    fig.suptitle("Which clock is each line coherent on?", color=TEXT_PRIMARY,
                 fontsize=13, fontweight="bold", y=0.99)
    fig.text(0.5, 0.930,
             f"{cap.path.rsplit('/', 1)[-1]} at {fs:.2f} Hz, "
             f"cycle count {cap.cycle_count}",
             color=TEXT_SECONDARY, fontsize=9, ha="center")
    fig.text(0.5, 0.897,
             "Ratio is host peak over chip peak.  Diagnostic only — no other "
             "figure uses the host time base.",
             color=TEXT_SECONDARY, fontsize=9, ha="center")
    fig.tight_layout(rect=[0, 0, 1, 0.88])
    fig.savefig(path, facecolor=SURFACE)
    plt.close(fig)
    return path


def trimmed(cap, seconds):
    """Drop `seconds` from both ends of a capture. Returns (capture, note).

    A capture typically opens while the sensor is still settling after power-up
    and closes while it is being handled for whatever comes next, and neither
    end is a noise measurement. Trimming is never silent -- the note says what
    went -- and never automatic, because discarding data the caller did not ask
    to discard is how a figure ends up quietly describing a subset.

    Refuses rather than trims when there would be too little left to mean
    anything, since a short deliberate capture is a legitimate thing to analyse.
    """
    if not seconds:
        return cap, ""
    if seconds < 0:
        raise capture.CaptureError(f"--trim {seconds} must not be negative")
    if cap.duration < 4 * seconds:
        return cap, (f"not trimming: {cap.duration:.1f} s is under 4x the "
                     f"{seconds:g} s requested from each end")
    kept = cap.restrict(seconds, cap.duration - seconds)
    return kept, (f"trimmed {seconds:g} s from each end: "
                  f"{len(cap.sample_index):,} -> {len(kept.sample_index):,} "
                  f"samples, {cap.duration:.1f} -> {kept.duration:.1f} s")


def main():
    ap = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("csv")
    ap.add_argument("-o", "--output", default=None)
    ap.add_argument("--start", type=float, default=0.0,
                    help="ignore samples before this elapsed time (s)")
    ap.add_argument("--end", type=float, default=None,
                    help="ignore samples after this elapsed time (s)")
    ap.add_argument("--trim", type=float, default=0.0, metavar="SECONDS",
                    help="drop this many seconds from BOTH ends -- settling at "
                         "the start, handling at the end (default: %(default)s)")
    ap.add_argument("--spectrogram", type=float, default=SPECTROGRAM_SECONDS,
                    metavar="SECONDS",
                    help="spectrogram window length; frequency resolution is "
                         "its reciprocal, time resolution is it "
                         "(default: %(default)s)")
    ap.add_argument("--overlap", type=float, default=SPECTROGRAM_OVERLAP,
                    metavar="FRACTION",
                    help="spectrogram segment overlap, 0 to <1; higher gives "
                         "more columns, not more information "
                         "(default: %(default)s)")
    ap.add_argument("--max-freq", type=float, default=None, metavar="HZ",
                    help="spectrogram frequency ceiling; the default is each "
                         "capture's own Nyquist. Set it the same for every "
                         "capture to make them pixel-for-pixel comparable, at "
                         "the cost of blank space on the slower ones")
    ap.add_argument("--colormap", default=SPECTROGRAM_COLORMAP,
                    help="spectrogram colormap; any matplotlib name, or "
                         "'sequential' for the single-hue ramp "
                         "(default: %(default)s)")
    args = ap.parse_args()

    # --trim and --start/--end both choose a window, and silently letting one
    # win would make the reported figures depend on argument order.
    if args.trim and (args.start or args.end is not None):
        ap.error("--trim sets the window from both ends; it cannot be combined "
                 "with --start/--end, which set it explicitly")

    try:
        cap = capture.load(args.csv)
        if args.start or args.end is not None:
            cap = cap.restrict(args.start or None, args.end)
        whole = cap.duration
        cap, trim_note = trimmed(cap, args.trim)
        trim_applied = cap.duration < whole
    except (OSError, capture.CaptureError) as exc:
        sys.exit(str(exc))

    # Total field is derived from the three axes; its noise is dominated by
    # whichever axis carries most of the field, weighted by direction cosines,
    # so it tracks that axis rather than being an independent measurement.
    data = dict(cap.axes())
    t = cap.elapsed
    # Rate measured against the host clock, not the nominal table value: a 6%
    # error would put every spectral feature 6% off.
    fs = cap.true_rate_hz
    lsb_nt = cap.lsb_nt
    # Host read latency, kept apart from the sample grid it does not affect.
    latency = np.diff(cap.system_time)  # local read spacing, not the fit residual

    # Normalising by |B| makes noise comparable across captures whose scale
    # differs -- a gain change carries the noise with it, so absolute sd alone
    # will read as a noise difference when only the scale moved.
    field = data["total"].mean()
    band = band_for(fs)
    # Table 3-1's noise figure for this cycle count -- what `white sd` is the
    # measurement of, so the two belong in the same table.
    spec_nt = rm3100.expected_noise_nt(cap.cycle_count)

    print(cap.summary())
    if trim_note:
        print(f"  {trim_note}")
    # How far the oscillator moved during the run. Imported here rather than at
    # module scope: calibrate.py imports this module, so the dependency only
    # runs one way at import time.
    import calibrate
    drift = calibrate.rate_stability(cap)
    # Reported, not flagged. The sensor oscillator is an LR relaxation circuit
    # -- coil inductance, a resistor and comparator thresholds, none of them
    # temperature-stable -- so a tenth of a percent across a run is ordinary,
    # not a fault. What it is, is the error bar on this capture's rate and on
    # anything derived from it.
    print(f"  oscillator drifted {drift * 1e6:.0f} ppm across the run "
          f"({drift * cap.dt_true * 1e6:.2f} us on a "
          f"{cap.dt_true * 1e3:.2f} ms period) -- the error bar on the rate")
    print()
    print(f"{'series':6s} {'sd':>9s} {'white sd':>9s} {'dither':>7s} "
          f"{'sd/|B|':>10s} {'p2p':>9s} {'vs spec':>8s} {'ASD':>7s}")

    fig, axs = plt.subplots(2, 2, figsize=(13.5, 9), dpi=150)
    fig.patch.set_facecolor(SURFACE)
    for ax in axs.flat:
        ax.set_facecolor(SURFACE)

    for key, label, color in SERIES:
        v = data[key]
        sd = v.std()
        # Drift-free: what Table 3-1 quotes, and what the dither margin is
        # really made of. The plain sd above it is drift on a drifting capture.
        wsd = white_sd(v)

        freqs, asd = welch_asd(v, fs)
        # Thin: four traces overlaid on a linear 0-100 Hz axis, and the detail
        # between the lines is what the panel is read for. The Allan panel below
        # keeps a heavier stroke -- four smooth curves, not four dense ones.
        axs[0, 0].semilogy(freqs, asd, color=color, linewidth=0.7,
                           label=label, alpha=0.85)

        taus, devs = allan_deviation(v, fs)
        axs[0, 1].loglog(taus, devs, color=color, linewidth=1.6, label=label)

        # Detrend before the histogram so slow drift does not masquerade as
        # a fat tail. Bin edges are locked to quarter-LSB steps: the residual
        # of an integer-valued signal is a comb at one LSB by construction, and
        # arbitrary bin edges beat against it into ragged spikes that look like
        # structure. Aligned bins render the comb as what it is.
        resid = v - np.polyval(np.polyfit(t, v, 3), t)
        step = lsb_nt / 4
        edge = np.ceil(max(abs(resid.min()), abs(resid.max())) / step) * step
        axs[1, 0].hist(resid, bins=np.arange(-edge, edge + step, step),
                       histtype="step", linewidth=1.4,
                       color=color, label=label, density=True)

        median_asd, _, _, _ = band_stats(v, fs, band)
        print(f"{label:6s} {sd:9.1f} {wsd:9.2f} {wsd/lsb_nt:7.2f} "
              f"{sd/field*1e6:9.1f}p {v.max()-v.min():9.1f} "
              f"{(wsd/spec_nt - 1) * 100:+7.1f}% {median_asd:7.2f}")

    print(f"       sd, white sd, p2p in nT; dither = white sd in LSB "
          f"(1 LSB = {lsb_nt:.2f} nT);")
    print(f"       sd/|B| in ppm of {field:,.0f} nT; vs spec against Table 3-1's "
          f"{spec_nt:.1f} nT at cycle count {cap.cycle_count};")
    print(f"       ASD = median nT/rtHz over {band[0]:g}-{band[1]:.1f} Hz.")

    # The assumption the cycle-count choice rests on, checked rather than
    # asserted: with dither under one LSB, does averaging still get below it?
    print("\ndither: averaging against the quantiser")
    for key, label, _ in SERIES:
        steps = dither_check(data[key], lsb_nt)
        if not steps:
            continue
        longest = steps[-1]
        print(f"  {label:6s} {steps[0][2]:.2f} LSB per sample -> "
              f"{longest[2]:.3f} LSB after {longest[0]:,} "
              f"({longest[0] / fs:.1f} s), {longest[3]:.1f}x the ideal 1/sqrt(n)")
    print("  Sub-LSB resolution is being recovered wherever the LSB column "
          "keeps falling.\n  A ratio climbing while it stalls is the "
          "quantiser losing its dither (section 5.1);\n  a ratio climbing "
          "while it falls is just drift, which averaging cannot remove.")

    # Structure locked to the sample index rather than to a frequency. Reported
    # per axis because it is not isotropic -- one axis here carries 5x another.
    lines = [(label, sample_locked_lines(data[key]))
             for key, label, _ in SERIES]
    if any(found for _, found in lines):
        print("\nlines at simple fractions of the sample rate:")
        for label, found in lines:
            for line in found:
                fraction = line.numerator / line.period
                print(f"  {label:6s} {line.numerator}/{line.period} of fs "
                      f"= {fraction * fs:8.3f} Hz   {line.amplitude:6.2f} nT "
                      f"amplitude  {line.sigma:5.1f} sigma")
        print("  These may be sample-locked -- moving with the rate rather than "
              "sitting at a\n  fixed frequency -- which would make them "
              "artefacts of measuring, not field.\n  compare.py decides it, by "
              "checking whether they hold the same fraction of fs\n  in a "
              "capture taken at a different rate.")

    a = axs[0, 0]
    # The band the quoted median comes from, so the number and the picture
    # cannot drift apart. Drawn under the traces, hence the low zorder.
    a.axvspan(band[0], band[1], color=GRID, alpha=0.55, zorder=0, linewidth=0)
    # Named in the legend rather than annotated on the line. The spec sits low
    # in the range, where on a linear axis the traces are densest, and floating
    # text disappears into them wherever it is put.
    a.axhline(SPEC_ASD_NT, color=REFERENCE, linestyle="--", linewidth=1.2,
              label=f"Table 3-1 spec, {SPEC_ASD_NT} nT/√Hz")
    a.set_title("Amplitude spectral density", loc="left",
                color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8)
    a.annotate(f"shaded: {band[0]:g}–{band[1]:.0f} Hz, the quoted band",
               xy=(0.02, 0.04), xycoords="axes fraction",
               color=TEXT_SECONDARY, fontsize=9)
    a.set_xlabel("frequency (Hz)"); a.set_ylabel("nT/√Hz")
    a.set_xlim(*ASD_HZ); a.set_ylim(*ASD_NT)
    decade_ticks(a, which="y")

    a = axs[0, 1]
    a.set_title("Allan deviation", loc="left",
                color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8)
    a.set_xlabel("averaging time τ (s)"); a.set_ylabel("σ (nT)")
    a.set_xlim(*ALLAN_TAU_S); a.set_ylim(*ALLAN_NT)
    decade_ticks(a)
    a.annotate("slope −½ = white noise; upturn = drift",
               xy=(0.02, 0.04), xycoords="axes fraction",
               color=TEXT_SECONDARY, fontsize=9)

    a = axs[1, 0]
    a.set_title("Residual distribution (cubic trend removed)", loc="left",
                color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8)
    a.set_xlabel("nT"); a.set_ylabel("density")
    a.set_xlim(*RESIDUAL_NT); a.set_ylim(*RESIDUAL_DENSITY)

    a = axs[1, 1]
    # Fixed bin edges as well as fixed limits: comparing two histograms whose
    # bins are different widths compares the binning as much as the data.
    a.hist(latency * 1000, bins=np.linspace(*LATENCY_MS, 121), color=REFERENCE)
    a.set_yscale("log")
    a.axvline(cap.dt_true * 1000, color=TEXT_PRIMARY, linestyle="--", linewidth=1.2)
    a.annotate(f"grid {cap.dt_true * 1000:.2f} ms", xy=(cap.dt_true * 1000, 1),
               xytext=(6, 0), textcoords="offset points",
               color=TEXT_PRIMARY, fontsize=9)
    a.set_title("Host read latency", loc="left",
                color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8)
    # The oscillator's own drift used to be annotated here, on a panel about
    # host-side comms that has nothing to do with it. It has its own figure now.
    a.annotate("spread is comms only;\nthe sample grid is exact",
               xy=(0.98, 0.94), xycoords="axes fraction", ha="right", va="top",
               color=TEXT_SECONDARY, fontsize=9)
    a.set_xlabel("interval between reads (ms)"); a.set_ylabel("count")
    a.set_xlim(*LATENCY_MS); a.set_ylim(*LATENCY_COUNT)
    decade_ticks(a, which="y")

    for ax in axs.flat:
        ax.grid(True, which="both", color=GRID, linewidth=0.7)
        ax.set_axisbelow(True)
        for side in ("top", "right"):
            ax.spines[side].set_visible(False)
        for side in ("left", "bottom"):
            ax.spines[side].set_color(GRID)
        ax.tick_params(colors=TEXT_SECONDARY, labelsize=9, length=0)
        ax.xaxis.label.set_color(TEXT_SECONDARY)
        ax.yaxis.label.set_color(TEXT_SECONDARY)
    # Explicit placement: "best" puts the ASD legend on top of the spec-line
    # annotation in the lower left.
    for ax, loc in ((axs[0, 0], "upper right"), (axs[0, 1], "upper right"),
                    (axs[1, 0], "upper right")):
        ax.legend(frameon=False, fontsize=9, labelcolor=TEXT_SECONDARY, loc=loc)

    fig.suptitle("RM3100 noise characterisation", color=TEXT_PRIMARY,
                 fontsize=15, fontweight="bold", y=0.985)
    fig.text(0.5, 0.945,
             f"{len(t):,} samples, {cap.duration:.1f} s at {fs:.2f} Hz measured "
             f"(nominal {cap.nominal_rate_hz:g} Hz, {cap.rate_error * 100:+.1f}%), "
             f"cycle count {cap.cycle_count}, 1 LSB = {lsb_nt:.2f} nT. "
             f"Drift {drift * 1e6:.0f} ppm — see the drift figure."
             + (f" {args.trim:g} s trimmed from each end." if trim_applied
                else ""),
             color=TEXT_SECONDARY, fontsize=10, ha="center")

    fig.tight_layout(rect=[0, 0, 1, 0.935])
    out = args.output or args.csv.rsplit(".", 1)[0] + "_noise.png"
    fig.savefig(out, facecolor=SURFACE)
    print(f"\n-> {out}")

    # Separate files rather than more panels. The spectrogram needs the width to
    # resolve anything; the other two answer questions the four panels do not --
    # whether a feature lasted, how far the chip's clock moved, and which clock
    # each line is coherent on. Each degrades to a message rather than aborting
    # the run, since a short capture can legitimately support none of them.
    for suffix, build in (("spectrogram",
                           lambda p: make_spectrogram(cap, p, args.spectrogram,
                                                      args.overlap,
                                                      args.colormap,
                                                      args.max_freq)),
                          ("drift", lambda p: make_drift(cap, p)),
                          ("timebase", lambda p: make_time_base(cap, p))):
        try:
            print(f"-> {build(sibling_path(out, suffix))}")
        except ValueError as exc:
            print(f"no {suffix}: {exc}")


if __name__ == "__main__":
    main()