diff --git a/NOTES.md b/NOTES.md index 8eee567..cd7aaeb 100644 --- a/NOTES.md +++ b/NOTES.md @@ -27,6 +27,9 @@ are no longer right. | "Pair cycle count with a TMRC step: 200/0x94, 400/0x95, 800/0x96" | the cycle-count rate model | Hold TMRC fast and use cycle count as the only rate knob. The pairing advice exists only because the rate was being set by TMRC | | "Abort the run on a missed measurement" | placeholder rows | The run continues and writes a `MISSED` row, keeping `sample_index` contiguous | | `~84,333 / cc / 3` for the sensor ceiling | fitting overhead separately | `1 / (3 × (cc/90,000 + 68.7 µs))`. The naive constant is not constant — it drifts from 84,429 at cc=100 to 89,191 at cc=1200 | +| "3.3 V vs 3.0 V — indistinguishable" | interleaved A/B at two cycle counts, 2026-08-23 | Indistinguishable in *noise*, yes. But the **scale differs by 7.1%** and the **sample rate by 4.94%**. The earlier entry could not see it because the sensor had moved and the two effects were confounded | +| "60 Hz is 25× the noise floor and dominates every sd figure" | the same session, at 285–299 Hz | True at ~35 Hz where it folds to 9 Hz. At the default rate mains is resolved in band at only **1.7–1.8×** the floor. At cc=400 the dominant line is not mains at all but an **fs/4 artefact** | +| "68.7 µs per-axis overhead" | solving the model per supply | Right in order, wrong in value for this unit: **40.6 µs at 3.006 V, 38.1 µs at 3.353 V**. The 90,000 counts/s divisor holds up; the overhead does not, and it moves with the supply | --- @@ -62,6 +65,10 @@ as Earth's field. ## Supply A/B: 3.3 V adapter vs 3.0 V LDO — indistinguishable +> **Superseded in the noise conclusion only, and wrong on scale.** The +> interleaved re-run below finds the scale differs by 7.1% and the sample rate by +> 4.94%. This entry could not separate either from the sensor having moved. + Two 60 s captures on the bench, cycle count 200, TMRC 0x96, 400 kHz. | | 3V3 adapter | 3.0 V LDO | @@ -122,6 +129,12 @@ magnitude difference to gain unless both movement checks pass. At ~35 Hz the largest spectral feature by far was a line near 9 Hz reaching 100–125 nT/√Hz against a 4–5 nT/√Hz broadband floor — **25× the noise floor**. +> **The 25× is specific to that rate and that site.** Measured at 285–299 Hz on +> 2026-08-23, with mains resolved in band rather than folded, the 60 Hz line +> reaches only **1.7–1.8×** the floor, and at cycle count 400 the dominant line +> is an fs/4 artefact rather than mains. The mechanism below is unchanged; the +> "dominant term" framing is not general. + It is 60 Hz aliasing, proven by making it move: | capture | fs | observed peak | predicted 60 Hz alias | @@ -435,14 +448,90 @@ the CH347's 3V3 rail ripple actually measures against the 50 mV limit. --- +## Interleaved supply A/B, 2026-08-23 — the scale factor, measured + +Four captures, LDO / 3V3 / LDO / 3V3 at cycle counts 100, 100, 400, 400, ~6 +minutes each, rails measured at the board: **3.006 V** and **3.353 V**. Full +write-up, figures and raw output in +[Noise_Floor_Testing/NOISE_FLOOR.md](Noise_Floor_Testing/NOISE_FLOOR.md); this is +the summary and the things that changed the tooling. + +**The oscillator moves with the supply.** +4.94% from 3.006 V to 3.353 V, +measured independently at two cycle counts and agreeing to **0.04%** — V^+0.44. +This is the cleanest number in the session because, unlike |B|, it cannot be +contaminated by the sensor moving. + +**|B| moves the other way, by −7.1%** (−3,245 nT), consistent across both cycle +counts to 0.74%. Exactly ratiometric with the rail would be V^−1; measured is +V^−0.65, and correcting by the rail overshoots by +3.5 to +4.3% against a +run-to-run reproducibility of 2.3–3.1%. Correcting by the **measured sample +rate** instead leaves −2.6 to −1.9%, inside that floor — and needs nothing +external, since the logger already calibrates the period every run. Applying both +corrections is clearly wrong (+8.7%): the rate change *is* the supply change seen +through the oscillator, so stacking them counts it twice. + +**The sensor still moved**, 2.8–23° of rotation between pairs. |B| survives a +rotation, so this does not invalidate the above, but translation through a +gradient would, and only a clamped re-run can bound that. + +**Cycle count 400 carries an artefact that cycle count 100 does not.** Tones at +exactly fs/4 and fs/2 of each capture's *own* rate — the two cc=400 runs differ +by 4.9% in rate and each peaks at its own fs/4, so it is locked to the sampling +rather than to any frequency. On one axis it reaches 2.43 nT at 25.7σ. Mechanism +unresolved; it is not a host-side period, since 4 samples is 52.1 ms at one rate +and 54.7 ms at the other. + +It also **grows with the supply**: comparing the rotation-invariant vector +magnitude of the tone, 3V3 is **1.217 ± 0.086×** the LDO, agreeing between the +fs/4 and fs/2 harmonics. That excludes a fixed digital pattern (which would be +1.000, +2.5σ) and external pickup (0.935, +3.3σ), and sits between V and V², +which two rail voltages cannot separate. + +**Filtering and decimating a fast capture matches sampling slowly**, which the +README had been asserting without evidence. Decimating a cc=100 capture by 4 +changes its own broadband floor by **+0.3 to +0.5%** — a comparison with no +second run in it, so no gain or position difference to confound it. It is also +strictly better: the decimated path carries none of the fs/4 artefact, and a +proper anti-alias filter removes 60 Hz before it can fold, which sampling slowly +cannot. + +**The dither question is settled at cycle count 100.** 0.65–0.81 LSB per sample, +and averaging 1,024 samples reaches 0.035–0.050 LSB — within 1.5–2.1× of ideal. +A comb in the residual histogram is not evidence against this: the residual of an +integer-valued signal is a comb at one LSB by construction, and its sharpness +tracks only how far that axis drifted. + +**What this changed in the tooling.** Four estimator bugs surfaced, all of which +would have biased a supply comparison: + +- `compare.py` scored each capture over its own band, so a 299 Hz run was judged + over 3–135 Hz and a 73 Hz one over 3–33. One common band now, chosen to sit + inside the slowest path compared *including decimated ones*. +- Welch segments were a fixed sample count, so captures at different rates got + different numbers of segments; the median of a Welch estimate is biased low by + an amount that depends on that count. A fixed 20 s segment duration now. +- Broadband noise was reported only in absolute nT/√Hz, which credits a low gain + as a low noise floor. Reported fractionally too, and the figures are fractional. +- `characterize.py` autoscaled every panel, so a quiet capture and a noisy one + drew identical-looking plots. Fixed axes now. + +--- + ## Open questions -- Re-run the LDO comparison with the sensor **clamped**, to separate gain from - movement. This is the one experiment the A/B entry above failed to perform. +- Re-run the LDO comparison with the sensor **clamped**, and in an order that is + not degenerate. The 2026-08-23 session used L,3,L,3 against cycle counts + 100,100,400,400, which makes `t = supply + 2·cc` exactly — supply, cycle count + and a linear drift cannot all three be fitted. L,3,3,L breaks it. +- Identify the **fs/4 artefact** at cycle count 400. A cc=200 run at both + supplies would say whether it scales with cycle count or appears at a threshold. +- Fix `rm3100.AXIS_OVERHEAD_S`: 68.7 µs against 38–41 µs measured on this unit, + which puts the cc100:cc400 rate ratio 1.9% out. +- Compare |B| against IGRF/WMM for the capture site (39.661262° N, + 82.604173° W, 255 m). Not a calibration standard — indoor hard and soft iron + are unknown and nothing here is corrected — but it is the only absolute + reference available, and it would say which rail reads closer to right. - Deal with 60 Hz coupling at source — standoff distance, then a repeat characterisation — rather than only sampling fast enough to see it. -- Check the dither assumption at cycle count 100 against a genuinely quiet site. - `characterize.py` prints `sd/LSB` per axis; the concern is that a sensor - quieter than spec has less dither than the 0.79 LSB the spec figure implies. - Reliability over extended runs on bare metal, where the VM ceilings above do not apply. diff --git a/README.md b/README.md index ff64bc5..05f231d 100644 --- a/README.md +++ b/README.md @@ -89,12 +89,29 @@ requirement. §4.2.3 adds a placement rule: *"Keep capacitors, especially tantal capacitors, far away from the sensor coils"* — a magnetic concern, not an electrical one. -**This rig runs at 3.3 V, but every number in Table 3-1 is quoted at 3.0 V** -(footnote 1: *"performance will differ"*). 3.3 V is comfortably in range, but the -gain fit behind `tesla_per_count` is strictly a 3.0 V figure, so the scale factor -is an extrapolation off the calibration point. A 3.0 V LDO was built and compared -— see [NOTES.md](NOTES.md), which also covers what the manual does and does not -say about ripple. +**Every number in Table 3-1 is quoted at 3.0 V** (footnote 1: *"performance will +differ"*). 3.353 V is comfortably in range, but the gain fit behind +`tesla_per_count` is strictly a 3.0 V figure, so running there is an +extrapolation off the calibration point. **This has now been measured**, on +interleaved captures at 3.006 V and 3.353 V: + +| | effect of +0.347 V | as a power of V | +|---|---|---| +| **Sample rate** | **+4.94%** (repeats to 0.04%) | V^+0.44 | +| **Total field \|B\|** | **−7.1%**, about −3,245 nT | V^−0.65 | +| Broadband noise | none resolvable, ~4% fractionally | — | + +So **supply quality does not set the noise floor, but supply *voltage* sets the +scale**, at roughly −20 %/V. A ±1% regulator at 3.0 V holds scale to ~0.2%; a +USB-derived rail free to move ±5% holds it to ~1%, on top of the ~7% fixed offset +from sitting off the calibration point. + +The chip's oscillator carries the same dependence, which is useful: the +calibrated period every capture already records tracks the rail without a +voltmeter, and correcting |B| by the measured rate removes most of the offset. +Full analysis, including what is *not* settled, in +[Noise_Floor_Testing/NOISE_FLOOR.md](Noise_Floor_Testing/NOISE_FLOOR.md). +[NOTES.md](NOTES.md) covers what the manual does and does not say about ripple. ## Setup @@ -194,23 +211,45 @@ host clock. Both appear in the capture header. ./.venv/bin/python plot.py capture.csv # -> capture.png: X, Y, Z, |B| ./.venv/bin/python plot.py capture.csv --smooth 0 # raw trace only ./.venv/bin/python characterize.py capture.csv # -> capture_noise.png + summary -./.venv/bin/python characterize.py capture.csv --start 60 --end 360 +./.venv/bin/python characterize.py capture.csv --trim 30 # drop settling/handling ./.venv/bin/python compare.py --group note *.csv # A/B two conditions +./.venv/bin/python compare.py --trim 30 --supply LDO=3.006 --supply 3V3=3.353 \ + LDO/cc100=a.csv 3V3/cc100=b.csv LDO/cc400=c.csv 3V3/cc400=d.csv -o cmp.png ./.venv/bin/python sweep.py --rates 10,50,150 # measure what each rate delivers ``` - **`plot.py`** — four stacked panels with independent y-scales, since the three axes sit at very different DC offsets. - **`characterize.py`** — amplitude spectral density against the 1.2 nT/√Hz of - Table 3-1, Allan deviation, residual distribution, and host read latency. The - `sd/LSB` column is what settles whether the quantiser is properly dithered. -- **`compare.py`** — reports noise **fractionally** (ppm of |B|), because a gain - change carries the noise with it and an absolute comparison reads a pure scale - change as a noise difference. It also refuses to attribute a magnitude - difference to gain unless the sensor demonstrably held still. + Table 3-1, Allan deviation, residual distribution, and host read latency, all + on **fixed axes** so two runs can be laid side by side. Reports `white sd` + (`sd(diff)/√2`, which rejects drift) beside the plain `sd`, and a **dither + check** — whether averaging still gets below one quantiser step. Flags any line + sitting at a simple fraction of the sample rate. Also writes a + `_spectrogram.png`: |B| in time and frequency at 0.025 Hz resolution, with + fs/4 and mains marked — where mains has folded, the line is drawn at the alias + and labelled as one. Drawn in **ppm of |B| per √Hz** on a fixed colour scale, + with native bins averaged to a fixed 0.1 Hz step so every capture gets the + same degrees of freedom per cell; both are needed before one scale across + several captures means anything. `--spectrogram`, `--overlap`, `--max-freq` + and `--colormap` tune it; the default `turbo` makes narrow lines legible, + `viridis` is the lightness-monotonic option. +- **`compare.py`** — everything cross-capture. Reports noise **fractionally** + (ppm of |B|), because a gain change carries the noise with it and an absolute + comparison reads a pure scale change as a noise difference. Also: one common + band for every capture; the timing model solved per condition; whether a line + is locked to the sampling or to a frequency; whether filtering and decimating a + fast capture matches a natively slow one; and |B| under each candidate + correction. Label captures `CONDITION/variant=path` — the part before the `/` + is the thing under test, and `--supply CONDITION=VOLTS` attaches a rail to it. - **`sweep.py`** — runs `logger.py` once per target rate and tabulates what each configuration actually achieved. +`--trim SECONDS` drops that much from *both* ends of a capture. A run usually +opens while the sensor is still settling and closes with a hand on the rig, and +neither end is a noise measurement. Nothing is dropped unless asked, and a +capture shorter than 4× the trim is analysed whole with a note rather than gutted. + A capture is analysable whether the run ended on `--duration`, on **Ctrl-C**, or early — it is valid up to wherever it stopped, so the analysis tools can be pointed at it either way. @@ -353,6 +392,12 @@ how *often* one is started. | cycle count | `3 × (cc/90,000 + 68.7 µs)` per measurement | continuous | | TMRC | how often a measurement starts | factor-of-two steps | +> The 90,000 counts/s is the specified figure and holds up; **the 68.7 µs +> overhead does not.** Solving the model against two cycle counts on this unit +> gives 40.6 µs at 3.0 V and 38.1 µs at 3.35 V, so the predicted cc100:cc400 rate +> ratio is 1.9% out. The overhead also moves with the supply, because the same +> oscillator times it — see below. + Leaving TMRC faster than the cycle-count ceiling makes the sensor free-run at ~100% duty and renders TMRC irrelevant. Setting it slower makes the sensor idle, which costs sensitivity for nothing: noise after filtering scales as `1/√duty`, @@ -367,11 +412,17 @@ cycle count alone. So the default in every branch is the fastest TMRC, and | 50 | 534 Hz | 267 Hz | 50.37 | 0.58 | 89.0% | spectrum, thin dither | | **100** | **283 Hz** | **141 Hz** | **26.17** | **0.79** | **94.2%** | **default** | | 200 | 145 Hz | 73 Hz | 13.35 | 1.10 | 97.0% | resolution, 60 Hz only | -| 400 | 74 Hz | 37 Hz | 6.74 | 1.54 | 98.5% | resolution, mains aliases | +| 400 | 74 Hz | 37 Hz | 6.74 | 1.54 | 98.5% | mains aliases, fs/4 artefact | *dither* is the sensor's own noise in LSB. Below roughly 0.2 LSB the quantiser stops being dithered and averaging no longer recovers sub-LSB resolution. +**Measured, not assumed:** at cc=100 the dither came out at 0.65–0.81 LSB and +averaging 1,024 samples reached **0.035–0.050 LSB** — a factor of 20 below the +step, within 1.5–2.1× of the ideal `1/√n`. The quantiser is not stalling. +`characterize.py` prints this for any capture; see +[Noise_Floor_Testing/NOISE_FLOOR.md](Noise_Floor_Testing/NOISE_FLOOR.md) §2. + **The default is cycle count 100 at 750 kHz**, which runs the sensor at its own ~283 Hz ceiling. It sits deliberately between the two things pulling in opposite directions: @@ -382,10 +433,20 @@ directions: and its second harmonic are in band and can be notched. - **Against cc=50** it gives up half the spectrum and buys 36% more dither margin. cc=50 is right when something above 141 Hz needs identifying; it is not - the right default, because its dither margin is the one assumption here not yet - checked against a quiet-site capture. + the right default, because its dither margin is the thinnest here and the only + one still unmeasured. -Fall back to cycle count 200 if `characterize.py` reports `sd/LSB` below ~0.35. +Fall back to cycle count 200 if `characterize.py`'s dither check shows the LSB +column flattening rather than continuing to fall. + +**Prefer decimating a fast capture over sampling slowly.** Measured on a cc=100 +capture decimated by 4 against a natively-recorded cc=400 one: decimation changes +a capture's own broadband floor by **+0.3 to +0.5%**, so the two are equivalent +for noise — and decimation is strictly better on everything else. The natively +slow capture folds 60 Hz irrecoverably to 13–17 Hz and carries sample-locked +lines at fs/4 and fs/2 that the decimated path does not have at all. `compare.py` +prints this comparison whenever two captures differ by an integer cycle-count +factor. **Bus speed is independent of the rate** — it appears in neither mechanism. What it sets is latency: how long a read takes, hence how tightly DRDY can be @@ -469,6 +530,7 @@ than merely "data arrived". | `compare.py` | A/B captures — fractional noise, plus the checks that separate a gain change from a moved sensor | | `diagnose-comms.py` | walks USB → bus → identity → registers → BIST → live read, stopping at the first failure | | `tests/` | pytest suite; needs no hardware | +| `Noise_Floor_Testing/` | captures, figures and [NOISE_FLOOR.md](Noise_Floor_Testing/NOISE_FLOOR.md) — the supply and rate-scaling analysis. Untracked: `.gitignore` excludes it | The adapter/sensor split is deliberate: `ch347.py` is a general I²C master usable with any device, and `rm3100.py` needs only a bus object exposing @@ -564,8 +626,18 @@ inserting anything. - [x] Calibrated period, configuration warnings, cycle-count rate model - [x] 3.0 V LDO built and compared — noise indistinguishable once normalised - [x] Test suite covering everything that does not need the adapter +- [x] Interleaved supply A/B at two cycle counts — **scale factor quantified at + −7.1% for +0.347 V**, rate at +4.94% +- [x] Dither margin at cc=100 checked against real data — averaging reaches + 0.035 LSB, so the quantiser is not stalling +- [x] Filter-and-decimate shown equivalent to sampling slowly (+0.5% on the + broadband floor), and better on aliasing and artefacts - [ ] Re-test the LDO with the sensor **clamped**, to separate gain from movement -- [ ] Deal with 60 Hz coupling at source — the dominant term + — every pair in the A/B still shows 2.8–23° of rotation +- [ ] Identify the **fs/4 artefact** that appears at cc=400 and not at cc=100; + a cc=200 run at both supplies would say whether it scales with cycle count +- [ ] Fix `rm3100.AXIS_OVERHEAD_S`, which is ~1.7× too large for this unit +- [ ] Deal with 60 Hz coupling at source - [ ] Reliability testing over extended runs on bare metal Captures written before the current header format are unreadable and must be diff --git a/characterize.py b/characterize.py index 4211b4f..ef69003 100644 --- a/characterize.py +++ b/characterize.py @@ -44,6 +44,7 @@ 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 @@ -84,6 +85,68 @@ 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 @@ -330,6 +393,62 @@ def sample_locked_lines(v, max_period=8, threshold=LINE_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? @@ -391,6 +510,118 @@ def decade_ticks(ax, which="both", minimum_decades=1.0): 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). @@ -428,6 +659,25 @@ def main(): 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 @@ -555,7 +805,10 @@ def main(): # 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), + # 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", @@ -640,6 +893,20 @@ def main(): 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() diff --git a/compare.py b/compare.py index 9db0cac..269608c 100644 --- a/compare.py +++ b/compare.py @@ -708,22 +708,6 @@ def binned_asd(freqs, asd, bins=140): return np.array(centres), np.array(values) -def clip_to_band(ax, series, band, headroom=3.0, floor=0.35): - """Set y-limits from what is inside the band, not from filter stopbands. - - An anti-alias filter falls away by six decades past its corner, and on a log - axis that decides the scale for everything else. Only the band being - compared should. - """ - inside = [a[(f >= band[0]) & (f <= band[1])] for f, a in series] - inside = [a for a in inside if len(a)] - if not inside: - return - lo = min(float(a.min()) for a in inside) - hi = max(float(a.max()) for a in inside) - ax.set_ylim(floor * lo, headroom * hi) - - def make_figure(records, band, supplies, path): conditions = conditions_in_order(records) variants = sorted({variant_of(l) for l, _, _ in records}) @@ -742,7 +726,6 @@ def make_figure(records, band, supplies, path): a = axs[0, 0] a.axvspan(band[0], band[1], color=ch.GRID, alpha=0.55, zorder=0, linewidth=0) - drawn = [] for label, _, s in records: freqs, asd = binned_asd(*ch.welch_asd(dict(s["cap"].axes())["total"], s["fs"], nperseg=ch.segment_length(s["fs"]))) @@ -750,11 +733,10 @@ def make_figure(records, band, supplies, path): # nT/rtHz a condition reading 7% smaller plots 7% lower while being no # quieter, and the chart would say the opposite of the truth. asd = asd / s["field"] * 1e6 - drawn.append((freqs, asd)) a.loglog(freqs, asd, color=color[condition_of(label)], linestyle=dash[variant_of(label)], linewidth=1.4, alpha=0.9, label=label) - clip_to_band(a, drawn, band) + a.set_xlim(*ch.ASD_HZ); a.set_ylim(*ch.FRACTIONAL_ASD_PPM) ch.decade_ticks(a) title(a, "Spectral density of |B|, as a fraction of |B|") a.set_xlabel("frequency (Hz)"); a.set_ylabel("ppm of |B| per √Hz") @@ -771,6 +753,7 @@ def make_figure(records, band, supplies, path): title(a, "Fast capture, filtered and decimated") a.set_xlabel("frequency (Hz)"); a.set_ylabel("ppm of |B| per √Hz") if plotted: + a.set_xlim(*ch.ASD_HZ); a.set_ylim(*ch.FRACTIONAL_ASD_PPM) ch.decade_ticks(a) a.legend(frameon=False, fontsize=9, labelcolor=ch.TEXT_SECONDARY, loc="upper right") @@ -786,6 +769,7 @@ def make_figure(records, band, supplies, path): a.loglog(taus, devs / s["field"] * 1e6, color=color[condition_of(label)], linestyle=dash[variant_of(label)], linewidth=1.6, label=label) + a.set_xlim(*ch.ALLAN_TAU_S); a.set_ylim(*ch.FRACTIONAL_ALLAN_PPM) ch.decade_ticks(a) title(a, "Allan deviation of |B|, as a fraction of |B|") a.set_xlabel("averaging time τ (s)"); a.set_ylabel("σ (ppm of |B|)") @@ -814,7 +798,6 @@ def _plot_decimation(ax, records, band, color, dash): ax.axvspan(band[0], band[1], color=ch.GRID, alpha=0.55, zorder=0, linewidth=0) plotted = False - drawn = [] for condition, (_, fast), (_, slow), k in decimation_pairs( [(label, s) for label, _, s in records]): try: @@ -826,20 +809,17 @@ def _plot_decimation(ax, records, band, color, dash): nperseg=ch.segment_length(fast["fs"] / k))) # Fractional again: the two cycle counts do not share a gain either. asd = asd / fast["field"] * 1e6 - drawn.append((freqs, asd)) ax.loglog(freqs, asd, color=hue, linestyle="--", linewidth=1.5, label=f"{condition} cc{fast['cycle_count']} ÷{k}") freqs, asd = binned_asd(*ch.welch_asd( dict(slow["cap"].axes())["total"], slow["fs"], nperseg=ch.segment_length(slow["fs"]))) asd = asd / slow["field"] * 1e6 - drawn.append((freqs, asd)) ax.loglog(freqs, asd, color=hue, linestyle="-", linewidth=1.5, alpha=0.7, label=f"{condition} cc{slow['cycle_count']} native") plotted = True if plotted: - clip_to_band(ax, drawn, band) ax.annotate("dashed: fast capture decimated. solid: natively slow.", xy=(0.02, 0.05), xycoords="axes fraction", color=ch.TEXT_SECONDARY, fontsize=9)