#!/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.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 -- 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) # 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. FRACTIONAL_ASD_PPM = (1e0, 1e3) FRACTIONAL_ALLAN_PPM = (1e0, 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 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 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] 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"{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 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}") 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) # Anchored to the left margin of the fixed scale, which no capture reaches: # the lowest bin a 20 s segment can produce is 0.05 Hz. a.annotate(f"Table 3-1 spec {SPEC_ASD_NT} nT/√Hz", xy=(ASD_HZ[0] * 1.2, 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}") # The spectrogram is a second file rather than a fifth panel: it needs the # width to resolve anything, and it answers a different question -- whether # a feature was there throughout or only for part of the run. spectrogram_path = out.rsplit(".", 1)[0] + "_spectrogram.png" if spectrogram_path.endswith("_noise_spectrogram.png"): spectrogram_path = spectrogram_path.replace("_noise_spectrogram", "_spectrogram") try: print(f"-> {make_spectrogram(cap, spectrogram_path, args.spectrogram, args.overlap, args.colormap, args.max_freq)}") except ValueError as exc: print(f"no spectrogram: {exc}") if __name__ == "__main__": main()