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.
  Sample interval             whether the timing supports spectral analysis
                              at all.

Timing caveat: in continuous measurement mode the sensor samples on its own
internal schedule, so the true sample instants are near-uniform even when our
reads are jittery. The spectral estimates assume uniform spacing at the mean
observed rate. That assumption holds only if no samples were missed or read
twice -- which is exactly what the sample-interval panel is there to check.
"""

import argparse
import csv
import sys

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

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")]

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


def load(path):
    t, cols, gains = [], {a: [] for a, _, _ in AXES}, []
    with open(path, newline="") as fh:
        for row in csv.DictReader(fh):
            t.append(float(row["elapsed_s"]))
            for a, _, _ in AXES:
                ut = float(row[f"{a}_uT"])
                cols[a].append(ut)
                if abs(ut) > 1.0:
                    gains.append(float(row[f"{a}_raw"]) / ut)
    if len(t) < 64:
        sys.exit(f"{path} has too few samples to characterise")
    return (np.array(t),
            {a: np.array(v) * 1000.0 for a, v in cols.items()},   # work in nT
            float(np.median(gains)))


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 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)")
    args = ap.parse_args()

    t, data, gain = load(args.csv)
    keep = t >= args.start
    if args.end is not None:
        keep &= t <= args.end
    t, data = t[keep], {a: v[keep] for a, v in data.items()}
    if len(t) < 64:
        sys.exit("selected window has too few samples")

    duration = t[-1] - t[0]
    fs = (len(t) - 1) / duration
    intervals = np.diff(t)
    lsb_nt = 1000.0 / gain

    print(f"{args.csv}: {len(t):,} samples over {duration:.2f} s")
    print(f"  mean rate {fs:.1f} Hz   interval {intervals.mean()*1000:.3f} ms "
          f"+/- {intervals.std()*1000:.3f} ms   max {intervals.max()*1000:.1f} ms")
    print(f"  gain {gain:.2f} LSB/uT -> 1 LSB = {lsb_nt:.1f} nT")
    print(f"  Nyquist {fs/2:.1f} Hz\n")
    print("axis      sd       p2p    sd/LSB   white-noise ASD   median ASD")

    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 AXES:
        v = data[key]
        sd = v.std()
        # A flat (white) spectrum of this sd would sit at this level.
        implied = sd / np.sqrt(fs / 2)

        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.
        resid = v - np.polyval(np.polyfit(t, v, 3), t)
        axs[1, 0].hist(resid, bins=120, histtype="step", linewidth=1.4,
                       color=color, label=label, density=True)

        print(f"{label:4s} {sd:8.1f} {v.max()-v.min():9.1f} nT "
              f"{sd/lsb_nt:7.2f}   {implied:9.2f} nT/rtHz "
              f"{np.median(asd):9.2f} nT/rtHz")

    a = axs[0, 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.set_xlabel("frequency (Hz)"); a.set_ylabel("nT/√Hz")

    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.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 = axs[1, 1]
    a.hist(intervals * 1000, bins=120, color=REFERENCE)
    a.set_yscale("log")
    a.axvline(1000 / fs, color=TEXT_PRIMARY, linestyle="--", linewidth=1.2)
    a.annotate(f"mean {1000/fs:.2f} ms", xy=(1000 / fs, 1),
               xytext=(6, 0), textcoords="offset points",
               color=TEXT_PRIMARY, fontsize=9)
    a.set_title("Sample interval", loc="left",
                color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8)
    a.set_xlabel("ms"); a.set_ylabel("count")

    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)
    for ax in (axs[0, 0], axs[0, 1], axs[1, 0]):
        ax.legend(frameon=False, fontsize=9, labelcolor=TEXT_SECONDARY)

    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, {duration:.1f} s at {fs:.0f} Hz, "
             f"1 LSB = {lsb_nt:.1f} nT. Spectra assume uniform sampling at the "
             f"mean rate (see interval panel).",
             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()