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.

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.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 -- and are deliberately whole
# decades so the tick labels are powers of ten.
ASD_HZ = (1e-2, 1e3)            # 0.05 Hz is the lowest bin a 20 s segment gives
ASD_NT = (1e-1, 1e2)            # floor ~1 nT/rtHz, lines to ~15
ALLAN_TAU_S = (1e-3, 1e3)       # one sample period to a 1000 s average
ALLAN_NT = (1e-1, 1e2)
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)

# 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


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 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 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)")
    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}")
    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)
        axs[0, 0].loglog(freqs, asd, color=color, linewidth=1.2,
                         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)
    a.axhline(SPEC_ASD_NT, color=REFERENCE, linestyle="--", linewidth=1.2)
    a.annotate(f"Table 3-1 spec {SPEC_ASD_NT} nT/√Hz", xy=(freqs[1], SPEC_ASD_NT),
               xytext=(0, 5), textcoords="offset points",
               color=REFERENCE, fontsize=9)
    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)

    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.annotate(f"bins are {lsb_nt / 4:.2f} nT, a quarter LSB.\n"
               f"Integer counts make this a comb at {lsb_nt:.1f} nT however\n"
               f"well dithered; how sharp it looks only tracks how far\n"
               f"the axis drifted. Whether averaging still beats the\n"
               f"quantiser is the dither figure printed above.",
               xy=(0.02, 0.97), xycoords="axes fraction", va="top",
               color=TEXT_SECONDARY, fontsize=9)
    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)
    a.annotate("spread is host-side 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"Sample grid is exact."
             + (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}")


if __name__ == "__main__":
    main()