From 0c053e06f1cc5ccabc5c545b3e9f1d5b281bd4f4 Mon Sep 17 00:00:00 2001 From: Jeremy Karst Date: Wed, 2 Sep 2026 13:20:10 -0400 Subject: [PATCH] Several corrections and improvements to the code based on report feedback. Cleaned up and consolidated in-code documentation. --- .gitignore | 4 +- NOTES.md => Noise_Floor_Testing/NOTES.md | 307 ++++++-- README.md | 588 ++++------------ calibrate.py | 309 +++----- capture.py | 79 +-- ch347.py | 44 +- characterize.py | 859 ++++++++++------------- compare.py | 340 ++++----- diagnose-comms.py | 12 +- logger.py | 381 ++++------ plot.py | 50 +- rm3100.py | 100 +-- sweep.py | 42 +- tests/test_analysis.py | 97 +++ tests/test_calibrate.py | 241 +++---- tests/test_capture.py | 6 +- 16 files changed, 1502 insertions(+), 1957 deletions(-) rename NOTES.md => Noise_Floor_Testing/NOTES.md (62%) diff --git a/.gitignore b/.gitignore index 158db13..5f3ec0d 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,6 @@ __pycache__ *.png *.csv -/Noise_Floor_Testing \ No newline at end of file +# Captures, figures and reports stay out of the repo; the lab notebook does not. +/Noise_Floor_Testing/* +!/Noise_Floor_Testing/NOTES.md \ No newline at end of file diff --git a/NOTES.md b/Noise_Floor_Testing/NOTES.md similarity index 62% rename from NOTES.md rename to Noise_Floor_Testing/NOTES.md index d8f344e..d364d7b 100644 --- a/NOTES.md +++ b/Noise_Floor_Testing/NOTES.md @@ -1,8 +1,8 @@ # Measurement notes The lab notebook for this rig: what was measured, when, on what, and what it -changed. [README.md](README.md) documents how the code works *now*; this file is -the record of how that was arrived at, and it is append-mostly — an entry is +changed. [README.md](../README.md) documents how the code works *now*; this file is +the record of how that was arrived at, and it is append-mostly, an entry is corrected rather than deleted when a later measurement supersedes it. Every capture cited here predates the current header format and has been @@ -18,7 +18,7 @@ labelled as such. ## Superseded conclusions -Read this first — these appear as firm recommendations in earlier entries and +Read this first: these appear as firm recommendations in earlier entries and are no longer right. | Earlier conclusion | Superseded by | Now | @@ -26,15 +26,15 @@ are no longer right. | "400 kHz measured best; use 400" | three-run repetition at each speed | 750 kHz is fastest and is the default. The 4% edge for 400 kHz sat inside the run-to-run spread of a single short run | | "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**. And the divisor does not hold either — solving both freely gives **88,546 counts/s** at 3.006 V and **92,889** at 3.353 V. Holding it at 90,000 and solving for the overhead alone gives 58.9 µs at cc=100 and 113.6 µs at cc=400, fitting neither. Both terms are per-unit and supply-dependent; `calibrate.py` carries them | +| `~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 artifact** | +| "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**. And the divisor does not hold either: solving both freely gives **88,546 counts/s** at 3.006 V and **92,889** at 3.353 V. Holding it at 90,000 and solving for the overhead alone gives 58.9 µs at cc=100 and 113.6 µs at cc=400, fitting neither. Both terms are per-unit and supply-dependent; `calibrate.py` carries them | | "`× (1 + rate_error)` corrects gain, cutting the spread to 3.1%" | referencing the correction to a measured oscillator | **The 3.1% was two errors cancelling.** `rate_error` is measured against the *nominal* model, whose error is cycle-count dependent (+0.84% at cc=100, −0.99% at cc=400, same supply), and that partly cancelled the cycle-count gain difference. Referenced properly the figure is **4.9%**. The matched-cycle-count contrasts are unaffected either way | -| "gain ∝ V^−0.65, i.e. rate^−1.47; mechanism unresolved" | freeing the gain offset | **The 1.47 was absorbing a gain-offset error.** Table 3-1's `0.3671·cc + 1.5` implies gain ∝ (cc + 4.086) counts; measured on this unit it is 0.90 (LDO) and −0.09 (3V3). With the offset free, an exponent of exactly **1** is fully consistent — which is what an oscillator-clocked collection interval predicts. Correcting both leaves **2.48%**, inside the 2.3–3.1% movement floor | -| "the per-axis overhead is constant in counts (0.7σ) and changes as a time (2.8σ)" | propagating the thermal drift properly | Those σ came from a half-drift error bar chosen without justification. Propagated, the overhead is 3.60 ± 0.30 counts against 3.54 ± 0.15, and neither reading is excluded. The clean test is the supply rate ratio across cycle counts: oscillator-clocked sits **0.3σ** from it, fixed-time **1.2σ**. Favoured, not proven — drift is the limit | -| "the within-run drift could be used to correct gain in post-processing" | the coupling test, then σ_y(τ) | **Closed — do not reopen.** Gain does not follow the within-run rate: the measured slopes are **+0.08 to +0.35** where gain ∝ 1/f demands **−1**. The mechanism is that gain is set by the *ratio* of the MagI2C internal clock to f_LR, both on the same die, so a thermal shift is common-mode and divides out. That is what PNI's "stable over temperature" claim in §2 actually protects, and it is why 1241 ppm of rate drift produces no measurable gain drift. The drift is real and well resolved; it simply does not reach the gain | -| "the rate measurement is limited by host read jitter, so a longer window is always better" | σ_y(τ) on the four captures | Wrong in both halves. White jitter predicts σ/rate = √12·σ_t/(√n·T), or sub-ppm at 30 s; measured, σ_y bottoms out at **8–14 ppm around τ = 2–10 s** and then *rises*. The floor is correlated host-side noise, 10–30× the naive figure, and past ~10 s a longer window makes the estimate worse, not better — it starts measuring the drift instead of the rate | +| "gain ∝ V^−0.65, i.e. rate^−1.47; mechanism unresolved" | freeing the gain offset | **The 1.47 was absorbing a gain-offset error.** Table 3-1's `0.3671·cc + 1.5` implies gain ∝ (cc + 4.086) counts; measured on this unit it is 0.90 (LDO) and −0.09 (3V3). With the offset free, an exponent of exactly **1** is fully consistent, which is what an oscillator-clocked collection interval predicts. Correcting both leaves **2.48%**, inside the 2.3–3.1% movement floor | +| "the per-axis overhead is constant in counts (0.7σ) and changes as a time (2.8σ)" | propagating the thermal drift properly | Those σ came from a half-drift error bar chosen without justification. Propagated, the overhead is 3.60 ± 0.30 counts against 3.54 ± 0.15, and neither reading is excluded. The clean test is the supply rate ratio across cycle counts: oscillator-clocked sits **0.3σ** from it, fixed-time **1.2σ**. Favoured, not proven: drift is the limit | +| "the within-run drift could be used to correct gain in post-processing" | the coupling test, then σ_y(τ) | **Closed. Do not reopen.** Gain does not follow the within-run rate: the measured slopes are **+0.08 to +0.35** where gain ∝ 1/f demands **−1**. The mechanism is that gain is set by the *ratio* of the MagI2C internal clock to f_LR, both on the same die, so a thermal shift is common-mode and divides out. That is what PNI's "stable over temperature" claim in §2 actually protects, and it is why 1241 ppm of rate drift produces no measurable gain drift. The drift is real and well resolved; it simply does not reach the gain | +| "the rate measurement is limited by host read jitter, so a longer window is always better" | σ_y(τ) on the four captures | Wrong in both halves. White jitter predicts σ/rate = √12·σ_t/(√n·T), or sub-ppm at 30 s; measured, σ_y bottoms out at **8–14 ppm around τ = 2–10 s** and then *rises*. The floor is correlated host-side noise, 10–30× the naive figure, and past ~10 s a longer window makes the estimate worse, not better: it starts measuring the drift instead of the rate | --- @@ -56,7 +56,7 @@ are no longer right. - **Magnitude.** 54.7 µT sits inside the 25–65 µT band for Earth's field, which is the real confirmation that the gain and the int24 decode are right. - **Noise** is well above the 15 nT Table 3-1 quotes for cycle count 200, and is - not a driver problem — most likely proximity to a PC and the USB adapter. + not a driver problem, most likely proximity to a PC and the USB adapter. - **Drift.** All three axes fell monotonically over the 10 s (Z most, ~1.3 µT), so |B| went from 55.5 to 54.3 µT. Consistent with thermal settling after AVDD is first powered; coil DC resistance moves 0.4 %/°C (Table 3-3). Allow a @@ -68,7 +68,7 @@ as Earth's field. --- -## Supply A/B: 3.3 V adapter vs 3.0 V LDO — indistinguishable +## 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 @@ -84,8 +84,8 @@ Two 60 s captures on the bench, cycle count 200, TMRC 0x96, 400 kHz. **Compare noise fractionally, not absolutely.** The measured scale moved 5.9% between the runs, and a scale change carries the noise with it, so most of the -absolute sd difference is an artefact. Normalised, the two supplies differ by -2.5% — indistinguishable on single 60 s captures. +absolute sd difference is an artifact. Normalised, the two supplies differ by +2.5%, indistinguishable on single 60 s captures. That is the expected outcome. Table 3-5 allows 50 mVPP of supply ripple and the LP5907's ~6.5 µV RMS sits roughly 1000× inside it, as did the @@ -100,7 +100,7 @@ translation through a field gradient, so movement and gain are confounded here. There is still reason to expect 3.0 V to read *more accurately*, but it comes from the datasheet rather than these captures: Table 3-1 footnote 1 establishes -every specification — including the gain fit behind `tesla_per_count` — at +every specification, including the gain fit behind `tesla_per_count`, at 3.0 V, noting that performance "will differ" elsewhere. At 3.3 V the gain is an extrapolation off the calibration point. @@ -132,12 +132,12 @@ magnitude difference to gain unless both movement checks pass. ## 60 Hz mains folds into the band 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**. +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 +> is an fs/4 artifact 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: @@ -149,11 +149,11 @@ It is 60 Hz aliasing, proven by making it move: | long run | 34.5783 Hz | 9.1764 Hz | 9.1566 Hz | A real 9 Hz magnetic signal would sit still. This one moves with the sample -rate: `fs` changed by 0.20 Hz and the peak moved 0.43 Hz — the ×2 sensitivity of +rate: `fs` changed by 0.20 Hz and the peak moved 0.43 Hz, the ×2 sensitivity of a second-order fold, exactly as predicted. -The RM3100 has **no anti-alias filter** — cycle-count integration is a boxcar -with a poor stopband — so out-of-band interference folds in freely. +The RM3100 has **no anti-alias filter**: cycle-count integration is a boxcar +with a poor stopband, so out-of-band interference folds in freely. **Sampling faster resolves it instead.** At 150 Hz (Nyquist 69.3 Hz) mains lands in band and the strongest line sits at **59.923 Hz**, observed directly. At @@ -167,22 +167,22 @@ filtering can do to an alias. ## Overnight run: 13.2 h, 1.63 M samples -Cycle count 200, TMRC 0x96. **Acquisition integrity was perfect** — sample index +Cycle count 200, TMRC 0x96. **Acquisition integrity was perfect**: sample index contiguous across all 1,635,055 rows, system clock monotonic with no NTP steps, no malformed rows, no int24 saturation, no missed measurements (worst DRDY bracket 11.0 ms of 27.03 ms, 41% of margin). Both problems in the data were outside the logger. **1. The first ~1.5 h and the last ~1 h are contaminated by handling.** |B| -reaches 198,268 nT — five times Earth's field — in bursts around t+0.44 h. +reaches 198,268 nT, five times Earth's field, in bursts around t+0.44 h. Hourly |B| sd tells the story: | hours | \|B\| sd | usable | |---|---|---| -| 0–1 | 3,750 nT | no — handling | +| 0–1 | 3,750 nT | no, handling | | 2–10 | **34–40 nT** | yes | | 11–12 | 41–68 nT | marginal | -| 13 | 247 nT | no — activity resumed | +| 13 | 247 nT | no, activity resumed | A 5-minute slice from the middle of the quiet window gave |B| = 42,132 nT with sd 32.5 nT (771 ppm). @@ -197,19 +197,19 @@ show the oscillator speeding up monotonically: | 6 | 28.9794 ms | −102 ppm | | 12 | 28.9187 ms | −2195 ppm | -That is ~2,500 ppm of rate variation across the quiet hours — thermal, as +That is ~2,500 ppm of rate variation across the quiet hours, thermal, as expected of an RC oscillator. So **long captures cannot be given a single frequency axis**; spectra must be computed on windows short enough that the rate is constant, via `Capture.restrict()`. This is what `drift_limited` and the separate reporting of read jitter and fit residual exist for. Over 5 minutes the residual is 3.7 ms; over 8.5 hours it is -2,356 ms — while read jitter stays 0.21 ms in both. Conflating the two would +2,356 ms, while read jitter stays 0.21 ms in both. Conflating the two would have read as a 20,000× worse host. **Drift-limited beyond a second.** A 769 s window of the same run: per-axis sd rises to 214–284 nT while the broadband floor barely moves (5.1–5.7 nT/√Hz). The -extra spread is low-frequency drift, not broadband noise — the Allan deviation +extra spread is low-frequency drift, not broadband noise: the Allan deviation bottoms at ~13 nT near τ = 0.12 s and climbs steadily after. --- @@ -220,8 +220,8 @@ Two fixes roughly doubled the achievable rate. ### 1. Two USB round trips per register read → one -Reading a register as a pointer write plus a separate read — the form the manual -draws (§4.5.2, §5.8.4) — costs two round trips. The part accepts a **repeated +Reading a register as a pointer write plus a separate read, the form the manual +draws (§4.5.2, §5.8.4), costs two round trips. The part accepts a **repeated START**, so `ch347.write_read()` does it in one. Measured at 400 kHz: | | separate | combined | @@ -244,7 +244,7 @@ the default GIL switch interval *is* the jitter floor. ### 3. Splitting the writer off the sampling thread The same 60 s capture yielded 22,558 samples (376 Hz) with inline writing and -28,864 (481 Hz) with the writer split off — about **28%**. +28,864 (481 Hz) with the writer split off, about **28%**. ### Measured ceilings (VM, cycle count 50, 400 kHz) @@ -256,7 +256,7 @@ The same 60 s capture yielded 22,558 samples (376 Hz) with inline writing and | 0x95 | 75 Hz | clean, 69.2 Hz | 2.63–5.16 ms of 13.33 ms (20–39%) | | 0x96 | 37 Hz | clean, 34.7 Hz | 3.30–8.56 ms of 27.03 ms (12–32%) | -> **These ceilings are host-specific and pessimistic** — VM, virtualised USB, +> **These ceilings are host-specific and pessimistic**: VM, virtualised USB, > shared cores, load average 2.76 at the time of the one 150 Hz failure. The > *shape* carries over (the limit is host round-trip latency, not I²C bandwidth > or the sensor), but the numbers must be re-measured on the machine that will @@ -265,7 +265,7 @@ The same 60 s capture yielded 22,558 samples (376 Hz) with inline writing and ### The bus is never the constraint -Per sample the bus carries a DRDY poll plus a results read — 150 bit-times — +Per sample the bus carries a DRDY poll plus a results read, 150 bit-times, against a host cost 2–3× larger. | | 400 kHz | 100 kHz | @@ -273,10 +273,10 @@ against a host cost 2–3× larger. | bus time per sample | 0.375 ms | 1.50 ms | | as % of a 28.86 ms period (37 Hz) | **1.3%** | 5.2% | | as % of a 7.21 ms period (150 Hz) | **5.2%** | 20.8% | -| measured host cost | 0.89–0.98 ms | — | +| measured host cost | 0.89–0.98 ms | n/a | That puts the bus-limited ceiling near 1100 Hz, roughly double the sensor's own -562 Hz three-axis maximum at cycle count 50 — so absent host stalls, zero misses +562 Hz three-axis maximum at cycle count 50, so absent host stalls, zero misses is always achievable. Every miss observed has been a host stall. ### Bus speed, three runs each at cycle count 228 @@ -298,7 +298,7 @@ difference was inside the run-to-run spread and does not survive repetition. ## Timing: how the miss detector was arrived at **The read-to-read interval is not a usable threshold.** On the VM it reaches -37 ms against a 28.8 ms period — a 29% overshoot from host stalls alone — while +37 ms against a 28.8 ms period, a 29% overshoot from host stalls alone, while the DRDY bracket stays under 10 ms. The first implementation thresholded the interval and aborted a perfectly healthy capture at sample 88. @@ -310,7 +310,7 @@ out on this unit: |---|---|---| | nominal (TMRC table) | 27.03 ms | −6.4% | | **calibrated (1 s)** | **28.8849 ms** | **0.068%** | -| whole-run fit (25 s) | 28.8653 ms | — | +| whole-run fit (25 s) | 28.8653 ms | n/a | **Counting needs the interval, detection needs the bracket.** Fit residual after adopting the two-signal scheme: @@ -324,7 +324,7 @@ adopting the two-signal scheme: The load case is the one that shows the design paying off: the bracket reached 119% of threshold, so a stall genuinely occurred, but the interval showed only -one period had elapsed — nothing was lost. The older `int(bracket/period)` +one period had elapsed, nothing was lost. The older `int(bracket/period)` estimator would have inserted a spurious placeholder and slipped the index from there on. At 300 Hz the ambiguous counter fired 19 times. @@ -349,7 +349,7 @@ is sound over at least a minute. ``` 256 Hz is where the VM gives out; everything below it is exact. Noise past cycle -count ~400 is extrapolation — Table 3-1 ends there. +count ~400 is extrapolation, Table 3-1 ends there. **Cycle count sets the ceiling, not TMRC.** Measured at 400 kHz: @@ -364,7 +364,7 @@ count ~400 is extrapolation — Table 3-1 ends there. Table 3-1 gives 440 Hz *single-axis* and says to divide by 3, i.e. ~147 Hz. Measurement backs Table 3-1 at all three cycle counts. -**When the cycle count governs, TMRC is irrelevant** — cc=400 gave 73.85 Hz at +**When the cycle count governs, TMRC is irrelevant**: cc=400 gave 73.85 Hz at TMRC 0x92 and 73.86 Hz at 0x94. Varying only cc at TMRC 0x94: | cc | measured | governed by | duty | @@ -381,7 +381,7 @@ TMRC 0x92 and 73.86 Hz at 0x94. Varying only cc at TMRC 0x94: | TMRC-governed (cc=50, 0x94) | 138.5 Hz | 50 | **23%** | 3.42 nT/√Hz | | cycle-count-governed (cc=228) | 128.5 Hz | 228 | **98%** | 2.39 nT/√Hz | -Same chip, same field, same per-sample physics — the difference is idle time, +Same chip, same field, same per-sample physics, the difference is idle time, and it is why the current model holds TMRC fast and uses cycle count alone. **Equal duty means equal noise.** Two configurations at ~100% duty and very @@ -392,11 +392,11 @@ cc200 @150Hz 138.77 Hz LSB 13.35 nT ASD 2.98 nT/rtHz (66% bracket margin) cc800 @37Hz 34.72 Hz LSB 3.39 nT ASD 3.26 nT/rtHz (10% bracket margin) ``` -Near-identical spectral density, as the equal-duty argument predicts — so the +Near-identical spectral density, as the equal-duty argument predicts: so the choice between them is bandwidth versus resolution, not noise. **Low cycle counts are dominated, not fast.** cc=50 reached the same 138.8 Hz as -cc=200 on the VM — both host-limited, not sensor-limited — while idling 75% of +cc=200 on the VM: both host-limited, not sensor-limited, while idling 75% of the time and quantising 3.8× more coarsely. Its extra sensor speed is capability the host cannot collect. @@ -415,13 +415,13 @@ at t ≈ 41 s, so statistics are split around that: | \|B\| | 369 nT | 358 nT | +313 nT | - **The step was the sensor being physically moved**, confirmed at the time, not - an electrical artefact. A norm shift does not by itself imply an external + an electrical artifact. A norm shift does not by itself imply an external source changed: |B| is preserved under rotation in a uniform field, but a translation through a gradient samples a different local field. Near a PC the gradients are steep (§4.2.2: field falls off as 1/distance³), so a few centimetres is ample for the +313 nT seen here. Treat a |B| change as "not a pure rotation", nothing more. -- **Noise is strongly axis-dependent** — Y is 4.6× X (388 vs 84 nT) and stable +- **Noise is strongly axis-dependent**, Y is 4.6× X (388 vs 84 nT) and stable across the event, so it is not a consequence of the disturbance. The coils are nominally identical, so this points at orientation relative to a local noise source rather than a sensor fault. X at 84 nT is within ~3× of the 30 nT @@ -437,12 +437,12 @@ at t ≈ 41 s, so statistics are split around that: Table 3-5 gives a flat 50 mVPP ripple limit on AVDD and DVDD with **no frequency qualifier**, and no transfer function from ripple to field error -— no µT-per-mV, no PSRR, no ripple-vs-noise curve. +, no µT-per-mV, no PSRR, no ripple-vs-noise curve. Slow drift should matter less than fast ripple, though this is inference from the described architecture rather than an explicit claim. §4.1 describes the output as "the difference in the time to complete the measurement for each -bias" — a differential measurement across forward and reverse coil bias — and §2 +bias": a differential measurement across forward and reverse coil bias: and §2 claims measurements are "stable over temperature and inherently free from offset drift". Anything common to both half-measurements largely cancels, so supply variation slower than one measurement cycle should be rejected far better than @@ -453,16 +453,15 @@ the CH347's 3V3 rail ripple actually measures against the 50 mV limit. --- -## Interleaved supply A/B, 2026-08-23 — the scale factor, measured +## 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 +write-up, figures and raw output in [NOISE_FLOOR.md](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. +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. @@ -470,7 +469,7 @@ contaminated by the sensor moving. 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 +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. @@ -479,8 +478,8 @@ through the oscillator, so stacking them counts it twice. 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 +**Cycle count 400 carries an artifact 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 @@ -492,7 +491,7 @@ 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. -**Every capture drifts through its own thermal warm-up** — 186 to 1241 ppm, +**Every capture drifts through its own thermal warm-up**, 186 to 1241 ppm, monotonically downward, measured by re-fitting the rate on sub-windows (`calibrate.rate_stability`). That is the error bar under every timing figure here, and the reason the collection-interval mechanism is favoured rather than @@ -507,26 +506,26 @@ continuously across the session; this resets at each `logger.py` invocation. The **field** drift within a run is separate and probably is diurnal: |B| moves −1.4 to −1.8 nT/min in three of the four runs, in a consistent direction, which -is the right order for Sq variation. It cannot be the 2.48% residual — that is +is the right order for Sq variation. It cannot be the 2.48% residual: that is 1,166 nT against ~40 nT of diurnal across the whole 26-minute session. **The sensor moved during a run, not only between them.** LDO/cc400 rotates **0.49°** across its own run, ten times the quietest, and that is the same run -carrying 485 nT of Z drift at nearly constant |B| — which is rotation, not field. +carrying 485 nT of Z drift at nearly constant |B|, which is rotation, not field. The rig was not restrained and had to be handled to move the 3 V jumper between runs, so this is expected rather than surprising; it is why the residual is attributed to movement. **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 +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 +strictly better: the decimated path carries none of the fs/4 artifact, 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. +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. @@ -547,21 +546,197 @@ would have biased a supply comparison: --- +## Host setup: why the udev rule looks the way it does + +Moved here from the README, which now states the rule and not its history. + +The CH347's USB node defaults to `root:root 0664`, so it needs a rule. The one +`setup.sh` installs grants access two ways, because no single mechanism covers +every distro: + +``` +# /etc/udev/rules.d/60-ch347.rules +SUBSYSTEM=="usb", ATTRS{idVendor}=="1a86", ATTRS{idProduct}=="55db", TAG+="uaccess" +SUBSYSTEM=="usb", ATTRS{idVendor}=="1a86", ATTRS{idProduct}=="55db", GROUP="plugdev", MODE="0660" +``` + +`uaccess` has systemd-logind put an ACL on the node for whoever holds the local +seat, no group, no logout, and is the only mechanism that can work on atomic +Fedora. `GROUP`/`MODE` is the fallback for ssh sessions and seatless systems. + +**Two details are load-bearing**, and getting either wrong produces a rule that +installs cleanly and grants nothing: + +- **The prefix must sort below 73.** udev merges every rules directory into one + lexicographic sequence, and the only thing that acts on the tag is + `TAG=="uaccess", ENV{MAJOR}!="", RUN{builtin}+="uaccess"` in systemd's + `73-seat-late.rules`. At `99-` the tag is added after that line has already + run: set, never read, no ACL. systemd keeps its own uaccess rules in + `70-uaccess.rules` for this reason. +- **`TAG+="uaccess"` and `GROUP=` must be on separate lines.** udev discards a + whole rule line whose `GROUP=` it cannot resolve, and since systemd 258 that + includes any group that exists but is not a *system* group. Sharing a line + means an unusable group silently takes the uaccess tag down with it. + +Together those explain a failure that looks distro-specific but is not: on +Debian and Ubuntu the `plugdev` group grants access on its own, masking a +uaccess tag that never fired. Fedora removed `plugdev` years ago, so on +**Bazzite** neither half worked and the adapter stayed inaccessible. + +Two more atomic-Fedora traps, which `setup.sh` prints fixes for: + +- `usermod -aG dialout $USER` fails with *"group 'dialout' does not exist"* even + though `getent` finds it, because `nss_altfiles` lets `getent` read + `/usr/lib/group` while `usermod` writes `/etc/group` alone. Copy the line + across first, then re-run `usermod` and log back in. +- Bazzite is known not to reload `/etc/udev/rules.d` when it switches to the + final rootfs ([ublue-os/bazzite#2516](https://github.com/ublue-os/bazzite/issues/2516)), + so a correct rule can sit inert until + `sudo udevadm control --reload-rules && sudo udevadm trigger` is run once. + +`setup.sh` checks the node itself at the end rather than guessing, and says +*which* mechanism granted access, `uaccess ACL` or `group membership`, since +the two fail in different ways and only one is available on Bazzite. + +--- + +## Design rationale, moved out of the README + +The README states what the code does; this is why it does it that way. + +### How the time base and miss detection work + +**The chip samples on its own internal schedule**, unaffected by bus traffic +(§5.8.2, §5.7.2). Two consequences drive the design: read jitter is not sample +jitter, so `sample_index` is a grid coordinate; and a late read returns the +*newer* measurement, so an unnoticed miss would silently compress the time axis. + +A lost measurement is therefore **recorded, not dropped**: a placeholder row +keeps `sample_index` contiguous and makes the gap explicit. Double-counting is +structurally impossible: `HSHAKE` DRC1=1 makes DRDY clear only on a results read, +and `configure()` verifies the readback rather than assuming it. If the writer +queue ever fills the run **stops** rather than dropping a row, since a hole in +`sample_index` is worse than a short capture. + +**The period is calibrated before recording starts**, because counting how many +grid points passed unseen needs the real period and the TMRC table value is 6–9% +out on this unit. `calibrate_period()` fits one second of loss-free samples to +~0.07%, then keeps refining from clean intervals to follow thermal drift. If no +clean stretch can be found the run aborts before recording: that is the honest +signal that the requested rate is not sustainable. + +**Counting lost measurements uses two signals**, because the two questions have +different best answers. *Did we lose any?*, the **DRDY bracket**, which is exact +where the read-to-read interval is only a heuristic (measured here the interval +overshoots a 28.8 ms period by 29% on host stalls alone, while the bracket stays +under 10 ms). *How many?*: the **interval since the previous accepted sample**, +rounded to whole periods, since every accepted read sits on a grid point. + +**Clock drift is compensated in post-processing, not baked into the file.** +`capture.py` regresses the host clock on `sample_index` to recover the true +period: short-term regularity from the chip, long-term rate calibration from the +host. Over hours no single slope fits; `capture.py` reports that as +`drift_limited` and says to analyse shorter windows via `Capture.restrict()`. + +**Threading.** The sampling loop does I²C and nothing else, pushing raw counts +onto a bounded queue; a writer thread does unit conversion, CSV formatting, +flushing and the console display. The sampler releases the GIL inside each USB +transfer and must re-acquire it, so Python's default 5 ms switch interval becomes +the jitter floor, `sys.setswitchinterval(0.0005)` cut the worst bracket from +9.96 ms to 5.65 ms. The writer is reniced out of the way, which needs no +privilege; `--high-priority` additionally raises the sampler, which does. + +### Choosing a cycle count + +*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. At +cc=100 this is **measured, not assumed**: dither came out at 0.65–0.81 LSB and +averaging 1,024 samples reached 0.035–0.050 LSB, within 1.5–2.1× of the ideal +`1/√n`. `characterize.py` prints this for any capture; fall back to cycle count +200 if its dither check shows the LSB column flattening rather than falling. + +**The default is cycle count 100 at 750 kHz**, which runs the sensor at its own +~283 Hz ceiling. It costs 1.5% in post-filter noise against cc=200 and buys 1.9× +the spectrum, which is worth taking because **aliased interference cannot be +filtered out afterwards at any cycle count**, at 283 Hz both mains and its +second harmonic are in band and can be notched. cc=50 is right when something +above 141 Hz needs identifying, but its dither margin is the thinnest here. + +**Prefer decimating a fast capture over sampling slowly.** 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, since a natively +slow capture folds 60 Hz irrecoverably and carries sample-locked lines the +decimated path does not have at all. `compare.py` prints this comparison whenever +two captures differ by an integer cycle-count factor. + +### CH347 access + +The adapter presents three USB interfaces. Interfaces 0 and 1 are CDC-ACM and the +kernel binds them as `/dev/ttyACM0` (the UART). **Interface 2** is vendor-class, +carries I2C/SPI/JTAG with no kernel driver bound, and libusb can claim it +directly. `ch347.py` deliberately never calls `set_configuration()`: the device +is already configured, and re-setting it would disturb the CDC-ACM interfaces. + +There is no usable off-the-shelf option: no CH347 kernel driver exists (the +in-tree `spi-ch341` is for `1a86:5512`, a different chip), the PyPI `ch347` +package wraps a Windows DLL, and `ch347api` supports only HID mode (`55dc`), not +this vendor-bulk Mode 1. The framing follows the `aystarik/ch347-i2c-spi-gpio` +Linux driver, whose id-table entry +`USB_DEVICE_INTERFACE_NUMBER(0x1a86, 0x55db, 0x02)` matches this device exactly. + +Wire format, bulk OUT `0x06`, bulk IN `0x86`, max 63 bytes per transfer: + +| Purpose | Bytes out | Bytes back | +|---|---|---| +| Set speed | `AA 6 00` | 0 | +| Write n | `AA 74 (80\|n+1) (addr<<1) 75 00` | n+1 | +| Read n | `AA 74 81 ((addr<<1)\|1) [C0\|(n-1)] C0 75 00` | n+1 | +| Write n then read m | `AA 74 (80\|n+1) (addr<<1) 74 81 ((addr<<1)\|1) [C0\|(m-1)] C0 75 00` | n+2+m | +| Probe | `AA 74 81 (addr<<1) 75 00` | 1 | + +**Every returned byte must be `1`**: that is the per-byte ACK. On a read, the +leading bytes are address ACKs and the rest is payload. `tests/test_ch347.py` +pins all five forms byte for byte. + +### RM3100 quirks + +Each of these is a place where the part or the manual is not what you would +guess; `rm3100.py` carries the detail at the constant that encodes it. + +- **Register reads use the plain address, not `|0x80`.** §5 describes the SPI + convention of adding 0x80, but §5.8.4's I2C diagram writes `0x24` literally. +- **A register read uses a repeated START**, not the STOP-then-START the manual + draws (§4.5.2, §5.8.4). The part accepts it and it halves the USB round trips, + which is what sets the sample-rate ceiling. `rm3100.read_reg()` falls back to + the manual's form on a bus that cannot do it. +- **CMM = 0x79** for all three axes. Table 5-1 calls bit 3 reserved-zero, but the + manual's own examples (§5.7.2, §5.8.3) set it. The examples win. +- **HSHAKE is set to 0x0A** (`DRC0=0`, `DRC1=1`) during init. The 0x1B default has + `DRC0=1`, meaning *any* register write clears DRDY, including the pointer write + that reading STATUS itself requires, so polling could never observe DRDY set. +- **Gain** is `0.3671 × cycle_count + 1.5` LSB/µT, a linear fit to Table 3-1. +- **The rate model divisor is specified, not fitted:** a 180 kHz circuit + oscillation (Table 3-1) over two bias directions (§4.1) gives 90,000 counts/s. +- Measurements are three **24-bit big-endian two's-complement** values, read as + 9 bytes from 0x24 using the sensor's register auto-increment. + +--- + ## Open questions - 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 + 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 +- Identify the **fs/4 artifact** 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 + 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. +- Deal with 60 Hz coupling at source: standoff distance, then a repeat + characterisation, rather than only sampling fast enough to see it. - Reliability over extended runs on bare metal, where the VM ceilings above do not apply. diff --git a/README.md b/README.md index 735f038..e9f3712 100644 --- a/README.md +++ b/README.md @@ -9,8 +9,11 @@ Reference throughout: *RM3100 & RM2100 Sensor Suite User Manual*, Doc 1017252 R07 (`PNI Sensor - RM3100-Sensor-Suite-User-Manual-R07-1.pdf`, in this directory). Section numbers below refer to it. -[NOTES.md](NOTES.md) holds the measurement notebook — what was measured on this -rig, and why the design ended up as it is. +This file documents how the code works *now*. How that was arrived at, every +measurement plus a table of conclusions later measurements overturned, is in +[Noise_Floor_Testing/NOTES.md](Noise_Floor_Testing/NOTES.md), with the supply +and rate-scaling analysis in +[Noise_Floor_Testing/NOISE_FLOOR.md](Noise_Floor_Testing/NOISE_FLOOR.md). ## Quickstart @@ -19,12 +22,7 @@ rig, and why the design ended up as it is. ./.venv/bin/python diagnose-comms.py # walk the signal chain, stop at the first fault ./.venv/bin/python logger.py --duration 60 # capture to rm3100_.csv ./.venv/bin/python characterize.py rm3100_*.csv # noise floor, spectrum, stability -``` - -The test suite needs no hardware: - -```bash -./.venv/bin/python -m pytest +./.venv/bin/python -m pytest # the suite needs no hardware ``` ## Hardware @@ -35,43 +33,36 @@ The RM3100 is on a breakout board. All connections below are confirmed working. | Adapter | RM3100 | Notes | |---|---|---| | I2C VCC (3V3) | DVDD | digital supply | -| I2C VCC (3V3) | VDD (= AVDD) | analog supply — **required for measurements** | +| I2C VCC (3V3) | VDD (= AVDD) | analog supply, **required for measurements** | | I2C VCC (3V3) | I2CEN | **required to select I2C over SPI** | | I2C GND | DVSS | digital ground | | I2C SCL | SCK/SCL | pin 27, shared with SPI SCLK | | I2C SDA | SI/SDA | pin 1, shared with SPI MOSI | All three 3V3 connections share the adapter's single I2C VCC pin. Current draw is -negligible — ~260 µA per axis at 24 Hz (Table 3-1). +negligible, ~260 µA per axis at 24 Hz (Table 3-1). -Four pins are worth knowing about; two of them cause failures that are easy to -misdiagnose. +Four pins matter; two of them cause failures that are easy to misdiagnose. -**I2CEN (pin 22) must be tied HIGH.** `HIGH = I2C, LOW = SPI` (§4.3.1). SDA/SCL -are shared with SPI MOSI/SCLK and I2CEN selects between them. Left floating or -low, the chip stays in SPI mode and **never ACKs its I2C address** — the bus scan -finds nothing at all, even though the adapter is perfectly healthy. At -DVDD = 3.3 V it needs `VIH = 0.7 × DVDD` = 2.31 V, so tie it directly to 3V3. - -**AVDD/AVSS (pins 4, 5) must be powered to take measurements.** §4.3.1: *"AVDD -can be turned off when not making a measurement to conserve power, since all -other operations are supported with DVDD."* With only DVDD connected the part is -**half alive** in a way that reads as working: the address responds, REVID reads -`0x22`, registers write and read back — but the analog section driving the coils -is dead, so DRDY never asserts and the measurement registers stay at zero. BIST -reporting `XOK=YOK=ZOK=0` is the unambiguous tell. Both rails are needed, and -AVDD is not a substitute for DVDD: DVDD must come up **before or with** AVDD, -never after, and the two must stay within **0.1 V** while running (§4.3.1, -Table 3-5) — which sharing one 3V3 rail satisfies for free. - -**SA0 (pin 3) / SA1 (pin 28) set the address.** The top five bits are fixed at -`0b01000`, giving 0x20–0x23 (§4.5). On this breakout both straps are high, so the -sensor answers at **0x23**. The pins are shared with SPI SSN and MISO, so a -breakout may label them `SSN/SA0` and `MISO/SA1`. `logger.py` scans all four -addresses, so no configuration is needed. - -**DRDY (pin 23) is not connected.** Not required — the driver polls the STATUS -register instead, which the manual explicitly offers as an alternative (§5.4.1). +- **I2CEN (pin 22) must be tied HIGH.** `HIGH = I2C, LOW = SPI` (§4.3.1), and + SDA/SCL are shared with SPI MOSI/SCLK. Left floating or low the chip stays in + SPI mode and **never ACKs its I2C address**, so the bus scan finds nothing at + all on a perfectly healthy adapter. At DVDD = 3.3 V it needs + `VIH = 0.7 × DVDD` = 2.31 V, so tie it directly to 3V3. +- **AVDD/AVSS (pins 4, 5) must be powered to take measurements.** With only DVDD + the part is **half alive** in a way that reads as working: the address + responds, REVID reads `0x22`, registers write and read back, while the analog + section is dead, DRDY never asserts and the measurement registers stay zero. + BIST reporting `XOK=YOK=ZOK=0` is the unambiguous tell. DVDD must come up + **before or with** AVDD, never after, and the two must stay within **0.1 V** + while running (§4.3.1, Table 3-5), which one shared 3V3 rail gives for free. +- **SA0 (pin 3) / SA1 (pin 28) set the address**, top five bits fixed at + `0b01000`, giving 0x20–0x23 (§4.5). On this breakout both straps are high, so + the sensor answers at **0x23**. `logger.py` scans all four, so no configuration + is needed. The pins are shared with SPI SSN and MISO, so a breakout may label + them `SSN/SA0` and `MISO/SA1`. +- **DRDY (pin 23) is not connected.** The driver polls the STATUS register, + which the manual explicitly offers as an alternative (§5.4.1). ### Power @@ -84,37 +75,18 @@ register instead, which the manual explicitly offers as an alternative (§5.4.1) Reference decoupling (Figures 4-1, 4-2) is a **10 µF bulk capacitor in parallel with 0.1 µF ceramic**, rail to ground, feeding AVDD, both DVDD pins and I2CEN. -That is the entire filtering specification — no ferrite, no split rails, no LDO +That is the entire filtering specification: no ferrite, no split rails, no LDO requirement. §4.2.3 adds a placement rule: *"Keep capacitors, especially tantalum -capacitors, far away from the sensor coils"* — a magnetic concern, not an -electrical one. +capacitors, far away from the sensor coils"*, which is a magnetic concern. **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 | -|---|---| -| **Sample rate** | **+4.94%** (repeats to 0.04% across two cycle counts) | -| **Total field \|B\|** | **−7.1%**, about −3,245 nT | -| 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. - -Both effects are the same thing seen twice. The chip's measurement interval is -clocked by its own oscillator — the per-axis overhead is a fixed ~3.6 *counts*, -not a fixed time — so the same number of counts is collected over a window whose -duration goes as 1/f. Gain follows integration time, giving `gain ∝ (cc + n)/f` -with **no free exponent**. That makes the calibrated period every capture already -records a proxy for the rail, correctable without a voltmeter. 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. +differ"*), so the gain fit behind `tesla_per_count` is strictly a 3.0 V figure. +Measured on interleaved captures at 3.006 V and 3.353 V, +0.347 V moves the +sample rate by **+4.94%** and |B| by **−7.1%**, with no resolvable change in +broadband noise. So **supply quality does not set the noise floor, but supply +*voltage* sets the scale**, at roughly −20 %/V. Both effects are one thing seen +twice, and the calibrated period every capture already records is a proxy for the +rail; see [Calibration](#calibration). ## Setup @@ -125,7 +97,8 @@ including what is *not* settled, in Idempotent. It checks the adapter is present, installs a udev rule, creates `.venv` with pyusb, numpy, matplotlib and pytest, then verifies the device node is actually writable and prints whichever fix applies if it is not. `sudo` is -needed for the udev rule and nothing else. +needed for the udev rule and nothing else. Python here is PEP 668 +externally-managed, hence the venv. The rule is needed because the CH347's USB node defaults to `root:root 0664`. It grants access two ways, because no single mechanism covers every distro: @@ -136,47 +109,12 @@ SUBSYSTEM=="usb", ATTRS{idVendor}=="1a86", ATTRS{idProduct}=="55db", TAG+="uacce SUBSYSTEM=="usb", ATTRS{idVendor}=="1a86", ATTRS{idProduct}=="55db", GROUP="plugdev", MODE="0660" ``` -`uaccess` has systemd-logind put an ACL on the node for whoever holds the local -seat — no group, no logout — and is the only mechanism that can work on atomic -Fedora. `GROUP`/`MODE` is the fallback for ssh sessions and seatless systems, -using whichever of `plugdev` or `dialout` exists as a system group. - -**Two details in that file are load-bearing**, and getting either wrong produces -a rule that installs cleanly and grants nothing: - -- **The prefix must sort below 73.** udev merges every rules directory into one - lexicographic sequence, and the only thing that acts on the tag is - `TAG=="uaccess", ENV{MAJOR}!="", RUN{builtin}+="uaccess"` in systemd's - `73-seat-late.rules`. At `99-` the tag is added after that line has already - run: set, never read, no ACL. systemd keeps its own uaccess rules in - `70-uaccess.rules` for this reason. -- **`TAG+="uaccess"` and `GROUP=` must be on separate lines.** udev discards a - whole rule line whose `GROUP=` it cannot resolve, and since systemd 258 that - includes any group that exists but is not a *system* group. Sharing a line - means an unusable group silently takes the uaccess tag down with it. - -Together those explain a failure that looks distro-specific but is not: on -Debian and Ubuntu the `plugdev` group grants access on its own, masking a -uaccess tag that never fired. Fedora removed `plugdev` years ago, so on -**Bazzite** neither half worked and the adapter stayed inaccessible. - -Two more atomic-Fedora traps `setup.sh` prints fixes for: - -- `usermod -aG dialout $USER` fails with *"group 'dialout' does not exist"* even - though `getent` finds it, because `nss_altfiles` lets `getent` read - `/usr/lib/group` while `usermod` writes `/etc/group` alone. Copy the line - across first, then re-run `usermod` and log back in. -- Bazzite is known not to reload `/etc/udev/rules.d` when it switches to the - final rootfs ([ublue-os/bazzite#2516](https://github.com/ublue-os/bazzite/issues/2516)), - so a correct rule can sit inert until - `sudo udevadm control --reload-rules && sudo udevadm trigger` is run once. - -`setup.sh` checks the node itself at the end rather than guessing, and says -*which* mechanism granted access — `uaccess ACL` or `group membership` — since -the two fail in different ways and only one of them is available on Bazzite. - -Python here is PEP 668 externally-managed, hence the venv rather than a -system-wide `pip install`. +The `60-` prefix, the two separate lines and the choice of group are all +load-bearing; getting any of them wrong installs cleanly and grants nothing. +[NOTES.md](Noise_Floor_Testing/NOTES.md) explains why, and covers the +atomic-Fedora traps `setup.sh` prints fixes for. `setup.sh` checks the node +itself at the end and says *which* mechanism granted access: `uaccess ACL` or +`group membership`. ## Usage @@ -188,31 +126,32 @@ system-wide `pip install`. | Flag | Default | Meaning | |---|---|---| -| `--rate` | — | target Hz; derives the cycle count and TMRC. Mutually exclusive with `--cycle-count` | -| `--cycle-count` | 100 | sets both the rate and the LSB — see [Configuration](#configuration) | +| `--rate` | none | target Hz; derives the cycle count and TMRC. Mutually exclusive with `--cycle-count` | +| `--cycle-count` | 100 | sets both the rate and the LSB; see [Configuration](#configuration) | | `--tmrc` | `0x92` | rate register. The default is the fastest, letting the cycle count set the rate; give one only to sample *slower* | | `--duration` | 0 | seconds, 0 = until Ctrl-C | | `--output` | timestamped | CSV path | | `--address` | autodetect | skip the bus scan | | `--bus-speed` | 750 | I²C kHz. 100 kHz would spend 42% of each period on the bus at the default cycle count; 750 spends 6% | | `--calibrate` | 1.0 | seconds of loss-free samples used to measure the true period before recording. The run aborts if no clean stretch can be found | -| `--note` | — | free text recorded in the capture header, e.g. the supply under test | +| `--note` | none | free text recorded in the capture header, e.g. the supply under test | | `--high-priority` | off | raise the sampling thread to nice −10. Needs `CAP_SYS_NICE`, so run under `sudo`; output files are handed back to the invoking user | -| `--scan-only` | — | scan and exit | +| `--scan-only` | none | scan and exit | -**Exit status is non-zero if the capture is compromised in any way** — a lost +**Exit status is non-zero if the capture is compromised in any way**: a lost measurement, an ambiguous gap, a truncated run, or an undrained writer. Each prints an explanation on stderr, and the file is always valid as far as it goes. +A capture is analysable whether the run ended on `--duration`, on **Ctrl-C**, or +early. -Before recording, `logger.py` prints how the configuration was derived, so it can -be checked rather than trusted, then measures the true sample period against the -host clock. Both appear in the capture header. +Before recording, `logger.py` prints how the configuration was derived, so it +can be checked, then measures the true sample period against the host clock. +Both appear in the capture header. ### Analysis tools ```bash ./.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 --trim 30 # drop settling/handling ./.venv/bin/python compare.py --group note *.csv # A/B two conditions @@ -221,60 +160,40 @@ host clock. Both appear in the capture header. ./.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, 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. - Two further figures come out beside it: - `_drift.png` — the chip's oscillator measured against the host clock, as rate - against elapsed time with error bars, and beside it the two-sample deviation - σ_y(τ) that supplies those bars. The minimum of that curve is the best the - rate can be known and the averaging time worth using; the branch rising ∝ τ to - its right is real drift, not measurement noise. These parts run 186–1241 ppm - of warm-up against a floor near 10 ppm at τ ≈ 2–10 s. - `_timebase.png` — for each detected line, its coherent amplitude computed on - the chip's sample grid and again on the host clock. A line at a fixed - frequency sharpens on the host clock; one locked to the sampling is wrecked by - it. That settles **from a single capture** what `compare.py` otherwise needs - two different rates to decide. Diagnostic only: no other figure uses the host - time base, because neither choice is right for both kinds of line. -- **`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. `--calibration-out` fits the timing model - across the sweep and writes it as a calibration. -- **`calibrate.py`** — converts a raw capture to a nanotesla CSV, applying a - calibration when one is named and nothing when one is not. See below. +Each tool's `--help` carries its full reasoning; in one line each: -`--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. +- **`plot.py`**: four stacked panels with independent y-scales, since the three + axes sit at very different DC offsets. +- **`characterize.py`**: spectral density against Table 3-1, Allan deviation, + residual distribution and host read latency on **fixed axes**, plus a dither + check and any line at a simple fraction of the sample rate. Writes three more + figures beside it: `_spectrogram.png` (whether a feature lasted or was a + burst), `_drift.png` (the chip's oscillator against the host clock, with + σ_y(τ) supplying the error bars) and `_timebase.png` (which clock each line is + coherent on, which settles from one capture what `compare.py` needs two rates + to decide). +- **`compare.py`**: everything cross-capture: fractional noise in ppm of |B|, + one common band for every capture, the timing model per condition, whether a + line is locked to the sampling or to a frequency, whether decimating a fast + capture matches a natively slow one, and |B| under each candidate correction. + Label captures `CONDITION/variant=path`; `--supply CONDITION=VOLTS` attaches a + rail. +- **`sweep.py`**: runs `logger.py` once per target rate and tabulates what each + configuration achieved. `--calibration-out` fits the timing model across the + sweep and writes it as a calibration. +- **`calibrate.py`**: converts a raw capture to a nanotesla CSV, applying a + calibration when one is named and nothing when one is not. + +`--trim SECONDS` drops that much from *both* ends of a capture, since a run +usually opens while the sensor is settling and closes with a hand on the rig. Nothing is +dropped unless asked, and a capture shorter than 4× the trim is analysed whole +with a note. ### Calibration **Nothing is corrected by default.** `calibrate.py` without `--calibration` changes units and nothing else, and the analysis tools read raw captures as they -always did. That is deliberate: the corrections below are real but not yet -settled enough to apply silently. +always did. The corrections are real but not settled enough to apply silently. ```bash ./.venv/bin/python sweep.py --rates 32,64,128,256 \ @@ -285,49 +204,32 @@ settled enough to apply silently. A calibration is **per unit and per supply**, because the chip times everything from one oscillator and that oscillator is neither the specified frequency nor -fixed. Solving `period = 3 × (cc/C + overhead)` freely against one unit gave -C = 88,546 with a 40.6 µs overhead at 3.006 V, against the nominal 90,000 and -68.7 µs — and C moved to 92,889 at 3.353 V. +fixed. It is worth having even if you never correct a field, because it makes the +timing model right: with a measured pair the predicted rate matches the measured +one at both cycle counts, where the nominal model is 1.9% out between them. Two +corrections then become available, neither on by default: -That is worth having even if you never correct a field, because it makes the -timing model right: with the measured pair, the predicted rate matches the -measured one to −0.00% at both cycle counts, where the nominal model is 1.9% out -between them. Once it is right, **one capture recovers the oscillator** — -`C = cc / (period/3 − overhead)` — to about 0.09% at cycle count 100. +| | what it is | +|---|---| +| oscillator gain | `(C / reference) ^ gain_exponent`, where the exponent is **1**, which an oscillator-clocked measurement interval predicts | +| cycle-count gain | `gain ∝ (cc + gain_offset_counts)`; Table 3-1 implies n = 4.086, measurement puts it near zero | -Two corrections then become available, and neither is on by default: +Every correction arrives with a propagated standard deviation, from the timing +fit and from `rate_stability()`, which measures how far the oscillator drifted +during a run and reads **186 to 1241 ppm** here. Let the rig settle first. -| | what it is | status | -|---|---|---| -| oscillator gain | `(C / reference) ^ gain_exponent` | −7.1% for +0.347 V. The exponent is **1**, which an oscillator-clocked measurement interval predicts and the measured rates support at 0.3σ against 1.2σ for the alternative | -| cycle-count gain | `gain_offset_counts` | `gain ∝ (cc + n)`. Table 3-1 implies n = 4.086; measured on one unit it is 0.90 and −0.09. Together the two corrections take a 10.3% disagreement to 2.5%, against a 2.3–3.1% movement floor | - -Every correction arrives with a propagated standard deviation — 0.3–0.4% on -these captures, well under the residual, which is what says the residual is real -rather than calibration slop. It comes from the timing fit and from -`rate_stability()`, which measures how far the oscillator drifted during a run: -**186 to 1241 ppm** here, monotonically downward. Let the rig settle first. - -The output is a different format and `capture.py` will refuse to read it, on -purpose — raw captures stay the one source of truth. Its header carries the -source's SHA-256 so a derived file that has drifted from its source is -detectable rather than merely plausible. Full analysis in -[Noise_Floor_Testing/NOISE_FLOOR.md](Noise_Floor_Testing/NOISE_FLOOR.md). +The output is a different format and `capture.py` refuses to read it, which +keeps raw captures the one source of truth. Its header carries the source's +SHA-256, so a derived file can be tied back to what produced it. > `Capture.rate_error` is measured against the **nominal** model, whose error is -> cycle-count dependent — +0.84% at cc=100 and −0.99% at cc=400 on one unit, same -> supply. It is a configuration diagnostic, not a gain measurement, and -> multiplying a field by `1 + rate_error` is not a correction. `calibrate.py` -> does that against a measured reference instead. - -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. +> cycle-count dependent, which makes it a configuration diagnostic. +> `calibrate.py` corrects gain against a measured per-unit reference. ## Capture format -A capture stores **only irreducible facts** — the chip's sample count, the host -clock, and the raw counts — behind a header of configuration. Chip time, elapsed +A capture stores **only irreducible facts**, the chip's sample count, the host +clock and the raw counts, behind a header of configuration. Chip time, elapsed time, tesla and magnitude are all reconstructed on load by `capture.py`, so there is exactly one source of truth for each. @@ -351,14 +253,13 @@ sample_index,system_time_unix,x_raw,y_raw,z_raw,warning ``` `tesla_per_count` is written with `repr()` so it round-trips through float64 -exactly, and is expressed per-count rather than the datasheet's LSB/µT so -conversion is a multiply: `tesla = count * tesla_per_count`. All times are unix -epoch seconds; there is no ISO-8601 anywhere in the data. +exactly, and is expressed per-count, so conversion from the datasheet's LSB/µT +is a multiply. All times are unix epoch seconds. The `warning` column carries space-separated flags, empty when the row is fine, and generalises to future flags. Because it already marks a row as having no -data, placeholders carry **zeros rather than blanks**, so `x/y/z` parse as -integers on every row without special-casing. +data, placeholders carry **zeros**, which keeps `x/y/z` parsing as integers on +every row. | flag | meaning | |---|---| @@ -366,91 +267,16 @@ integers on every row without special-casing. | `AMBIGUOUS` | the gap ending at this row was of uncertain length, so the index may have slipped from here on | **The two are independent.** A gap measuring 1.35 periods rounds to one, so no -placeholder is written — yet it sits far enough from an integer to distrust, and +placeholder is written, yet it sits far enough from an integer to distrust, and that case carries `AMBIGUOUS` on the *real* sample ending the gap, which keeps -its data. A placeholder inside an uncertain gap carries both. - -The flag is a confidence measure, not a claim that 1.35 and 1.65 are equally -likely — 1.35 probably is one period and 1.65 probably two. Deciding which needs -neighbouring timestamps and assumptions and can still be wrong, which is exactly -why these rows are flagged rather than silently resolved. +its data. A placeholder inside an uncertain gap carries both. The flag is a +confidence measure: 1.35 probably is one period and 1.65 probably two, but +deciding which needs neighbouring timestamps and assumptions, so the row carries +the flag and the reader decides. `capture.py` interpolates `MISSED` rows on load so the uniform grid the spectra depend on survives, and reports how many were substituted so it is never silent. -## How it works - -**The chip samples on its own internal schedule**, unaffected by bus traffic -(§5.8.2 *"this can run in the background"*; §5.7.2 *"This will not affect the -measurement process"*). Two consequences drive the whole design: - -- **Read jitter is not sample jitter.** The measurement grid stays uniform even - when host reads are late, so `sample_index` is a grid coordinate. -- **A late read returns the *newer* measurement**, not a delayed one. An - unnoticed miss therefore skips a grid point and silently compresses the time - axis. - -So a lost measurement is **recorded, not dropped**: a placeholder row keeps -`sample_index` contiguous and makes the gap explicit. The run continues and -reports the total at the end. Double-counting is structurally impossible — -`HSHAKE` DRC1=1 makes DRDY clear only on a results read, and `configure()` -verifies the readback rather than assuming it. - -**The period is calibrated before recording starts.** Counting how many grid -points passed unseen needs the real period, and the TMRC table value is 6–9% out -on this unit — enough to insert the wrong number of placeholders and slip the -index. `calibrate_period()` takes a least-squares slope over one second of -loss-free samples, which pins it to ~0.07%, then keeps refining it from clean -intervals so it follows the oscillator's thermal drift. If no clean stretch can -be found the run aborts before recording anything: that is the honest signal that -the requested rate is not sustainable. - -**Counting lost measurements uses two different signals**, because the two -questions have different best answers: - -- *Did we lose any?* — the **DRDY bracket**, the span between the last poll - showing DRDY clear and the poll showing it set. This is exact rather than - heuristic: if DRDY reads clear at t꜀ then every earlier measurement has already - been read, measurements complete one period apart, so a bracket narrower than a - period can contain at most one completion — and DRDY going high proves it - contained at least one. -- *How many?* — the **interval since the previous accepted sample**, rounded to - whole periods. Every accepted read sits on a grid point, so that interval is - nearly an exact multiple. The bracket cannot say, since it only reaches back to - the last poll that saw DRDY clear and so discards where the grid is. - -The read-to-read interval alone cannot make the first claim: measured here it -reaches 37 ms against a 28.8 ms period — a 29% overshoot from host stalls alone — -while the bracket stays under 10 ms. Thresholding the interval flags healthy -captures as lossy. - -**Clock drift is compensated in post-processing, not baked into the file.** The -RC oscillator is *regular* but only accurate to ±7% (§5.2.1). `capture.py` -regresses the host clock on `sample_index` to recover the true period: - -``` -elapsed_nominal = N * dt_nominal # uniform, but wrongly scaled -elapsed = N * dt_true # uniform and correctly scaled -``` - -Short-term regularity comes from the chip, long-term rate calibration from the -host. `characterize.py` uses `dt_true` for its frequency axis, since a 6% error -would displace every spectral feature by 6%. Over hours the oscillator drifts -enough that no single slope fits; `capture.py` reports that as `drift_limited` -and says to analyse shorter windows. - -**Threading.** The sampling loop does I²C and nothing else, pushing raw counts -onto a bounded queue; a writer thread does unit conversion, CSV formatting, -flushing and the console display. The sampler releases the GIL inside each USB -transfer and must re-acquire it, so Python's default 5 ms switch interval becomes -the jitter floor — `sys.setswitchinterval(0.0005)` cut the worst bracket from -9.96 ms to 5.65 ms. The writer is also reniced out of the way, which needs no -privilege; `--high-priority` additionally raises the sampler, which does. - -If the queue ever fills, the run **stops** rather than dropping a row: a missing -row would leave a hole in `sample_index`, and a truncated capture is worth more -than a longer one no tool can load. - ## Configuration **Cycle count is the rate knob, not TMRC.** Two ceilings compete and the slower @@ -462,18 +288,15 @@ 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. +> The 90,000 counts/s is specified and holds up. **The 68.7 µs overhead does +> not**: solving the model on this unit gives 40.6 µs at 3.0 V and 38.1 µs at +> 3.35 V, and it moves with the supply because the same oscillator times it. 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`, -and a measured 23% duty cost 1.43× the noise ASD against the same rate reached by -cycle count alone. So the default in every branch is the fastest TMRC, and -`--rate` derives the cycle count from there. +~100% duty. Setting it slower makes the sensor idle, which costs sensitivity for +nothing: noise after filtering scales as `1/√duty`, and a measured 23% duty cost +1.43× the noise ASD against the same rate reached by cycle count alone. So the +default in every branch is the fastest TMRC, and `--rate` derives the cycle count. **Cycle count also sets resolution**, which is the real trade: @@ -482,43 +305,15 @@ 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% | mains aliases, fs/4 artefact | +| 400 | 74 Hz | 37 Hz | 6.74 | 1.54 | 98.5% | mains aliases, fs/4 artifact | -*dither* is the sensor's own noise in LSB. Below roughly 0.2 LSB the quantiser +*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. +`characterize.py` prints the measured figure for any capture, and cycle count +200 is the fallback where its dither check shows the LSB column flattening out. The notebook has the reasoning behind the default and behind preferring +to decimate a fast capture over sampling slowly. -**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: - -- **Against cc=200** it costs 1.5% in post-filter noise for 1.9× the spectrum - (141 Hz of Nyquist against 73 Hz). Worth taking, because **aliased interference - cannot be filtered out afterwards at any cycle count** — at 283 Hz both mains - 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 thinnest here and the only - one still unmeasured. - -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 +**Bus speed is independent of the rate**, appearing in neither mechanism. What it sets is latency: how long a read takes, hence how tightly DRDY can be timestamped and how much margin there is against a stall. Host cost is `bus time + ~0.6 ms` of fixed USB round trip, so faster is simply better. @@ -536,14 +331,14 @@ bounds are enforced by `rm3100.plan()`. ### Configuration warnings `rm3100.plan()` resolves a configuration, `logger.py` prints its derivation, then -checks it against six known traps. Each is a *silent* failure — the capture -completes, the numbers look plausible, and the defect only shows up afterwards. +checks it against six known traps. Each is a *silent* failure: the capture +completes, the numbers look plausible, and the defect shows up afterwards. So each is reported on stderr and **nothing is fixed automatically**: changing a setting that was asked for would hide the problem behind a configuration change. | Warning | Trigger | Why it matters | |---|---|---| -| cycle count below the recommended 50 | `cc < RECOMMENDED_MIN_CYCLE_COUNT` | dither thins toward the ~0.2 LSB where averaging stops recovering sub-LSB resolution | +| cycle count below the recommended 50 | `cc < RECOMMENDED_MIN_CYCLE_COUNT` | duty falls, and post-filter noise goes as `1/√duty` | | TMRC governs and the sensor idles | TMRC-governed, idle > 20% | idle time buys nothing; measured 1.43× the ASD at 23% duty | | rate differs from the one requested | \|error\| > 2% | the run silently uses the ceiling, and every derived figure moves with it | | Nyquist below 60 Hz | rate < 120 Hz | mains folds onto signal and no later filter undoes it | @@ -557,29 +352,29 @@ the line still visible. ## Diagnostics -`diagnose-comms.py` walks the chain — USB → bus → identity → registers → BIST → -live read — and stops at the first failure, so a fault points at a specific wire. +`diagnose-comms.py` walks the chain, USB → bus → identity → registers → BIST → +live read, and stops at the first failure, so a fault points at a specific wire. Otherwise, symptoms map onto causes: | Symptom | Cause | |---|---| -| `Cannot claim CH347 interface 2: Access denied` | udev rule missing or not yet applied to this node — replug the adapter, then run `./setup.sh`, which reports whether `uaccess` or group membership granted access and diagnoses whichever fell through. On Bazzite also try `sudo udevadm control --reload-rules && sudo udevadm trigger`. Do **not** reach for `sudo`: it works as a normal user, and masking a permissions problem with root only defers it | +| `Cannot claim CH347 interface 2: Access denied` | udev rule missing or not yet applied to this node; replug the adapter, then run `./setup.sh`. On Bazzite also try `sudo udevadm control --reload-rules && sudo udevadm trigger`. Do **not** reach for `sudo`: it works as a normal user, and masking a permissions problem with root only defers it | | `No CH347 adapter found` | not plugged in, or not in Mode 1 | -| `Could not find N s of loss-free samples` | the host cannot sustain this rate — lower it (`--rate`, or a higher `--tmrc`) or raise `--bus-speed` | +| `Could not find N s of loss-free samples` | the host cannot sustain this rate; lower it (`--rate`, or a higher `--tmrc`) or raise `--bus-speed` | | `N measurement(s) were lost` | host stalls during the run. Same fixes; the capture is still usable, with explicit gaps | -| `N gap(s) could not be counted confidently` | the index may have slipped — re-record before doing spectral work | -| `no capture header found` | a capture predating the header format — re-record it | +| `N gap(s) could not be counted confidently` | the index may have slipped; re-record before doing spectral work | +| `no capture header found` | a capture predating the header format; re-record it | | `sample_index is not contiguous` | the file was truncated mid-row, damaged, or hand-edited | | `HSHAKE did not take` | the I²C write is unreliable; try a lower `--bus-speed` | | Adapter opens, `--scan-only` finds nothing | **I2CEN not tied high** (most likely), or SDA/SCL swapped, or no bus pull-ups | -| Found at 0x23, REVID `0x22`, registers fine, but DRDY never sets and results are all zero | **AVDD/VDD not powered** — confirm with BIST | -| Cycle-count read-back mismatch | bus integrity — try a lower speed | +| Found at 0x23, REVID `0x22`, registers fine, but DRDY never sets and results are all zero | **AVDD/VDD not powered**; confirm with BIST | +| Cycle-count read-back mismatch | bus integrity; try a lower speed | **BIST is the definitive test for the analog side** (§5.6.1). Write `0x8F` to BIST (STE=1, max timeout and periods), write `0x70` to POLL, wait, then read BIST back: bits 4/5/6 are XOK/YOK/ZOK, and `1` means that axis's LR oscillator ran. -All zeros means the coils are not oscillating, pointing at AVDD or the REXT -timing resistor rather than anything on the I²C side. +All zeros means the coils are not oscillating, which points at AVDD or the REXT +timing resistor. A healthy total field magnitude is roughly **25–65 µT** (Earth's field). Near zero, railed, or wildly out of range means the decode or the gain is wrong rather @@ -590,84 +385,28 @@ than merely "data arrived". | File | Role | |---|---| | `setup.sh` | udev rule, venv, dependencies, node permission check | -| `ch347.py` | CH347 I²C transport — USB only, knows nothing about the sensor | -| `rm3100.py` | RM3100 driver and configuration model — knows nothing about USB | +| `ch347.py` | CH347 I²C transport; USB only, knows nothing about the sensor | +| `rm3100.py` | RM3100 driver and configuration model; knows nothing about USB | | `logger.py` | CLI: scan → identify → configure → calibrate → log, recording lost measurements as explicit placeholders | -| `capture.py` | the only capture reader — parses the header, rebuilds tesla and both time bases | +| `capture.py` | the only capture reader; parses the header, rebuilds tesla and the time base | | `plot.py` | four-panel plot of a capture: X, Y, Z and the norm | -| `characterize.py` | noise floor: spectral density, Allan deviation, residuals, read latency | +| `characterize.py` | noise floor and stability: spectral density, Allan deviation, residuals, read latency, oscillator drift | | `sweep.py` | sweep target rates, reporting measured rate, LSB, noise, duty and bus use per point | -| `compare.py` | A/B captures — fractional noise, plus the checks that separate a gain change from a moved sensor | +| `compare.py` | A/B captures: fractional noise, plus the checks that separate a gain change from a moved sensor | | `calibrate.py` | convert a raw capture to nanotesla, applying a per-unit calibration if one is given; also holds the calibration model | | `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 | +| `Noise_Floor_Testing/` | the lab notebook, captures, figures and the supply analysis. Only `NOTES.md` is tracked | 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 -`write(addr, data)` and `read(addr, count)` — optionally `write_read()`, which it +`write(addr, data)` and `read(addr, count)`, optionally `write_read()`, which it prefers when available. -## Implementation notes - -Things that cost time to work out, recorded so they don't have to be again. - -### CH347 access - -The adapter presents three USB interfaces. Interfaces 0 and 1 are CDC-ACM and the -kernel binds them as `/dev/ttyACM0` (the UART). **Interface 2** is vendor-class, -carries I2C/SPI/JTAG with no kernel driver bound, and libusb can claim it -directly with nothing to detach. `ch347.py` deliberately never calls -`set_configuration()` — the device is already configured, and re-setting it would -disturb the CDC-ACM interfaces driving the UART. - -There is no usable off-the-shelf option: no CH347 kernel driver exists (the -in-tree `spi-ch341` is for `1a86:5512`, a different chip), the PyPI `ch347` -package wraps a Windows DLL, and `ch347api` supports only HID mode (`55dc`), not -this vendor-bulk Mode 1. The framing follows the `aystarik/ch347-i2c-spi-gpio` -Linux driver, whose id-table entry -`USB_DEVICE_INTERFACE_NUMBER(0x1a86, 0x55db, 0x02)` matches this device exactly. - -Wire format — bulk OUT `0x06`, bulk IN `0x86`, max 63 bytes per transfer: - -| Purpose | Bytes out | Bytes back | -|---|---|---| -| Set speed | `AA 6 00` | 0 | -| Write n | `AA 74 (80\|n+1) (addr<<1) 75 00` | n+1 | -| Read n | `AA 74 81 ((addr<<1)\|1) [C0\|(n-1)] C0 75 00` | n+1 | -| Write n then read m | `AA 74 (80\|n+1) (addr<<1) 74 81 ((addr<<1)\|1) [C0\|(m-1)] C0 75 00` | n+2+m | -| Probe | `AA 74 81 (addr<<1) 75 00` | 1 | - -**Every returned byte must be `1`** — that is the per-byte ACK. On a read, the -leading bytes are address ACKs and the rest is payload. `tests/test_ch347.py` -pins all five forms byte for byte. - -### RM3100 quirks - -- **Register reads use the plain address, not `|0x80`.** §5 describes the SPI - convention of adding 0x80, but the I2C diagram in §5.8.4 writes `0x24` - literally. Only 7 bits are decoded, so both work; the plain form is used here. -- **A register read uses a repeated START**, not the STOP-then-START the manual - draws (§4.5.2, §5.8.4). The part accepts it — verified against REVID and the - measurement registers — and it halves the USB round trips, which is what sets - the sample-rate ceiling. `rm3100.read_reg()` falls back to the manual's form on - a bus that cannot do it. -- **CMM = 0x79** for all three axes. Table 5-1 describes bit 3 as reserved-zero, - but the manual's own examples (§5.7.2, §5.8.3) set it. The examples win. -- **HSHAKE is set to 0x0A** (`DRC0=0`, `DRC1=1`) during init. The 0x1B default has - `DRC0=1`, meaning *any* register write clears DRDY — including the pointer write - that reading STATUS itself requires, so polling could never observe DRDY set. - With `DRC0=0`, DRDY clears only on a results read, which is what makes - exactly-once sampling possible. -- **Gain** is `0.3671 × cycle_count + 1.5` LSB/µT, a linear fit to Table 3-1 - (50→20, 100→38, 200→75; reproduces all three within a count). -- **The rate model divisor is specified, not fitted.** Table 3-1 gives a 180 kHz - circuit oscillation and §4.1 measures each cycle count in *both* bias - directions, so one count costs two oscillations — hence 90,000 counts/s. The - 68.7 µs per-axis overhead has no specified value and was measured; it is why - the naive `rate × cc` "constant" drifts from 84,429 at cc=100 to 89,191 at 1200. -- Measurements are three **24-bit big-endian two's-complement** values, read as - 9 bytes from 0x24 using the sensor's register auto-increment. +The CH347 wire format, the RM3100 register quirks and the reasoning behind the +time base are in +[NOTES.md](Noise_Floor_Testing/NOTES.md#design-rationale-moved-out-of-the-readme); +`tests/test_ch347.py` pins every framing form byte for byte. ## Tests @@ -686,38 +425,11 @@ failures are invisible downstream: a stall of 1.1 periods must lose nothing, 2.1 must record one placeholder, and 1.35 must flag the sample as `AMBIGUOUS` without inserting anything. -## Status - -- [x] udev rule, venv and dependency setup -- [x] CH347 vendor protocol — verified byte-for-byte against the kernel driver -- [x] I²C bus scan, sensor identified at 0x23, REVID `0x22` -- [x] Analog section verified — BIST reports `XOK=YOK=ZOK=1` -- [x] Exactly-once sampling, chip-grid time base, explicit gaps for lost samples -- [x] Verified on hardware: exact grid, −6.11% drift recovered, 0.017% run-to-run -- [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 - — 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 -- [x] Per-unit timing and gain calibration — `calibrate.py`, written by - `sweep.py --calibration-out`, nothing corrected by default -- [ ] Measure a calibration on the bench and check the oscillator against a - third rail voltage, which would separate V from V² in the gain scaling -- [ ] 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 -re-recorded — `capture.py` says so explicitly rather than guessing. Figures -quoted in [NOTES.md](NOTES.md) came from those files and stand as historical -measurements only. +re-recorded; `capture.py` says so explicitly. Figures quoted in the notebook +came from those files and stand as historical measurements only. A proper Python API is planned; these modules are the initial-communication milestone, deliberately kept simple. diff --git a/calibrate.py b/calibrate.py index 6714e7a..dc69429 100644 --- a/calibrate.py +++ b/calibrate.py @@ -7,51 +7,30 @@ Without `--calibration` this only changes the format: counts become nanotesla using the `tesla_per_count` the capture already carries, and nothing is -corrected. That is the default deliberately -- see "What is not corrected". +corrected. That is deliberate: the corrections are real but rest on a small +number of captures taken with a sensor that moved between runs, so they are +applied only when asked for, and every factor carries the standard deviation it +earned. See Noise_Floor_Testing/NOTES.md for the measurements behind them. -**Why a calibration is per-unit.** The chip times everything from one oscillator, -and `rm3100`'s model of it, `period = 3 x (cycle_count / C + overhead)`, uses a -specified C and an overhead measured on one rig. Neither fits every part. Solving -the model against one unit's two cycle counts gave C = 88,546 with a 40.6 us -overhead at 3.006 V, against the model's 90,000 and 68.7 us -- and C moved to -92,889 when the supply moved to 3.353 V, because the oscillator does. - -That matters beyond timing. Gain moves with the oscillator, so a capture taken at -a different supply reads a different field for the same magnet. A calibration -records the timing pair and a reference oscillator frequency, and the correction -is how far this capture's oscillator sits from that reference. +**A calibration is per unit and per supply.** The chip times everything from one +oscillator, and `rm3100`'s model of it, `period = 3 x (cycle_count / C + +overhead)`, uses a specified C and a nominal overhead that fit no particular +part. A calibration records a measured pair plus a reference oscillator +frequency; the correction is how far this capture's oscillator sits from that +reference. Gain moves with the oscillator, so a capture taken at another supply +reads a different field for the same magnet. Once the timing pair is right, **one capture is enough** to recover the -oscillator: rearranging the model gives `C = cycle_count / (period/3 - overhead)`, -and a 1 us error in the overhead moves C by only 0.09% at cycle count 100. With -the model's 28 us error it moves by 2.5%, and by a *different* amount at each -cycle count, which is what makes the uncalibrated figures unusable for this. +oscillator: `C = (cycle_count + n) / (period/3)`, where n is the overhead in +counts. `gain_exponent` is 1 because the measurement interval is clocked +entirely by that oscillator. The overhead is a fixed *count* of cycles, so the +same number of counts is collected over a window whose length goes as 1/f, and +gain follows integration time. `oscillator_hz` rests on the same premise, which +keeps the two consistent. -**Why the exponent is 1.** The measurement interval is clocked entirely by that -oscillator: the per-axis overhead measures 3.60 counts at 3.006 V and 3.54 at -3.353 V while its *duration* moves 6%, so the same number of counts is always -collected over a window whose length goes as 1/f. Gain then follows integration -time, giving `gain ~ (cc + n) / f` -- exponent exactly 1. Measured rates agree: -the supply rate ratio is the same at both cycle counts to 0.04%, which a -fixed-time overhead would not produce. It is favoured rather than proven, at 0.3 -sigma against 1.2 for a fixed-time overhead; thermal drift is what limits it. - -**Why there is a gain offset.** Table 3-1's `0.3671*cc + 1.5` LSB/uT implies the -gain goes as `cc + 4.086` in count units. Solving the same-supply cycle-count -contrast on one unit gives 0.90 and -0.09 -- near zero, so the overhead counts -cost time without integrating field. Leaving it at the datasheet value accounts -for only half of what the correction can remove. - -**What is not corrected, and why nothing is by default.** All of it was measured -on four captures, two per supply, with a sensor that rotated between runs by 2.8 -to 23 degrees, and every capture drifts 186 to 1241 ppm through its own thermal -warm-up. So it is applied only when a calibration is named, every factor arrives -with the standard deviation it earned, and the output records what was done. - -The output is a different format from a raw capture and `capture.py` will refuse -to read it, which is intended: raw captures stay the one source of truth, and a -derived file that drifts from its source should be detectable rather than -plausible. The header carries the source's SHA-256 for that reason. +The output is a different format from a raw capture and `capture.py` refuses to +read it, which keeps raw captures the one source of truth. The header carries +the source's SHA-256, so a derived file can be tied back to what produced it. """ import argparse @@ -70,7 +49,7 @@ import characterize as ch import rm3100 # Bumped when either file format changes incompatibly, so a reader can refuse -# rather than misinterpret. +# a file it would misread. CALIBRATION_VERSION = 2 CALIBRATED_VERSION = 1 @@ -78,11 +57,6 @@ CALIBRATED_VERSION = 1 # no reason to hold one in memory twice. HASH_BLOCK = 1 << 20 -# Sub-windows used to measure how much the oscillator moved during a capture. -# Eight over a five-minute run is ~40 s each, long enough to fit a rate against -# and short enough to see a thermal transient. -STABILITY_WINDOWS = 8 - class CalibrationError(Exception): pass @@ -91,10 +65,10 @@ class CalibrationError(Exception): class Estimate(namedtuple("Estimate", "value sd")): """A number and its standard deviation, carried together. - Every corrected field here rests on a chain of measurements, and a residual - only means something against the uncertainty of the thing that was supposed - to remove it. Returning the two together is what makes that comparison - possible at the point of use rather than in a comment. + Every corrected field rests on a chain of measurements, and a residual only + means something against the uncertainty of the thing that was supposed to + remove it. Carrying the two together puts that comparison at the point of + use. """ __slots__ = () @@ -103,9 +77,6 @@ class Estimate(namedtuple("Estimate", "value sd")): def relative(self): return self.sd / abs(self.value) if self.value else float("inf") - def __format__(self, spec): - return f"{self.value:{spec or '.6g'}} +/- {self.sd:{spec or '.2g'}}" - Calibration = namedtuple("Calibration", [ "counts_per_second", # C in the timing model, this unit's own @@ -113,7 +84,7 @@ Calibration = namedtuple("Calibration", [ "reference_oscillator_hz", # C at which the gain is taken as correct "reference_cycle_count", # cycle count at which it is taken as correct "gain_exponent", # gain ~ (cc + offset) / f ** exponent - "gain_offset_counts", # the "+ n" -- see gain_factor() + "gain_offset_counts", # the "+ n"; see gain_factor() "counts_per_second_sd", "axis_overhead_s_sd", "gain_exponent_sd", @@ -133,14 +104,13 @@ def calibration(counts_per_second, axis_overhead_s, A calibration whose reference is its own oscillator frequency corrects nothing on the unit and supply it was taken at, and corrects other captures - relative to it. That is the useful default: the reference is a choice of - where "right" is, and the only defensible choice without an absolute field - standard is the condition the gain was characterised at. + relative to it. That is the useful default: without an absolute field + standard, the only defensible choice of where "right" is, is the condition + the gain was characterised at. - The uncertainties all default to zero, so a calibration that makes no claim - about its own accuracy still works -- and every field it corrects then - reports a standard deviation of zero, which is a visible claim rather than a - silent one. + The uncertainties default to zero, so a calibration making no claim about + its own accuracy still works. Every field it corrects then reports a + standard deviation of zero, which states that claim on the face of the file. """ return Calibration( counts_per_second=float(counts_per_second), @@ -163,8 +133,7 @@ def _finite(value, field): """A field that must parse as a real number. Booleans are not numbers here. JSON's `true` reaches Python as a bool, which floats to 1.0 without - complaint -- a silent way for a malformed calibration to become a plausible - one. + complaint, so booleans are rejected here. """ if isinstance(value, bool) or not isinstance(value, (int, float, str)): raise CalibrationError(f"{field}: {value!r} is not a number") @@ -207,7 +176,7 @@ def load_calibration(path): f"{path}: rm3100_calibration is {version!r}, this reads " f"{CALIBRATION_VERSION}. Version 1 carried a per-cycle-count gain " "table; version 2 replaces it with gain_offset_counts, which is " - "not derivable from it -- re-measure rather than convert.") + "not derivable from it; re-measure it.") missing = {"counts_per_second", "axis_overhead_s"} - raw.keys() if missing: raise CalibrationError(f"{path}: missing {', '.join(sorted(missing))}") @@ -267,89 +236,6 @@ def save_calibration(path, cal): return path -def window_rates(cap, seconds): - """Sample rate fitted independently on each consecutive window of a capture. - - Returns (centres, rates) -- seconds from the start of the capture, and Hz -- - with one point per window that held enough samples to fit. - - `Capture.dt_true` is a single straight line through a run whose rate is - moving. This is the same fit repeated on pieces short enough that the - movement within any one of them is small, which is what turns the drift from - an error on the whole-capture number into something measurable in its own - right. - """ - seconds = float(seconds) - start = float(cap.elapsed[0]) - if seconds <= 0 or cap.duration < 2 * seconds: - return np.array([]), np.array([]) - centres, rates = [], [] - for w in range(int(cap.duration // seconds)): - low = start + w * seconds - try: - rate = cap.restrict(low, low + seconds).true_rate_hz - except capture.CaptureError: - continue - rates.append(rate) - centres.append((w + 0.5) * seconds) - return np.array(centres), np.array(rates) - - -def rate_stability(cap, windows=STABILITY_WINDOWS): - """Fractional spread of the sample rate across sub-windows of a capture. - - The chip's RC oscillator warms up, and every capture here drifts downward - through its own run -- 186 ppm on the quietest, 1241 ppm on the worst. A - single fitted `dt_true` is the average over that transient, so this is the - honest error bar on it, and the input to every uncertainty downstream. - - It is not the same thing as `Capture.drift_limited`, which asks whether one - straight line describes the capture at all. This asks how far the rate moved - while it did. Nor is it `rate_allan`, which asks how much of that movement - is real: a spread this wide is only meaningful against the noise of the - measurement that produced it. - """ - if len(cap.sample_index) < 64 * windows: - return 0.0 - _, rates = window_rates(cap, cap.duration / windows) - if len(rates) < 2: - return 0.0 - return float((rates.max() - rates.min()) / rates.mean()) - - -def rate_allan(cap, taus): - """Two-sample (Allan) deviation of the measured sample rate, fractional. - - For each averaging time tau the rate is fitted on consecutive tau-second - windows and the deviation taken between neighbours: - - sigma_y(tau) = sqrt( mean( (y[k+1] - y[k])^2 ) / 2 ) - - Differences between neighbours rather than a spread about the mean, for the - same reason `characterize.white_sd` differences rather than taking a - standard deviation: on a series that drifts, the spread measures the drift. - Here both effects are wanted and the shape of the curve separates them. - Falling with tau is measurement noise averaging down; rising with tau is - drift the measurement has begun to resolve, and a rise proportional to tau - specifically is a deterministic frequency ramp rather than a random walk. - - The minimum is therefore the best the rate can be known, and the tau where - it falls is how long it is worth measuring for. - - Returns (taus, sigma_y) over the taus that yielded three or more windows. - Two windows give a single difference, which is not an estimate of anything. - """ - kept, devs = [], [] - for tau in taus: - _, rates = window_rates(cap, tau) - if len(rates) < 3: - continue - y = rates / rates.mean() - kept.append(float(tau)) - devs.append(float(np.sqrt(np.mean(np.diff(y) ** 2) / 2))) - return np.array(kept), np.array(devs) - - TimingFit = namedtuple("TimingFit", "counts_per_second axis_overhead_s " "counts_per_second_sd axis_overhead_s_sd " "residual_ppm points") @@ -359,17 +245,16 @@ def fit_timing(cycle_counts, periods, period_sd=None): """Solve `period = AXES * (cycle_count / C + overhead)` for both terms. `periods` are whole sample periods in seconds, not per-axis. `period_sd` is - the fractional uncertainty on each -- `rate_stability()` is where it comes - from. + the fractional uncertainty on each, which `rate_stability()` supplies. - Both terms are needed. Holding C at the specified 90,000 and solving for the - overhead alone gave 58.9 us at cycle count 100 and 113.6 us at 400 on the - same unit -- no single value fits, because the misfit is in both. + Both terms are solved for, since the misfit sits in both: holding C at the + specified 90,000 gives 58.9 us of overhead at cycle count 100 and 113.6 us + at 400 on the same unit. - **Two points determine both terms exactly and leave no residual**, so a bad - point cannot be detected: `residual_ppm` is 0 and means nothing. That case - requires `period_sd`, so the uncertainty is at least stated rather than - fabricated from a perfect fit. Three or more points carry their own. + **Two points determine both terms exactly and leave no residual**, so + `residual_ppm` is 0 and carries no information there. That case requires + `period_sd`, which states the uncertainty outright. Three or more points + carry their own. """ counts = np.asarray(cycle_counts, dtype=float) per_axis = np.asarray(periods, dtype=float) / rm3100.AXES @@ -418,33 +303,58 @@ def fit_timing(cycle_counts, periods, period_sd=None): float(np.abs(residual / per_axis).max() * 1e6), len(counts)) +def overhead_counts(cal): + """The calibration's per-axis overhead expressed in oscillator cycles. + + The overhead is stored as a duration, which is what `fit_timing` solves for, + and that duration belongs to the C it was measured with: the overhead is a + fixed *count* of cycles whose duration goes as 1/f. See `oscillator_hz`. + """ + return cal.axis_overhead_s * cal.counts_per_second + + def oscillator_hz(cap, cal, period_sd=None): """This capture's oscillator frequency, with its uncertainty. - Inverts `period = AXES * (cycle_count / C + overhead)`. The overhead comes - from the calibration because it cannot be recovered from a single cycle - count -- one equation, two unknowns -- and it is the term the nominal model - gets most wrong. + Inverts the timing model. The overhead comes from the calibration, a single + cycle count leaving one equation in two unknowns, and it is the term the + nominal model gets most wrong. + + **The overhead is carried as a count of cycles**, the same premise the gain + exponent rests on, so the model keeps it inside the count and closes in one + step: + + period / AXES = (cycle_count + overhead_counts) / C + + Written with the calibration's *duration* instead, the two forms disagree by + 0.23% on a figure whose propagated standard deviation is 0.30%. `period_sd` is the fractional uncertainty on the measured period; it defaults to measuring it from the capture with `rate_stability()`. """ - per_axis = cap.dt_true / rm3100.AXES - cal.axis_overhead_s - if per_axis <= 0: + overhead = overhead_counts(cal) + # The overhead is a handful of counts against a cycle count of 30 or more. + # One that rivals the cycle count would set the recovered frequency, so the + # calibration is refused. + if overhead >= cap.cycle_count: raise CalibrationError( - f"{cap.path}: a {cap.dt_true * 1e3:.3f} ms period leaves no time " - f"for {cap.cycle_count} counts once the calibration's " - f"{cal.axis_overhead_s * 1e6:.1f} us overhead is removed. The " - "calibration does not describe this capture.") + f"{cap.path}: the calibration's {cal.axis_overhead_s * 1e6:.1f} us " + f"overhead is {overhead:,.1f} counts against a " + f"{cap.cycle_count}-count measurement, so the overhead, not the " + "measurement, would set the answer. The calibration does not " + "describe this capture.") + counts = cap.cycle_count + overhead + per_axis = cap.dt_true / rm3100.AXES if period_sd is None: - period_sd = rate_stability(cap) - value = cap.cycle_count / per_axis + period_sd = ch.rate_stability(cap) + value = counts / per_axis # Two independent contributions: how well the period was measured, and how - # well the overhead subtracted from it is known. The second dominates at low - # cycle counts, where the overhead is a larger share of the per-axis time. - from_period = cap.dt_true * period_sd / rm3100.AXES - relative = math.hypot(from_period, cal.axis_overhead_s_sd) / per_axis - return Estimate(value, value * relative) + # well the overhead added to the count is known. The second dominates at low + # cycle counts, where the overhead is a larger share of the total count. + # Both go into fractional terms first, since they enter from opposite ends. + from_period = period_sd + from_overhead = cal.axis_overhead_s_sd * cal.counts_per_second / counts + return Estimate(value, value * math.hypot(from_period, from_overhead)) def gain_factor(cap, cal, period_sd=None): @@ -454,18 +364,20 @@ def gain_factor(cap, cal, period_sd=None): oscillator how far this capture's clock sits from the reference, raised to the fitted exponent. The measurement interval is clocked - entirely by that oscillator -- the per-axis overhead measures - 3.57 counts at both supplies, constant, while its *duration* - moves 6% -- so the same number of counts is always collected - over a window whose length goes as 1/f. That predicts an - exponent of exactly 1. + entirely by that oscillator, which predicts an exponent of + exactly 1. gain shape Table 3-1's `0.3671*cc + 1.5` LSB/uT implies the gain goes as - `cc + 4.086` in count units. Measured on one unit it is 0 to - 0.9, so the overhead counts cost time without integrating + `cc + 4.086` in count units, where measurement puts it near + zero, the overhead counts costing time while integrating no field. `gain_offset_counts` is that "+ n". reference the shape term is only known up to an overall scale, since `gain = A(cc + n)` leaves A free. `reference_cycle_count` pins it: the correction is exactly 1 there. + + **A matched-cycle-count contrast cannot measure the offset.** The shape term + depends only on the cycle count, so two conditions compared at the *same* + cycle count carry the same shape and it cancels from their ratio exactly. + Only a contrast spanning cycle counts says anything about `n`. """ oscillator = oscillator_hz(cap, cal, period_sd) ratio = oscillator.value / cal.reference_oscillator_hz @@ -562,10 +474,9 @@ def write_calibrated(cap, path, source, cal=None, factor=None, class Calibrated: """A capture already converted to nanotesla, read back from disk. - Deliberately not a `capture.Capture`: the raw counts and the host clock are - gone, so the time base cannot be re-fitted and `restrict()` would be a lie. - What survives is what a plot or a spectrum needs, plus the provenance of how - it got that way. + A separate type from `capture.Capture`, since the raw counts and the host + clock are gone and the time base is fixed at what was written. What survives + is what a plot or a spectrum needs, plus the provenance of how it got there. """ def __init__(self, path, meta, sample_index, elapsed, x, y, z, total): @@ -602,19 +513,8 @@ class Calibrated: def load_calibrated(path): """Read a file written by write_calibrated().""" - meta, rows = {}, [] with open(path, newline="") as handle: - while True: - position = handle.tell() - line = handle.readline() - if not line: - break - if not line.startswith("#"): - handle.seek(position) - break - key, _, value = line[1:].partition(":") - if value: - meta[key.strip()] = value.strip() + meta = capture.parse_header(handle) rows = list(csv.DictReader(handle)) version = meta.get("rm3100_calibrated") @@ -634,7 +534,9 @@ def load_calibrated(path): if not rows: raise CalibrationError(f"{path}: no samples") - column = lambda name: np.array([float(r[name]) for r in rows]) + def column(name): + return np.array([float(r[name]) for r in rows]) + return Calibrated( path, meta, np.array([int(r["sample_index"]) for r in rows], dtype=np.int64), @@ -645,9 +547,8 @@ def load_calibrated(path): def load_any(path): """Load a raw capture or a calibrated one, whichever this is. - Which format a file is, is a property of the file, not of the caller -- so - tools take a path and get back something with axes() and a time base either - way, rather than each having to know. + Which format a file is, is a property of the file, so tools take a path and + get back something with axes() and a time base either way. """ with open(path) as handle: head = handle.read(4096) @@ -672,7 +573,7 @@ def main(): args = ap.parse_args() # Check the destination before reading a capture that may be tens of - # megabytes, so a mistyped directory fails in a second rather than a minute. + # megabytes, so a mistyped directory fails in a second. parent = os.path.dirname(os.path.abspath(args.output)) or "." if not os.path.isdir(parent): sys.exit(f"{args.output}: {parent} is not a directory") @@ -694,9 +595,9 @@ def main(): print("no calibration given: converting units only, correcting nothing") else: measured = oscillator_hz(cap, cal) - stability = rate_stability(cap) + stability = ch.rate_stability(cap) print(f"calibration {args.calibration}" - + (f" -- {cal.note}" if cal.note else "")) + + (f": {cal.note}" if cal.note else "")) print(f" rate {cap.true_rate_hz:.4f} Hz, drifting " f"{stability * 1e6:.0f} ppm across the run") print(f" oscillator {measured.value:,.0f} +/- {measured.sd:,.0f} Hz " diff --git a/capture.py b/capture.py index 57194cd..a3a6ef4 100644 --- a/capture.py +++ b/capture.py @@ -1,18 +1,14 @@ """Read an RM3100 capture and reconstruct everything the file does not store. -A capture holds only irreducible facts -- sample count, host clock, raw counts -- -plus a header of configuration. This module turns that back into field values and -a time base, and is the single place either is derived. +A capture holds only irreducible facts, the sample count, host clock and raw +counts, plus a header of configuration. This module turns that back into field +values and a time base, and is the single place either is derived. -Two time bases come out of `load()`: - - elapsed_nominal N * dt_nominal, using the TMRC table rate - elapsed N * dt_true, using the rate actually observed - -The chip's RC oscillator is accurate to about +/-7% (manual section 5.2.1) and -our unit runs ~6% slow, so `elapsed_nominal` is uniform but wrongly scaled. -`dt_true` comes from regressing the host clock on the sample index: the chip -supplies short-term regularity, the host supplies long-term rate calibration. +`elapsed` is `sample_index * dt_true`, where `dt_true` comes from regressing the +host clock on the sample index: the chip supplies short-term regularity, the +host supplies long-term rate calibration. The header's nominal rate is a model +no particular part obeys, out by up to ~6% in either direction depending on +cycle count and supply, so the time base is measured instead. """ import csv @@ -22,7 +18,7 @@ import numpy as np import rm3100 # Flags logger.py writes in the `warning` column, space separated. Duplicated -# rather than imported so that reading a capture never pulls in the USB stack. +# here so that reading a capture never pulls in the USB stack. WARN_MISSED = "MISSED" # placeholder row: no data, keeps the index contiguous WARN_AMBIGUOUS = "AMBIGUOUS" # gap ending here was of uncertain length @@ -32,7 +28,7 @@ class CaptureError(Exception): class Capture: - """A loaded capture: raw counts, nanotesla axes, and two time bases.""" + """A loaded capture: raw counts, nanotesla axes and a drift-corrected clock.""" def __init__(self, path, meta, sample_index, system_time, counts, missed=None, ambiguous=None): @@ -60,35 +56,26 @@ class Capture: # Least-squares fit of host clock against grid coordinate, with both # origins removed first. polyfit rescales the Vandermonde columns but # never centres the ordinate, so fitting raw epoch seconds (~1.7e9) for - # a ~4 ms slope loses six significant digits to cancellation -- measured - # at 1.8e-6 of relative slope error against 3.8e-9 once centred. - # Subtracting a constant from each axis leaves the slope untouched and - # only moves the intercept, which nothing outside this fit uses. + # a ~4 ms slope loses six significant digits to cancellation. Shifting + # each axis leaves the slope untouched and only moves the intercept. index_0, time_0 = sample_index[0], system_time[0] slope, intercept = np.polyfit(sample_index - index_0, system_time - time_0, 1) self.dt_true = float(slope) # Error in the *rate*, so the sign matches the reported Hz: negative - # means the chip samples slower than the nominal table value. - # - # Against the *nominal* model, which no particular part obeys. Its own - # error is cycle-count dependent -- on one unit, +0.84% at cycle count - # 100 and -0.99% at 400, same supply -- so this is a diagnostic of the - # configuration, not a measurement of the oscillator, and multiplying a - # field by (1 + rate_error) is not a gain correction. calibrate.py does - # that against a measured per-unit reference. + # means the chip samples slower than the nominal table value. It is a + # diagnostic of the configuration against a nominal model. calibrate.py + # corrects gain against a measured per-unit reference. self.rate_error = self.dt_nominal / self.dt_true - 1.0 self.residuals = (system_time - time_0) - ( intercept + slope * (sample_index - index_0)) - self.elapsed_nominal = sample_index * self.dt_nominal self.elapsed = sample_index * self.dt_true - # Two different things live in the timing error, and conflating them is - # misleading. Read latency is local: the scatter of one read interval - # about the next. Fit residual is global: how far the whole capture - # departs from a single straight line, which on a long run is dominated - # by the oscillator's rate drifting with temperature, not by the host. + # Two different things live in the timing error. Read latency is local: + # the scatter of one read interval about the next. Fit residual is + # global: how far the whole capture departs from a single straight line, + # which on a long run is the oscillator drifting with temperature. self.read_jitter = float(np.diff(system_time).std()) if len(system_time) > 1 else 0.0 self.residual_sd = float(self.residuals.std()) self.residual_span = float(np.ptp(self.residuals)) @@ -112,9 +99,8 @@ class Capture: def restrict(self, start=None, end=None): """Return a new Capture covering a window of drift-corrected seconds. - Re-fits on the subset, so a window's rate is its own rather than - inherited -- which is what makes it usable for spotting drift across a - long run by comparing windows. + Re-fits on the subset, so each window carries its own rate, which is + what makes comparing windows a drift measurement. """ keep = np.ones(len(self.sample_index), dtype=bool) if start is not None: @@ -148,9 +134,6 @@ class Capture: f"(residual sd {self.residual_sd * 1e3:.0f} ms", f" against a {self.dt_true * 1e3:.1f} ms period)."] if self.missed.any(): - # Each lost measurement is an independent chance to insert one - # placeholder too many or too few, and the error accumulates in - # sample_index. That is the likelier cause here than thermal drift. lines += [ f" With {int(self.missed.sum()):,} lost measurement(s) the likely cause is " "miscounted placeholders:", @@ -179,7 +162,7 @@ class Capture: return "\n".join(lines) -def _parse_header(handle): +def parse_header(handle): """Consume leading '# key: value' lines, leaving the reader at the CSV.""" meta = {} while True: @@ -199,26 +182,23 @@ def _parse_header(handle): def load(path): """Load a capture written by logger.py.""" with open(path, newline="") as handle: - meta = _parse_header(handle) + meta = parse_header(handle) rows = list(csv.DictReader(handle)) if "rm3100_capture" not in meta: raise CaptureError( f"{path}: no capture header found. Files written before the header " - "format was introduced cannot be read -- re-record them.") + "format was introduced cannot be read; re-record them.") missing = {"cycle_count", "tesla_per_count", "nominal_rate_hz"} - meta.keys() if missing: raise CaptureError(f"{path}: header missing {', '.join(sorted(missing))}") if len(rows) < 64: raise CaptureError(f"{path}: only {len(rows)} samples, too few to analyse") - # logger.py writes "MISSED" in the count fields for a measurement it could - # not read, keeping sample_index contiguous so the chip-time grid stays - # valid across the gap. if "warning" not in rows[0]: raise CaptureError( f"{path}: no 'warning' column. Captures predating it cannot be " - "read -- re-record them.") + "read; re-record them.") flags = [set(r["warning"].split()) for r in rows] missed = np.array([WARN_MISSED in f for f in flags]) # Independent of missed: a gap can be uncertain yet round to zero losses, @@ -228,10 +208,9 @@ def load(path): raise CaptureError(f"{path}: every row is a lost measurement") sample_index = np.array([int(r["sample_index"]) for r in rows], dtype=np.int64) - # logger.py keeps the index contiguous even across lost measurements -- it - # writes MISSED placeholders rather than skipping grid points, and stops the - # run outright if it ever cannot. So a gap here means the file was truncated - # mid-row, damaged, or hand-edited, not merely cut short. + # logger.py writes a placeholder at every grid point it misses, and stops + # the run outright if it ever cannot, so a gap here means the file was + # truncated mid-row, damaged, or hand-edited. gaps = np.diff(sample_index) if np.any(gaps != 1): bad = int(sample_index[np.argmax(gaps != 1)]) @@ -245,8 +224,6 @@ def load(path): counts = {} good = ~missed for a in "xyz": - # Every row parses as an integer -- placeholders carry zeros, and the - # warning column is what marks them as having no data. v = np.array([int(r[f"{a}_raw"]) for r in rows], dtype=np.float64) if missed.any(): v[missed] = np.interp(np.flatnonzero(missed), np.flatnonzero(good), diff --git a/ch347.py b/ch347.py index 67a223a..1d30fc3 100644 --- a/ch347.py +++ b/ch347.py @@ -1,7 +1,7 @@ """I2C master over a CH347 USB adapter (Waveshare USB to UART/I2C/SPI/JTAG). -Knows nothing about any particular I2C device -- it only moves bytes to and -from a 7-bit address. +Knows nothing about any particular I2C device; it moves bytes to and from a +7-bit address. The CH347 in Mode 1 (UART1+I2C+SPI) presents three interfaces: two CDC-ACM interfaces that the kernel binds for the UART, and interface 2, a vendor-class @@ -37,18 +37,16 @@ CMD_IN = 0xC0 # read; low 6 bits are the byte count CMD_SET = 0x60 # set bus speed; low bits select the rate # Bus speeds, as the low nibble of CMD_SET. The fastest is the default: bus -# speed does not set the sample rate -- the sensor does -- but it does set how -# long a read takes, and so how tightly an event on the far side of the bus can -# be timestamped. Measured host cost per sample is bus time plus ~0.6 ms of -# fixed USB round-trip latency: 2.13 ms at 100 kHz, 0.92 at 400, 0.79 at 750. +# speed sets how long a read takes, and so how tightly an event on the far side +# of the bus can be timestamped. The sensor sets the sample rate. Host cost per +# sample is bus time plus ~0.6 ms of fixed USB round-trip latency. SPEED_20KHZ = 0 SPEED_100KHZ = 1 SPEED_400KHZ = 2 SPEED_750KHZ = 3 # The one place the rate constants are tied to their kHz figures. Callers that -# take a speed from a user or print one in a report go through this rather than -# keeping their own copy. +# take a speed from a user or print one in a report go through this. SPEEDS = { 20: SPEED_20KHZ, 100: SPEED_100KHZ, @@ -110,13 +108,13 @@ class CH347I2C: """Send one command stream, then read in_len status/data bytes. Returns (data, mono, wall), the clocks read the instant the reply - lands -- before any parsing or unwinding of the call stack. A caller - timing an event on the far side of the bus (DRDY going high, say) wants - that moment, not one several Python frames later: the frames add both a - systematic lag and jitter from whatever the interpreter does between. + lands, before any parsing or unwinding of the call stack. A caller + timing an event on the far side of the bus wants that moment; several + Python frames later carries a systematic lag and the interpreter's + jitter. - The stamp travels with its own reply rather than being left on the - instance, so it cannot go stale or be picked up by the wrong call. + The stamp travels with its own reply, which ties it to the transaction + that produced it. """ written = self._dev.write(EP_OUT, bytes(out_bytes), TIMEOUT_MS) if written != len(out_bytes): @@ -183,19 +181,15 @@ class CH347I2C: def write_read(self, addr, data, count): """Write data then read count bytes in one transaction. - Returns (data, mono, wall) -- see _xfer() on why the stamp is returned - rather than stored. + Returns (data, mono, wall); see _xfer() on the stamp. - Uses a repeated START rather than STOP-then-START, so the whole - exchange is a single USB round trip instead of two. Measured at - 400 kHz this halves a one-byte register read (0.65 -> 0.35 ms), which - matters because the host round trip, not the bus, sets the sample-rate - ceiling. + Uses a repeated START rather than STOP-then-START, so the exchange is a + single USB round trip, which halves a one-byte register read. The host + round trip is what sets the sample-rate ceiling. The RM3100's manual draws its register reads with a STOP between the - pointer write and the read (sections 4.5.2, 5.8.4), but the part - accepts a repeated START -- verified against REVID and the measurement - registers. + pointer write and the read (sections 4.5.2, 5.8.4), but the part accepts + a repeated START, verified against REVID and the measurement registers. """ data = bytes(data) if not data: @@ -213,7 +207,7 @@ class CH347I2C: packet += [CMD_IN, CMD_STO, CMD_END] # One status byte per clocked-out byte: the write address, the payload, - # and the read address -- then the payload itself. + # and the read address, then the payload itself. acks = len(data) + 2 reply, mono, wall = self._xfer(packet, acks + count) if any(byte != ACK for byte in reply[:acks]): diff --git a/characterize.py b/characterize.py index 2dbb164..6aafc38 100644 --- a/characterize.py +++ b/characterize.py @@ -5,54 +5,39 @@ 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. + Amplitude spectral density nT/sqrt(Hz) against Table 3-1's noise figure + expressed at this capture's own rate. Allan deviation where averaging stops helping and drift takes - over -- the honest measure of a noise floor. + 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. + Read latency host-side diagnostics. The measurement grid is + uniform whatever this shows. Three further figures are written beside it, each answering something the four panels collapse away: _spectrogram.png whether a feature lasted the run or was a burst. _drift.png how far the chip's oscillator moved against the - host clock, and -- from the two-sample deviation - of the rate -- how much of that is real. These - parts drift 186-1241 ppm against a floor near - 10 ppm, so it is a signal by a wide margin. - _timebase.png which clock each line is coherent on. A line at a - fixed frequency sharpens when the samples are - placed by the host clock; one locked to the - sampling is wrecked by it, which settles from a - single capture what compare.py needs two rates to - decide. Diagnostic only: every other figure here - uses the chip's grid, because neither time base is - right for both kinds of line. + host clock, and how much of that is real. + _timebase.png which clock each line is coherent on, settled + from a single capture where compare.py needs two + rates. Diagnostic only: every other figure uses + the chip's grid. -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. +Two noise figures are printed. `sd` is the plain standard deviation, which on a +drifting capture measures the drift. `white sd` is sd(diff)/sqrt(2), which +rejects drift, and is what Table 3-1's 208/sqrt(cycle_count) compares 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. +`--trim` drops equal time from both ends: settling at the start, handling at +the end. 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. +Frequencies use the rate measured against the host clock, never the nominal +table value, which is out by a cycle-count and supply dependent amount. 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. +a capture with many of them reads as artificially smooth at high frequency; +capture.py's summary reports how many. """ import argparse @@ -79,29 +64,16 @@ 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 field is derived from the three axes, so it takes neutral ink. TOTAL = ("total", "Total", TEXT_PRIMARY) SERIES = AXES + [TOTAL] -# Table 3-1: "Noise Density @ Max. Single-Axis Sample Rate". -SPEC_ASD_NT = 1.2 - -# Every panel is drawn on a fixed scale so two runs can be laid side by side and -# compared by eye. Autoscaling defeats that completely: a quieter capture simply -# redraws its own axis and looks identical to a noisy one. The ranges are wide -# enough for any configuration this driver offers -- cycle counts 50 to 800, so -# rates from ~2 to ~530 Hz and LSBs from 3 to 50 nT. Where an axis is -# logarithmic the range is whole decades, so the tick labels are powers of ten. +# Every panel is drawn on a fixed scale, so two runs can be laid side by side +# and compared by eye. The ranges cover cycle counts 50 to 800, and a logarithmic +# axis spans whole decades so the ticks are powers of ten. # -# The ASD frequency axis is the exception: linear, not logarithmic. A log axis -# gives the sub-hertz drift most of the width and crams every line worth naming -# -- mains, its harmonics, fs/4 -- into the last fifth. Linear spaces them -# evenly, which is how a spectrum is read for lines. It costs the low-frequency -# decades, which the Allan panel beside it covers better anyway, and it cuts off -# above 100 Hz: at cycle count 100 that leaves out 100-150 Hz, where nothing but -# the noise floor has ever appeared. +# The ASD frequency axis is linear, which spaces mains, its harmonics and fs/4 +# evenly. The Allan panel beside it covers the low-frequency decades. ASD_HZ = (0.0, 100.0) ASD_NT = (5e-1, 5e1) # floor ~1 nT/rtHz, lines to ~15 ALLAN_TAU_S = (1e-3, 1e3) # one sample period to a 1000 s average @@ -112,134 +84,80 @@ LATENCY_MS = (0.0, 20.0) # covers 1.9 ms at cc=50 to 13.7 ms at cc=400 LATENCY_COUNT = (0.8, 1e5) # compare.py plots the same two quantities normalised by |B|, so they get their -# own fixed ranges in ppm. Same reasoning: two invocations should overlay. Its -# frequency axis stays logarithmic: that figure overlays captures whose rates -# differ several-fold, and log is what puts them on comparable footing. +# own fixed ranges in ppm. Its frequency axis is logarithmic, which is what puts +# captures whose rates differ several-fold on comparable footing. FRACTIONAL_ASD_PPM = (1e0, 1e3) FRACTIONAL_ALLAN_PPM = (1e0, 1e3) FRACTIONAL_ASD_HZ = (1e-2, 1e3) -# Spectrogram segment length. Longer than SEGMENT_SECONDS because the job here -# is to resolve one line from another -- 40 s gives 0.025 Hz, enough to separate -# a mains line from a sample-locked one even when they land close together -- -# and a ~300 s capture still yields ~100 columns at this overlap. +# Spectrogram segment length. The job here is to resolve one line from another, +# and 40 s gives 0.025 Hz, enough to separate a mains line from a sample-locked +# one. 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. +# Overlap buys columns, so a feature lasting a few seconds spans several of them. 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. +# Two things let captures share one colour scale. They are drawn in ppm of |B| +# per sqrt(Hz), which divides out the gain, and SPECTROGRAM_BIN_HZ averages +# native bins onto one displayed step. The window is a fixed duration, so that +# step gives every capture the same dof per cell. `--max-freq` pins the frequency +# axis too, at the cost of blank space on the slower captures. 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). +# Turbo by default: a spectrogram is read for narrow lines against a broad +# floor, and its hue steps make those legible. `viridis` is the +# lightness-monotonic alternative, `sequential` the single-hue ramp below. 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. +# the spectrum is drift, which varies far more between runs than the floor does. 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. +# Upper edge, as a fraction of Nyquist. Half of it is held below +# DECIMATE_CUTOFF_FRACTION, the corner of the anti-alias filter, so that every +# compared path is scored on its flat band. The assert below checks it. 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. +# Welch segment duration. Frequency resolution is 1/this, so 20 s resolves +# 0.05 Hz 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. +# the statistic is Rayleigh(1), so 5 is a one-in-270,000 bin across the ~20 +# fractions tested. Real artifacts here reach 30. LINE_SIGMA = 5.0 # Drift figure. The left panel fits the rate on windows of DRIFT_WINDOW_S; the -# right sweeps that window length to show what it costs. 10 s is a compromise -# and not the optimum -- sigma_y bottoms out nearer 2 s on these parts -- but 2 s -# windows put 145 points on a 290 s capture, which draws as a band rather than a -# trend. At 10 s there are ~29, the floor is still only 13-34 ppm, and the drift -# being plotted is 186-1241 ppm, so the error bars stay a minor part of the span. +# right sweeps that window length to show what it costs. 10 s puts ~29 points on +# a 300 s capture, with a sigma_y floor small against the drift being plotted. DRIFT_WINDOW_S = 10.0 -# Log-spaced, from a window that holds a few hundred samples to a third of a -# 300 s capture -- past that there are fewer than three windows and `rate_allan` -# drops the point rather than estimate a deviation from one difference. +# Log-spaced, from a window holding a few hundred samples to a third of a 300 s +# capture, which is as far as `rate_allan` finds enough windows. RATE_TAUS = (1.0, 1.5, 2.0, 3.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 50.0, 70.0, 100.0) -DRIFT_PPM = (-800.0, 800.0) # the worst run here spans 1241 ppm peak to peak +DRIFT_PPM = (-800.0, 800.0) RATE_TAU_S = (1.0, 1e2) RATE_ALLAN_PPM = (1e0, 1e3) +# Sub-windows used by rate_stability. Eight over a five-minute run is ~40 s +# each, long enough to fit a rate against and short enough to see a thermal +# transient. +STABILITY_WINDOWS = 8 + # Time-base figure. The host clock is fitted against sample index by a # polynomial of this degree: high enough to follow a warm-up transient, low -# enough that it cannot chase the per-read scheduling jitter, which carries no -# information about the chip. Degree 1 is exactly the uniform grid already in -# use, and nothing above 3 changed any measured amplitude on these captures -- -# consistent with `calibrate.rate_allan`, which finds nothing above the noise -# below a ~2 s timescale. +# enough to leave the per-read scheduling jitter behind. Degree 1 is the uniform +# grid already in use, and nothing above 3 changed a measured amplitude here. TIMEBASE_DEGREE = 3 -# Half-width of the scan around each line, and how finely it is sampled. A -# 290 s capture resolves 1/T = 0.0034 Hz, so +/-0.05 Hz is ~30 resolution widths -# -- wide enough to show the shape of a peak and whether it moved, narrow enough -# that the peak is not a single pixel. +# Half-width of the scan around each line, and how finely it is sampled: ~30 +# resolution widths on a 300 s capture, so a peak has a visible shape. TIMEBASE_SPAN_HZ = 0.05 TIMEBASE_POINTS = 401 TIMEBASE_PANELS = 3 -# A verdict needs a peak that is actually a peak: this is how far the scan's -# maximum must stand above the scan's own median before the ratio of the two -# time bases means anything. -# -# Local contrast rather than significance against the broadband noise, because -# that is the question. A line can sit at a respectable sigma across the whole -# band and still have no resolvable peak in a 0.1 Hz window -- the 60 Hz alias -# in |B| at cycle count 400 does exactly that, ranking third at 9 sigma while -# its scan is visibly flat -- and comparing where two flat curves happen to peak -# is comparing noise. 3x is a little above what a smooth curve reaches by -# chance and well below the 4-6x a real line here shows. +# Local contrast the scan's maximum must reach against the scan's own median +# before the ratio of the two time bases carries a verdict. TIMEBASE_MIN_CONTRAST = 3.0 @@ -248,26 +166,41 @@ def band_for(fs): return BAND_LO_HZ, BAND_NYQUIST_FRACTION * fs / 2 +def spec_asd_nt(cycle_count, fs): + """Table 3-1's noise as a spectral density at the rate actually sampled. + + Table 3-1 gives two figures. This uses the per-sample noise column, which + `white_sd` is directly comparable against: white noise of sd + `expected_noise_nt` sampled at fs has a one-sided amplitude spectral density + of sd * sqrt(2 / fs). Its other column, a 1.2 nT/sqrt(Hz) noise density, is + quoted at three times this rate. + """ + return rm3100.expected_noise_nt(cycle_count) * math.sqrt(2.0 / fs) + + +def _power_segments(v, fs, nperseg, step): + """Yield (start, one-sided PSD) for each overlapping segment of v. + + Shared by welch_asd and spectrogram, so both get the same Hann window, the + same one-sided scaling and the same linear detrend, which takes out DC and + any slow ramp before it smears across the low-frequency bins. + """ + window = np.hanning(nperseg) + scale = 1.0 / (fs * (window ** 2).sum()) + index = np.arange(nperseg) + for start in range(0, len(v) - nperseg + 1, step): + 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 # fold negative frequencies + yield start, spectrum + + 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) + psd = np.mean([s for _, s in _power_segments(v, fs, nperseg, nperseg // 2)], + axis=0) freqs = np.fft.rfftfreq(nperseg, 1.0 / fs) return freqs[1:], np.sqrt(psd[1:]) # drop DC bin @@ -291,17 +224,84 @@ def allan_deviation(v, fs, points=40): return np.array(taus), np.array(devs) +def window_rates(cap, seconds): + """Sample rate fitted independently on each consecutive window of a capture. + + Returns (centres, rates), in seconds from the start of the capture and in + Hz, with one point per window that held enough samples to fit. + + `Capture.dt_true` fits one straight line through the whole run. This repeats + that fit on pieces short enough that the rate barely moves within any one of + them, which makes the drift itself measurable. + """ + seconds = float(seconds) + start = float(cap.elapsed[0]) + if seconds <= 0 or cap.duration < 2 * seconds: + return np.array([]), np.array([]) + centres, rates = [], [] + for w in range(int(cap.duration // seconds)): + low = start + w * seconds + try: + rate = cap.restrict(low, low + seconds).true_rate_hz + except capture.CaptureError: + continue + rates.append(rate) + centres.append((w + 0.5) * seconds) + return np.array(centres), np.array(rates) + + +def rate_stability(cap, windows=STABILITY_WINDOWS): + """Fractional spread of the sample rate across sub-windows of a capture. + + The chip's RC oscillator warms up and every capture drifts downward through + its own run, so a single fitted `dt_true` is the average over that transient. + This measures how far the rate moved across it: the error bar on `dt_true`, + and the input to every uncertainty downstream. + """ + if len(cap.sample_index) < 64 * windows: + return 0.0 + _, rates = window_rates(cap, cap.duration / windows) + if len(rates) < 2: + return 0.0 + return float((rates.max() - rates.min()) / rates.mean()) + + +def rate_allan(cap, taus): + """Two-sample (Allan) deviation of the measured sample rate, fractional. + + For each averaging time tau the rate is fitted on consecutive tau-second + windows and the deviation taken between neighbours: + + sigma_y(tau) = sqrt( mean( (y[k+1] - y[k])^2 ) / 2 ) + + Falling with tau is measurement noise averaging down; rising with tau is + drift the measurement has begun to resolve, and a rise proportional to tau + marks a deterministic frequency ramp. The minimum is the best the rate can + be known, at the averaging time worth using. + + Returns (taus, sigma_y) over the taus that yielded three or more windows. + """ + kept, devs = [], [] + for tau in taus: + _, rates = window_rates(cap, tau) + if len(rates) < 3: + continue + y = rates / rates.mean() + kept.append(float(tau)) + devs.append(float(np.sqrt(np.mean(np.diff(y) ** 2) / 2))) + return np.array(kept), 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. + sqrt(2) recovers the original sd, while any component slower than the sample + rate is differenced away. - 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. + Differencing has a +6 dB/octave response, so a narrowband line near Nyquist + is weighted about twice the broadband floor. `band_stats` gives a figure + that is flat across the band. """ if len(v) < 2: return 0.0 @@ -311,12 +311,9 @@ def white_sd(v): 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. + A fixed duration gives two captures at different rates the same number of + segments to average, which is what makes their Welch medians comparable: that + median is biased low by an amount set by the segment count. """ return max(64, int(round(fs * seconds))) @@ -324,9 +321,9 @@ def segment_length(fs, seconds=SEGMENT_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. + The median is the floor, robust to the handful of bins a line occupies. The + RMS integrates the whole band, lines included, so the two together say how + much of the band is line. """ lo, hi = band if band is not None else band_for(fs) freqs, asd = welch_asd(v, fs, nperseg=nperseg or segment_length(fs)) @@ -342,13 +339,13 @@ def band_stats(v, fs, band=None, nperseg=None): 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. + `cutoff` is in cycles per sample (0 to 0.5). Blackman buys a -74 dB peak + sidelobe for a transition band of about 5.5/ntaps, which is the trade that + matters: what the filter stops has to stay stopped. In the decimate-by-4 + case the realised response measures -92 dB just past the corner and -135 dB + at 60 Hz, so mains stays out when a fast capture is decimated. - Normalised to unit gain at DC, so decimating does not rescale the field. + Normalised to unit gain at DC, which keeps the field at its own scale. """ if not 0 < cutoff < 0.5: raise ValueError(f"cutoff {cutoff} must be in (0, 0.5) cycles/sample") @@ -361,11 +358,11 @@ def fir_lowpass(cutoff, ntaps): # 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. +# The assert below pins it between BAND_NYQUIST_FRACTION / 2, which keeps the +# compared band on the flat part, and 0.5, which keeps the stopband out of the +# fold. That leaves a transition band of (0.5 - this)/k; Blackman needs about +# 5.5/ntaps, so 128*k+1 taps carries roughly 2x margin and is odd by +# construction. DECIMATE_CUTOFF_FRACTION = 0.45 assert BAND_NYQUIST_FRACTION / 2 < DECIMATE_CUTOFF_FRACTION < 0.5 @@ -373,15 +370,13 @@ 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: + Two filters, answering two 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. + boxcar average k consecutive samples, which is what the sensor does over + a cycle count, poor stopband and all. Answers "would a slower + cycle count have given me this?" + fir a proper anti-alias filter, which removes out-of-band + interference. Answers "what does filtering in post give me?" Returns the decimated series; the new rate is fs/k. """ @@ -403,9 +398,8 @@ def decimate(v, k, method="fir"): 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. + # 'same' pads with zeros, so the taper reaches ntaps//2 into the result. A + # full filter length off each end clears it. return np.convolve(v, h, mode="same")[len(h):-len(h)][::k] @@ -419,21 +413,18 @@ def sample_locked_lines(v, max_period=8, threshold=LINE_SIGMA): 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. + This finds *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 moves with the sampling, so it is an + artifact of the measurement. compare.py settles it across two captures. - 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. + The transform is evaluated at exactly j/period, which puts the whole capture + behind one number and 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. + once. A cubic trend is removed first, which keeps drift out of the low + fractions. """ v = np.asarray(v, dtype=float) if len(v) < 4 * max_period: @@ -449,8 +440,8 @@ def sample_locked_lines(v, max_period=8, threshold=LINE_SIGMA): 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. + # Nyquist is a real alternation with no phase to fit: one + # degree of freedom, and so a tighter null. amplitude = abs(float((v * (-1.0) ** index).mean())) sigma = amplitude / (sd / np.sqrt(len(v))) else: @@ -467,32 +458,18 @@ def sample_locked_lines(v, max_period=8, threshold=LINE_SIGMA): def spectrogram(v, fs, seconds=SPECTROGRAM_SECONDS, overlap=SPECTROGRAM_OVERLAP): - """(times, freqs, ASD) -- how the spectrum of a series evolves. + """(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. + Keeping the time axis separates a line that ran throughout from a burst. """ 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 + for start, spectrum in _power_segments(v, fs, nperseg, step): columns.append(np.sqrt(spectrum[1:])) times.append((start + nperseg / 2) / fs) freqs = np.fft.rfftfreq(nperseg, 1.0 / fs)[1:] @@ -502,11 +479,9 @@ def spectrogram(v, fs, seconds=SPECTROGRAM_SECONDS, overlap=SPECTROGRAM_OVERLAP) 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. + Averages in *power*, squaring and rooting around the mean, which is what + adds degrees of freedom. A fixed step gives every capture the same dof per + cell. """ width = max(1, int(round(step / (freqs[1] - freqs[0])))) if width == 1: @@ -526,20 +501,19 @@ 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. + samples. `sd` is the plain standard deviation of the detrended series; the + ratio is against the ideal `white_sd`/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. + This is the check the cycle-count choice rests on. Section 5.1 warns that + quantisation "generally dictates working above a cycle count of ~30": 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. That shows here as a ratio climbing 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. + A slow trend is removed first, which leaves the quantiser as the only thing + that can stop the average shrinking. """ v = np.asarray(v, dtype=float) index = np.arange(len(v)) @@ -558,14 +532,38 @@ def dither_check(v, lsb, longest=1024): return out +def style(ax, grid=True): + """Apply the shared light-mode look to one axes. Used by compare.py too.""" + ax.set_facecolor(SURFACE) + if grid: + 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) + + +def title(ax, text, size=12, pad=8): + ax.set_title(text, loc="left", color=TEXT_PRIMARY, fontsize=size, + fontweight="bold", pad=pad) + + +def caption(cap): + """The one-line 'file, N samples at F Hz, cycle count C' figure header.""" + return (f"{cap.path.rsplit('/', 1)[-1]}: {len(cap.sample_index):,} samples " + f"at {cap.true_rate_hz:.2f} Hz, cycle count {cap.cycle_count}") + + 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. + than a decade, producing a dense scale reading 2x10^0, 3x10^0. Widening the + range to a whole decade restores the convention. """ axes = {"x": [ax.xaxis], "y": [ax.yaxis], "both": [ax.xaxis, ax.yaxis]}[which] @@ -586,10 +584,8 @@ def decade_ticks(ax, which="both", minimum_decades=1.0): 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. + Mains is drawn where it lands. Above Nyquist it has folded, so the line + goes at the folded position and carries an "alias" label. """ marks = [(fs / 4, "fs/4", False), (fs / 2, "fs/2 = Nyquist", False)] for mains in (60.0, 120.0): @@ -604,20 +600,16 @@ def reference_lines(fs): def coherent_amplitude(v, t, freqs): """Amplitude of a tone at each of `freqs`, for samples taken at times `t`. - Evaluates 2*|mean(v * exp(-2*pi*i*f*t))| directly. That is what a spectrum - computes, but without requiring the samples to be evenly spaced, which is - the entire reason this exists rather than a call to `welch_asd`: the - question here is which *time base* a line is coherent on, and one of the two - candidates is not uniform. + Evaluates 2*|mean(v * exp(-2*pi*i*f*t))| directly, which is what a spectrum + computes but accepts unevenly spaced samples. One of the two time bases + compared here is non-uniform, so it needs that. - A tone present at f throughout the record sums in phase and returns its own - amplitude; anything else averages towards zero as 1/sqrt(n). A sample placed - at the wrong time by a fraction d of a cycle contributes cos(2*pi*d) instead - of 1, so coherence is lost at a rate set by the tone's frequency -- which - makes this the sharpest test available of whether a time base is right, and - also the reason a low-frequency line can say nothing about one. + A sample placed at the wrong time by a fraction d of a cycle contributes + cos(2*pi*d), so coherence is lost at a rate set by the tone's frequency. + That makes this the sharpest available test of a time base, on a line high + enough in frequency to lose coherence within the record. - The mean is removed first, or a DC offset leaks into the lowest frequencies. + The mean is removed first, keeping a DC offset out of the lowest frequencies. """ v = np.asarray(v, dtype=float) v = v - v.mean() @@ -629,22 +621,17 @@ def coherent_amplitude(v, t, freqs): def host_time_base(cap, degree=TIMEBASE_DEGREE): """Sample times taken from the host clock, smoothed against sample index. - The alternative to `cap.elapsed`, which is `sample_index * dt_true` and so - asserts the rate never moved. It did: `calibrate.rate_stability` measures - 186 to 1241 ppm of movement across the captures here. + `cap.elapsed` is `sample_index * dt_true`, a rate held constant across the + run. The rate in fact moves by the hundreds of ppm `rate_stability` + measures, and raw host timestamps carry that movement. - Raw host timestamps carry that drift, but also the scheduling jitter of each - individual read, which says nothing about when the chip actually sampled. - A low-order polynomial keeps the first and discards the second. Degree 1 - reduces exactly to `cap.elapsed`, which makes the uniform grid the degree-1 - member of this family rather than a separate idea, and nothing above degree - 3 changed any measured amplitude on these captures. + They also carry the scheduling jitter of each read. A low-order polynomial + keeps the drift and leaves the jitter behind; degree 1 reduces to + `cap.elapsed`, making the uniform grid the degree-1 member of this family. - Returned relative to the first sample, like `cap.elapsed - cap.elapsed[0]`. - Shifting to that origin is what makes the degree-1 case *exactly* the - uniform grid rather than merely parallel to it: the fit carries a non-zero - intercept, since a least-squares line through a curve does not pass through - its first point. + Returned relative to the first sample, which is what makes that degree-1 + case *exactly* the uniform grid: a least-squares line through a curve + passes above or below its first point. """ index = cap.sample_index - cap.sample_index[0] host = cap.system_time - cap.system_time[0] @@ -673,8 +660,7 @@ def make_spectrogram(cap, path, seconds=SPECTROGRAM_SECONDS, fs = cap.true_rate_hz times, freqs, asd = spectrogram(cap.total, fs, seconds, overlap) native = freqs[1] - freqs[0] - import calibrate - rate_drift = calibrate.rate_stability(cap) + rate_drift = rate_stability(cap) freqs, asd = bin_frequency(freqs, asd) # Fractional, so a gain difference between captures is not read as a noise # difference when they share a colour scale. @@ -685,15 +671,13 @@ def make_spectrogram(cap, path, seconds=SPECTROGRAM_SECONDS, fig, ax = plt.subplots(figsize=(13.5, 7), dpi=150) fig.patch.set_facecolor(SURFACE) - ax.set_facecolor(SURFACE) + style(ax, grid=False) 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. + # imshow suits a grid regular in both axes and draws thousands of times + # fewer objects. 'antialiased' resamples when there are more rows than + # pixels, which keeps a one-bin line visible wherever it falls. image = ax.imshow(asd, origin="lower", aspect="auto", cmap=ramp, norm=mcolors.LogNorm(*SPECTROGRAM_PPM), extent=[times[0], times[-1], freqs[0], freqs[-1]], @@ -704,17 +688,14 @@ def make_spectrogram(cap, path, seconds=SPECTROGRAM_SECONDS, 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. + # clipped by the colorbar. Anything at the very top is the axis boundary + # itself, 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. + # Lines can land within a label's height of each other, so labels are + # staggered across the width, which leaves both lines in place. ax.annotate(f"{label} {hz:.2f} Hz", xy=(0.008 + 0.17 * (slot % 3), hz), xycoords=("axes fraction", "data"), @@ -723,34 +704,25 @@ def make_spectrogram(cap, path, seconds=SPECTROGRAM_SECONDS, 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) + title(ax, "Spectral density of |B| over time", size=13) 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. + # The blank above Nyquist is real missing bandwidth, and marking it + # keeps that difference between configurations visible. 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. + # Columns far outnumber the independent windows, because they overlap, so + # cell-to-cell scatter is the two degrees of freedom of a single + # periodogram. What is real 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}", + fig.text(0.5, 0.965, caption(cap), color=TEXT_PRIMARY, fontsize=10, ha="center") fig.text(0.5, 0.935, - f"oscillator drifted {rate_drift * 1e6:.0f} ppm across the run — " + f"oscillator drifted {rate_drift * 1e6:.0f} ppm across the run, " f"a line at fs/4 smears by {rate_drift * fs / 4:.4f} Hz. " f"{seconds:g} s windows, {native:.3f} Hz native binned to " f"{SPECTROGRAM_BIN_HZ:g} Hz ({dof} dof/cell), {overlap:.2%} overlap, " @@ -769,31 +741,25 @@ def make_drift(cap, path, window=DRIFT_WINDOW_S): Two views of one measurement. The left panel is the chip's rate through the run; the right is how well that rate can be known as a function of how long - it is measured for -- which is what decides whether the left panel shows a - signal or its own noise. On these parts it is emphatically a signal: 186 to - 1241 ppm of drift against a floor near 10 ppm. + it is measured for, which sets the error bar on the left. """ - import calibrate - centres, rates = calibrate.window_rates(cap, window) + centres, rates = window_rates(cap, window) if len(rates) < 3: raise ValueError(f"{cap.duration:.0f} s gives {len(rates)} window(s) of " f"{window:g} s, too few to show a trend") - taus, devs = calibrate.rate_allan(cap, RATE_TAUS) + taus, devs = rate_allan(cap, RATE_TAUS) if not len(taus): raise ValueError(f"{cap.duration:.0f} s is too short for any averaging " f"time in {RATE_TAUS[0]:g}-{RATE_TAUS[-1]:g} s") - # The error bar is the two-sample deviation at the window length actually - # plotted. Not the scatter of the points themselves, which is the drift this - # figure exists to show, and using it would beg the question. + # The error bar is the two-sample deviation at the window length plotted, + # which is independent of the point scatter this figure exists to show. floor = float(np.interp(window, taus, devs)) best = int(np.argmin(devs)) ppm = (rates / rates.mean() - 1.0) * 1e6 - span = calibrate.rate_stability(cap) * 1e6 + span = rate_stability(cap) * 1e6 fig, axs = plt.subplots(1, 2, figsize=(13.5, 5.5), dpi=150) fig.patch.set_facecolor(SURFACE) - for ax in axs: - ax.set_facecolor(SURFACE) a = axs[0] a.axhline(0.0, color=REFERENCE, linestyle="--", linewidth=1.0) @@ -803,16 +769,15 @@ def make_drift(cap, path, window=DRIFT_WINDOW_S): trend = np.polyfit(centres, ppm, 1) a.plot(centres, np.polyval(trend, centres), color=TEXT_PRIMARY, linewidth=1.5, label=f"trend {trend[0] * cap.duration:+,.0f} ppm across the run") - a.set_title("Chip oscillator against the host clock", loc="left", - color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) + title(a, "Chip oscillator against the host clock") a.set_xlabel("elapsed (s)") a.set_ylabel("rate error (ppm of the run mean)") # The y axis is fixed so two runs can be compared; the x axis is the # capture's own length, which pinning would simply clip on a longer one. a.set_ylim(*DRIFT_PPM) - a.annotate("the sample grid stays uniform in index —\nthis is its spacing " - "in seconds moving.\nGain does not follow: it is set by the " - "ratio\nof two on-die clocks, which drift together.", + a.annotate("the sample grid stays uniform in index;\nthis is its spacing " + "in seconds moving.\nGain holds: regressing |B| on rate\ngives " + "+0.08 to +0.35, where 1/f needs −1.", xy=(0.02, 0.05), xycoords="axes fraction", color=TEXT_SECONDARY, fontsize=9) @@ -832,9 +797,8 @@ def make_drift(cap, path, window=DRIFT_WINDOW_S): rising = taus[best:] a.plot(rising, devs[-1] * 1e6 * rising / taus[-1], linestyle=":", color=REFERENCE, linewidth=1.5, - label="∝ τ — a deterministic frequency ramp") - a.set_title("How well the rate can be measured", loc="left", - color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) + label="∝ τ, a deterministic frequency ramp") + title(a, "How well the rate can be measured") a.set_xlabel("averaging time τ (s)"); a.set_ylabel("σ_y (ppm)") a.set_xlim(*RATE_TAU_S); a.set_ylim(*RATE_ALLAN_PPM) decade_ticks(a) @@ -844,33 +808,21 @@ def make_drift(cap, path, window=DRIFT_WINDOW_S): color=TEXT_SECONDARY, fontsize=9) for ax in axs: - ax.grid(True, which="both", color=GRID, linewidth=0.7) - ax.set_axisbelow(True) - for side in ("top", "right"): - ax.spines[side].set_visible(False) - for side in ("left", "bottom"): - ax.spines[side].set_color(GRID) - ax.tick_params(colors=TEXT_SECONDARY, labelsize=9, length=0) - ax.xaxis.label.set_color(TEXT_SECONDARY) - ax.yaxis.label.set_color(TEXT_SECONDARY) + style(ax) ax.legend(frameon=False, fontsize=9, labelcolor=TEXT_SECONDARY, loc="upper right") - fig.text(0.5, 0.955, - f"{cap.path.rsplit('/', 1)[-1]} — {len(cap.sample_index):,} samples " - f"at {cap.true_rate_hz:.2f} Hz, cycle count {cap.cycle_count}", + fig.text(0.5, 0.955, caption(cap), color=TEXT_PRIMARY, fontsize=10, ha="center") - # Two figures for the drift, deliberately: the peak-to-peak spread is what - # `rate_stability` feeds to the uncertainty propagation and what the rest of - # the tooling quotes, while the fitted trend uses every window and is the - # better estimate of the warm-up itself. They differ by how much of the - # transient a coarse window averages over, and printing only one invites the - # other to be read off the legend and mistaken for a discrepancy. + # Both drift figures are quoted deliberately: the peak-to-peak spread is + # what `rate_stability` feeds to the uncertainty propagation, the fitted + # trend uses every window and is the better estimate of the warm-up. Showing + # one invites the other to be read off the legend as a discrepancy. fig.text(0.5, 0.915, f"drifted {span:.0f} ppm peak to peak, " f"{abs(trend[0] * cap.duration):.0f} ppm by the fitted trend, " f"against a floor of {devs[best] * 1e6:.1f} ppm at τ = " - f"{taus[best]:g} s — resolved {span / (devs[best] * 1e6):.0f}× over.", + f"{taus[best]:g} s, resolved {span / (devs[best] * 1e6):.0f}× over.", color=TEXT_SECONDARY, fontsize=9, ha="center") fig.tight_layout(rect=[0, 0, 1, 0.90]) fig.savefig(path, facecolor=SURFACE) @@ -884,27 +836,23 @@ def make_time_base(cap, path): A line at a fixed frequency lives in wall time, so it sharpens when the samples are placed by the host clock. A line locked to the sampling lives in sample index, and the same substitution wrecks it. Measuring both ways - therefore separates a real signal from an artefact of measuring *within a + therefore separates a real signal from an artifact of measuring *within a single capture*, where compare.py needs two at different rates to do it. - The catch is that the test only has force on a strong line: coherence is - lost in proportion to the tone's frequency, so a low-frequency or aliased - line says little, and a weak one returns noise. Both are marked as such - rather than being given a verdict they cannot support. + The test has force on a strong line, coherence being lost in proportion to + the tone's frequency. A low-frequency or aliased line says little and a weak + one returns noise, and both are labelled as such. """ fs = cap.true_rate_hz chip_t = (cap.sample_index - cap.sample_index[0]) * cap.dt_true host_t = host_time_base(cap) - # The same cubic detrend `sample_locked_lines` applies, and for the same - # reason: drift is not a line, and left in place it leaks into everything - # below a few tenths of a hertz. + # The same cubic detrend `sample_locked_lines` applies, which keeps drift + # out of everything below a few tenths of a hertz. order = np.arange(len(cap.total)) v = cap.total - np.polyval(np.polyfit(order, cap.total, 3), order) - # Against the quantiser rather than against zero. A capture whose |B| never - # moves by even a thousandth of an LSB has no spectrum -- every sample is - # the same integer -- but the cubic detrend still leaves ~1e-11 nT of - # floating-point residue, so an exact test would let it through and then - # find "lines" in rounding error. + # Against the quantiser: on a capture whose |B| never moves the detrend + # still leaves ~1e-11 nT of floating-point residue, and an exact test would + # find "lines" in it. if v.std() < cap.lsb_nt * 1e-3: raise ValueError("|B| is constant to within a thousandth of an LSB; " "no lines to test") @@ -928,14 +876,14 @@ def make_time_base(cap, path): if not scored: raise ValueError("no lines to compare time bases on") - # Floored rather than purely per-panel: a capture with one or two lines - # worth testing would otherwise be narrower than its own caption. + # Floored at 9.8 in, which keeps a one- or two-panel figure wider than its + # own caption. fig, axs = plt.subplots(1, len(scored), figsize=(max(9.8, 4.7 * len(scored)), 5.2), dpi=150, squeeze=False) fig.patch.set_facecolor(SURFACE) for ax, (sigma, hz, label, folded) in zip(axs[0], scored): - ax.set_facecolor(SURFACE) + style(ax) freqs = np.linspace(hz - TIMEBASE_SPAN_HZ, hz + TIMEBASE_SPAN_HZ, TIMEBASE_POINTS) chip = coherent_amplitude(v, chip_t, freqs) @@ -943,24 +891,23 @@ def make_time_base(cap, path): ratio = host.max() / chip.max() contrast = chip.max() / np.median(chip) if folded: - verdict = "aliased — folding assumes uniformity, not a clean test" + verdict = "aliased: folding already assumes a uniform grid" elif contrast < TIMEBASE_MIN_CONTRAST: - verdict = f"no resolved peak ({contrast:.1f}× its floor) — no verdict" + verdict = f"no resolved peak ({contrast:.1f}× its floor), no verdict" elif ratio > 1.05: - verdict = "coherent in wall time — external" + verdict = "coherent in wall time: external" elif ratio < 0.95: - verdict = "coherent in sample index — an artefact of measuring" + verdict = "coherent in sample index: an artifact of measuring" else: verdict = "no preference between the two" ax.plot(freqs - hz, chip, color=AXES[0][2], linewidth=1.3, - label=f"chip index grid — {chip.max():.3f} nT") + label=f"chip index grid, {chip.max():.3f} nT") ax.plot(freqs - hz, host, color=AXES[1][2], linewidth=1.3, - label=f"host clock, degree {TIMEBASE_DEGREE} — {host.max():.3f} nT") + label=f"host clock, degree {TIMEBASE_DEGREE}, {host.max():.3f} nT") # Padded so the verdict can sit between title and axes without either # landing on the traces, which reach the top of a tightly scaled panel. - ax.set_title(f"{label} at {hz:.3f} Hz {ratio:.2f}×", loc="left", - color=TEXT_PRIMARY, fontsize=11, fontweight="bold", pad=22) + title(ax, f"{label} at {hz:.3f} Hz {ratio:.2f}×", size=11, pad=22) ax.annotate(verdict, xy=(0.0, 1.01), xycoords="axes fraction", color=TEXT_SECONDARY, fontsize=8.5, va="bottom") ax.set_xlabel(f"offset from {hz:.3f} Hz (Hz)") @@ -969,15 +916,6 @@ def make_time_base(cap, path): # Headroom for the legend, which would otherwise sit on the peak the # panel exists to show. ax.set_ylim(0.0, max(chip.max(), host.max()) * 1.45) - ax.grid(True, which="both", color=GRID, linewidth=0.7) - ax.set_axisbelow(True) - for side in ("top", "right"): - ax.spines[side].set_visible(False) - for side in ("left", "bottom"): - ax.spines[side].set_color(GRID) - ax.tick_params(colors=TEXT_SECONDARY, labelsize=9, length=0) - ax.xaxis.label.set_color(TEXT_SECONDARY) - ax.yaxis.label.set_color(TEXT_SECONDARY) ax.legend(frameon=False, fontsize=8.5, labelcolor=TEXT_SECONDARY, loc="upper right") @@ -988,8 +926,8 @@ def make_time_base(cap, path): f"cycle count {cap.cycle_count}", color=TEXT_SECONDARY, fontsize=9, ha="center") fig.text(0.5, 0.897, - "Ratio is host peak over chip peak. Diagnostic only — no other " - "figure uses the host time base.", + "Ratio is host peak over chip peak. Diagnostic only: every other " + "figure uses the chip's grid.", color=TEXT_SECONDARY, fontsize=9, ha="center") fig.tight_layout(rect=[0, 0, 1, 0.88]) fig.savefig(path, facecolor=SURFACE) @@ -1000,14 +938,10 @@ def make_time_base(cap, 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. + Never silent: the note says what went. Never automatic either, so a figure + always describes the capture the caller named. Below 4x the requested trim + the capture is returned whole with a note, which keeps a short deliberate + capture analysable. """ if not seconds: return cap, "" @@ -1032,7 +966,7 @@ def main(): 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 " + 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", @@ -1055,8 +989,8 @@ def main(): "(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. + # --trim and --start/--end both choose a window, so asking for both is an + # error and the reported figures stay independent of 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") @@ -1071,43 +1005,31 @@ def main(): 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. + # Total field is derived from the three axes, so its noise tracks whichever + # axis carries most of the field. 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 + fs = cap.true_rate_hz # measured, never the nominal table value 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 + latency = np.diff(cap.system_time) # read-to-read spacing # 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. + # differs: a gain change carries the noise with it, so absolute sd tracks + # the scale as well as the noise. field = data["total"].mean() band = band_for(fs) - # Table 3-1's noise figure for this cycle count -- what `white sd` is the - # measurement of, so the two belong in the same table. spec_nt = rm3100.expected_noise_nt(cap.cycle_count) print(cap.summary()) if trim_note: print(f" {trim_note}") - # How far the oscillator moved during the run. Imported here rather than at - # module scope: calibrate.py imports this module, so the dependency only - # runs one way at import time. - import calibrate - drift = calibrate.rate_stability(cap) - # Reported, not flagged. The sensor oscillator is an LR relaxation circuit - # -- coil inductance, a resistor and comparator thresholds, none of them - # temperature-stable -- so a tenth of a percent across a run is ordinary, - # not a fault. What it is, is the error bar on this capture's rate and on - # anything derived from it. + # Reported without a threshold. The sensor oscillator is an LR relaxation + # circuit, none of it temperature-stable, so a tenth of a percent across a + # run is ordinary. It is the error bar on this capture's rate. + drift = rate_stability(cap) print(f" oscillator drifted {drift * 1e6:.0f} ppm across the run " f"({drift * cap.dt_true * 1e6:.2f} us on a " - f"{cap.dt_true * 1e3:.2f} ms period) -- the error bar on the rate") + f"{cap.dt_true * 1e3:.2f} ms period), the error bar on the rate") print() print(f"{'series':6s} {'sd':>9s} {'white sd':>9s} {'dither':>7s} " f"{'sd/|B|':>10s} {'p2p':>9s} {'vs spec':>8s} {'ASD':>7s}") @@ -1115,30 +1037,28 @@ def main(): 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) + style(ax) 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. + # Drift-free: what Table 3-1 quotes, and what the dither margin is made + # of. The plain sd beside it is drift on a drifting capture. wsd = white_sd(v) freqs, asd = welch_asd(v, fs) - # Thin: four traces overlaid on a linear 0-100 Hz axis, and the detail - # between the lines is what the panel is read for. The Allan panel below - # keeps a heavier stroke -- four smooth curves, not four dense ones. + # Thin, because the detail between the lines is what the panel is read + # for. The Allan panel keeps a heavier stroke for its smooth curves. axs[0, 0].semilogy(freqs, asd, color=color, linewidth=0.7, label=label, alpha=0.85) taus, devs = allan_deviation(v, fs) axs[0, 1].loglog(taus, devs, color=color, linewidth=1.6, label=label) - # Detrend before the histogram so slow drift does not masquerade as - # a fat tail. Bin edges are locked to quarter-LSB steps: the residual - # of an integer-valued signal is a comb at one LSB by construction, and - # arbitrary bin edges beat against it into ragged spikes that look like - # structure. Aligned bins render the comb as what it is. + # 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, and arbitrary edges + # beat against it into ragged spikes that look like structure. 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 @@ -1157,24 +1077,26 @@ def main(): 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? + # The check the cycle-count choice rests on: 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 -> " + print(f" {label:6s} {steps[0][2]:.2f} LSB per sample (plain sd) -> " 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.") + print(" A falling LSB column means sub-LSB resolution is being " + "recovered. A ratio climbing\n while it stalls is the quantiser " + "losing its dither (section 5.1); a ratio climbing\n while it " + "falls is drift, which stays put under averaging.") + print(" The per-sample figure is plain sd, holding whatever the detrend " + "left. The ratio is\n referenced to white sd, the dither column of " + "the table above.") - # 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. + # Reported per axis because this structure is not isotropic. lines = [(label, sample_locked_lines(data[key])) for key, label, _ in SERIES] if any(found for _, found in lines): @@ -1185,23 +1107,21 @@ def main(): 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.") + print(" A sample-locked line moves with the rate, which makes it an " + "artifact of\n measuring. compare.py decides it, by checking " + "whether these hold the same\n fraction of fs in a capture " + "taken at a different rate.") a = axs[0, 0] # The band the quoted median comes from, so the number and the picture # cannot drift apart. Drawn under the traces, hence the low zorder. a.axvspan(band[0], band[1], color=GRID, alpha=0.55, zorder=0, linewidth=0) - # Named in the legend rather than annotated on the line. The spec sits low - # in the range, where on a linear axis the traces are densest, and floating - # text disappears into them wherever it is put. - a.axhline(SPEC_ASD_NT, color=REFERENCE, linestyle="--", linewidth=1.2, - label=f"Table 3-1 spec, {SPEC_ASD_NT} nT/√Hz") - a.set_title("Amplitude spectral density", loc="left", - color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) + # Named in the legend, since the spec line sits low in the range where the + # traces are densest and floating text would disappear into them. + spec_asd = spec_asd_nt(cap.cycle_count, fs) + a.axhline(spec_asd, color=REFERENCE, linestyle="--", linewidth=1.2, + label=f"Table 3-1 noise at {fs:.0f} Hz, {spec_asd:.2f} nT/√Hz") + title(a, "Amplitude spectral density") 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) @@ -1210,8 +1130,7 @@ def main(): decade_ticks(a, which="y") a = axs[0, 1] - a.set_title("Allan deviation", loc="left", - color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) + title(a, "Allan deviation") a.set_xlabel("averaging time τ (s)"); a.set_ylabel("σ (nT)") a.set_xlim(*ALLAN_TAU_S); a.set_ylim(*ALLAN_NT) decade_ticks(a) @@ -1220,8 +1139,7 @@ def main(): 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) + title(a, "Residual distribution (cubic trend removed)") a.set_xlabel("nT"); a.set_ylabel("density") a.set_xlim(*RESIDUAL_NT); a.set_ylim(*RESIDUAL_DENSITY) @@ -1234,10 +1152,7 @@ def main(): a.annotate(f"grid {cap.dt_true * 1000:.2f} ms", xy=(cap.dt_true * 1000, 1), xytext=(6, 0), textcoords="offset points", color=TEXT_PRIMARY, fontsize=9) - a.set_title("Host read latency", loc="left", - color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) - # The oscillator's own drift used to be annotated here, on a panel about - # host-side comms that has nothing to do with it. It has its own figure now. + title(a, "Host read latency") a.annotate("spread is comms only;\nthe sample grid is exact", xy=(0.98, 0.94), xycoords="axes fraction", ha="right", va="top", color=TEXT_SECONDARY, fontsize=9) @@ -1245,16 +1160,6 @@ def main(): 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"), @@ -1267,7 +1172,7 @@ def main(): f"{len(t):,} samples, {cap.duration:.1f} s at {fs:.2f} Hz measured " f"(nominal {cap.nominal_rate_hz:g} Hz, {cap.rate_error * 100:+.1f}%), " f"cycle count {cap.cycle_count}, 1 LSB = {lsb_nt:.2f} nT. " - f"Drift {drift * 1e6:.0f} ppm — see the drift figure." + f"Drift {drift * 1e6:.0f} ppm, see the drift figure." + (f" {args.trim:g} s trimmed from each end." if trim_applied else ""), color=TEXT_SECONDARY, fontsize=10, ha="center") @@ -1277,11 +1182,9 @@ def main(): fig.savefig(out, facecolor=SURFACE) print(f"\n-> {out}") - # Separate files rather than more panels. The spectrogram needs the width to - # resolve anything; the other two answer questions the four panels do not -- - # whether a feature lasted, how far the chip's clock moved, and which clock - # each line is coherent on. Each degrades to a message rather than aborting - # the run, since a short capture can legitimately support none of them. + # Separate files, since the spectrogram needs the width. Each degrades to a + # message, which keeps a short capture analysable when it can support none + # of them. for suffix, build in (("spectrogram", lambda p: make_spectrogram(cap, p, args.spectrogram, args.overlap, diff --git a/compare.py b/compare.py index 0bd6882..675666e 100644 --- a/compare.py +++ b/compare.py @@ -1,5 +1,5 @@ #!/usr/bin/env python3 -"""Compare captures taken under different conditions -- e.g. two supplies. +"""Compare captures taken under different conditions, e.g. two supplies. ./.venv/bin/python compare.py a.csv b.csv c.csv ... ./.venv/bin/python compare.py --group note *.csv # label from header @@ -7,35 +7,34 @@ A label may be given as `LABEL=path`, which is the only way to identify a run whose capture was recorded without `--note`. Everything before the first `/` in -a label is the **condition** -- the thing under test -- and captures sharing one +a label is the **condition**, the thing under test, and captures sharing one are treated as the same configuration. `--supply LDO=3.0` then attaches a rail voltage to every capture in that condition. -Three lessons are built in, each learned by getting it wrong first. +Three rules are built into how it reports. **Compare |B|, not the axes.** Swapping a supply means touching the rig, and the -sensor does not go back exactly where it was. |B| survives a rotation; the -individual axes do not, so a per-axis difference measures the handling. The axes -appear here only as movement diagnostics. +sensor comes back somewhere close to where it was. |B| survives a rotation, so a +per-axis difference measures the handling. The axes appear here as movement +diagnostics. **Compare fractionally, or over a stated band.** A gain change carries the noise with it, so a run reading 6% larger also reads ~6% noisier while being physically -identical. And a broadband figure is only comparable over a band both captures -cover -- otherwise a fast capture is scored over five times the bandwidth of a -slow one. One band is chosen for the whole invocation and printed. +identical. And a broadband figure is comparable over a band both captures cover, +so one band is chosen for the whole invocation and printed. -**A magnitude difference is not automatically gain.** |B| is preserved under -rotation but *not* under translation through a field gradient, so a sensor that -moved can change magnitude on its own. Two diagnostics bound it: +**A magnitude difference can be movement.** |B| is preserved under rotation, and +translation through a field gradient changes it, so a sensor that moved can shift +magnitude on its own. Two diagnostics bound that: rotation angle the angle between mean field directions. gain residual what is left after the best single scale factor, as a fraction of |B|. Zero means one number relates the two captures, which is what a gain change looks like. -Neither can prove the sensor held still, so neither vetoes the comparison; they -size the doubt. The reproducibility floor from repeats of one condition is the -honest error bar, and this prints it when there are repeats to use. +Both size the doubt; the comparison runs either way. The reproducibility floor +from repeats of one condition is the honest error bar, and this prints it when +there are repeats to use. """ import argparse @@ -60,11 +59,9 @@ import rm3100 RATIO_SPREAD_OK = 0.01 # 1% ROTATION_OK_DEG = 0.5 -# An axis whose mean field is near zero has a ratio dominated by noise -- a few -# nT of drift on a 10 nT mean is a ratio of 1.3 with nothing moving. Such an axis -# is excluded from the spread test rather than allowed to dominate it. The -# rotation check needs no such guard: it works on the direction of the whole -# vector, where a small component simply contributes little. +# An axis whose mean field is near zero has a ratio dominated by noise, so the +# spread test drops it. The rotation check works on the whole vector's direction, +# where a small component contributes little, and needs no such guard. RATIO_AXIS_MIN_FRACTION = 0.05 # of |B| # Two captures whose rates differ by less than this are treated as the same @@ -74,9 +71,9 @@ RATE_DISTINCT = 0.01 # 1% # Condition and variant are separated by this in a label: `LDO/cc100`. CONDITION_SEP = "/" -# Hue carries the condition -- the thing under test -- and never the cycle -# count, so a chart never repaints when a run is added. Validated all-pairs in -# light mode against the #fcfcfb surface: worst CVD dE 24.7, normal-vision 33.6. +# Hue carries the condition, the thing under test, and the cycle count goes to +# line style below, so a chart keeps its colours when a run is added. Validated +# all-pairs in light mode against the #fcfcfb surface. CONDITION_COLORS = ["#2a78d6", "#eb6834", "#1baf7a"] # Cycle count is the secondary encoding, carried by line style and marker, so # identity never rests on colour alone. @@ -91,10 +88,9 @@ VARIANT_MARKERS = ["o", "s", "^", "D"] def split_label(argument): """Split `LABEL=path` into (label, path); (None, path) when unlabelled. - A path may itself contain '=', so a bare existing path always wins over the + A path may itself contain '=', so a bare existing path wins over the labelled reading. When neither reading names a file that exists, the - labelled one is returned so the error message quotes the path the user - most likely meant. + labelled one is returned, which puts the likelier path in the error message. """ label, separator, path = argument.partition("=") if not separator: @@ -160,9 +156,8 @@ def variant_of(label): def conditions_in_order(records): """Distinct conditions, in the order they were given on the command line. - Command-line order, not alphabetical: it is the one thing the caller - controls, so contrasts read in the direction they asked for and a colour - never moves between runs because a condition was renamed. + Command-line order is the one thing the caller controls, so contrasts read + in the direction they asked for and a colour stays put across runs. """ seen = [] for label, _, _ in records: @@ -186,15 +181,14 @@ def summarise(cap, band): for key, _, _ in ch.SERIES: v = d[key] median_asd, band_rms, peak, peak_hz = ch.band_stats(v, fs, band) + # The absolute figures are nanotesla, and nanotesla are what a gain + # change moves, so every cross-condition comparison uses the columns + # divided by that condition's own |B|. out[key] = { "sd": v.std(), "white": ch.white_sd(v), "ppm": v.std() / field * 1e6, "asd": median_asd, - # Every absolute figure above is in nanotesla, and nanotesla are - # what a gain change moves. A condition that reads 7% smaller reads - # 7% quieter in nT while being no quieter at all, so any comparison - # across conditions has to divide by that condition's own |B|. "asd_ppm": median_asd / field * 1e6, "white_ppm": ch.white_sd(v) / field * 1e6, "band_rms": band_rms, @@ -218,11 +212,9 @@ def axis_ratio_spread(a, b, ratio): """Spread of the per-axis scale ratios, or None if it cannot be measured. A pure gain change scales X, Y and Z identically, so the spread is the - movement test. It is only meaningful on axes that actually carry field: - dividing two near-zero means amplifies their noise without bound, and one - such axis is enough to report a stationary sensor as moved. Axes below - RATIO_AXIS_MIN_FRACTION of |B| in *either* capture are therefore dropped, - and with fewer than two left the test abstains. + movement test. Dividing two near-zero means amplifies their noise without + bound, so axes below RATIO_AXIS_MIN_FRACTION of |B| in *either* capture are + dropped, and with fewer than two left the test abstains. """ usable = [i for i in range(3) if abs(a["mean"][i]) >= RATIO_AXIS_MIN_FRACTION * a["field"] @@ -237,10 +229,10 @@ def gain_and_movement(a, b): """(scale, residual, rotation deg) relating two mean field vectors. `scale` is the single factor that best maps a onto b, and `residual` is what - that factor cannot explain, as a fraction of |b|. A pure gain change leaves - a residual of zero at any angle; a rotation leaves a large one. Unlike a - per-axis ratio this stays finite when an axis passes through zero, which is - exactly when the sensor has been re-oriented and the diagnostic is needed. + that factor leaves over, as a fraction of |b|. A pure gain change leaves a + residual of zero at any angle; a rotation leaves a large one. It stays finite + when an axis passes through zero, which is when a re-oriented sensor most + needs the diagnostic. """ ma, mb = a["mean"], b["mean"] scale = float(ma @ mb / (ma @ ma)) @@ -257,42 +249,36 @@ def fit_rate_model(caps): linear in cycle count with slope 1/C and intercept h. Two cycle counts determine both exactly; more are a least-squares fit. - Fitting it per condition is what separates an oscillator that changed from a - model that was never right: a supply that only scales the clock moves C and + Fitting it per condition separates an oscillator that changed from a model + that was wrong to begin with: a supply that only scales the clock moves C and h together by one factor, leaving the *shape* of the model alone. """ try: # The uncertainty on each period is how far the oscillator drifted # during that capture; the worst of them bounds the fit. - drift = max(calibrate.rate_stability(c) for c in caps) - fit = calibrate.fit_timing([c.cycle_count for c in caps], - [c.dt_true for c in caps], - period_sd=drift or 1e-9) - return fit + drift = max(ch.rate_stability(c) for c in caps) + return calibrate.fit_timing([c.cycle_count for c in caps], + [c.dt_true for c in caps], + period_sd=drift or 1e-9) except calibrate.CalibrationError as exc: - # One fit, in calibrate.py, so sweep.py writing a calibration and this - # reporting one cannot drift apart. The exception type is translated - # because callers here already handle ValueError. + # Translated because callers here already handle ValueError. The fit + # itself stays in calibrate.py, shared with sweep.py. raise ValueError(str(exc)) from exc def decimation_paths(fast, slow, k, band): """Noise of a fast capture decimated by k, against a natively slow one. - The question this answers is whether sampling fast and filtering down is - equivalent to having sampled slowly -- which is the whole basis for - recommending a low cycle count. Both decimations are shown because they - differ in exactly one way that matters: the boxcar is what the sensor's own - integration does, poor stopband and all, while the FIR is what a filter - designed for the job does, and only the second actually rejects what would - otherwise fold in. + Answers whether sampling fast and filtering down is equivalent to having + sampled slowly, which is the basis for recommending a low cycle count. Both + decimations are shown: the boxcar is what the sensor's own integration does, + poor stopband and all, and the FIR rejects what the boxcar folds in. - Compare the ASD column, not the sd column. Each path has a different - effective noise bandwidth, so their per-sample sd figures are not on the - same footing, while ASD over a common band is. + Compare the ASD column, which is quoted over a common band. Each path has + its own effective noise bandwidth, which the per-sample sd carries with it. """ - v_fast = dict(fast["cap"].axes())["total"] - v_slow = dict(slow["cap"].axes())["total"] + v_fast = fast["axes"]["total"] + v_slow = slow["axes"]["total"] fs_fast, fs_slow = fast["fs"], slow["fs"] rows = [("raw fast", v_fast, fs_fast, fast["field"])] for method in ("boxcar", "fir"): @@ -302,9 +288,9 @@ def decimation_paths(fast, slow, k, band): out = [] for name, v, fs, field in rows: median_asd, band_rms, _, _ = ch.band_stats(v, fs, band) + # The two cycle counts have their own gains, so the fractional column + # is the comparable one here too. out.append({"name": name, "fs": fs, "asd": median_asd, - # The two cycle counts do not share a gain either, so the - # fractional column is the comparable one here too. "asd_ppm": median_asd / field * 1e6, "band_rms": band_rms, "white": ch.white_sd(v), "lines": ch.sample_locked_lines(v)}) @@ -316,8 +302,8 @@ def decimation_pairs(entries): `entries` are (label, summary-or-capture-like) with .cycle_count and a rate. A pair qualifies when one cycle count is an integer multiple of the other, - since that is when averaging k samples of the fast run buys the same - integration time as one sample of the slow one. + which is when averaging k samples of the fast run buys the same integration + time as one sample of the slow one. """ by_condition = defaultdict(list) for label, item in entries: @@ -338,13 +324,12 @@ def decimation_pairs(entries): def comparable_rates(loaded): """Every rate that will be scored, decimated paths included. - The band has to fit inside the narrowest of these, not merely inside the - slowest capture: decimating by k lands at fs/k, and the chip's per-axis - overhead means that is a little *below* the natively-slow rate rather than - equal to it. Two percent of bandwidth is enough to score the decimated path - on its own filter rolloff and report it as the quieter one. + The band has to fit inside the narrowest of these. Decimating by k lands at + fs/k, and the per-axis overhead puts that a little *below* the natively-slow + rate; two percent of bandwidth is enough to score the decimated path on its + own filter rolloff. - Runs before the summaries exist, so it takes captures rather than records. + Runs before the summaries exist, so it takes captures. """ entries = [(label, {"fs": cap.true_rate_hz, "cycle_count": cap.cycle_count}) @@ -359,9 +344,8 @@ def locked_line_verdict(records): A field line sits at a fixed frequency, so it lands on a different fraction of fs when the rate changes. One that keeps the same fraction is locked to - the sampling and is an artefact of measuring rather than something measured. - Deciding this needs two captures at genuinely different rates, which is why - it lives here and not in characterize.py. + the sampling, and so is an artifact of measuring. Deciding this needs two + captures at genuinely different rates, which is why it lives here. """ seen = defaultdict(list) for label, _, s in records: @@ -388,7 +372,7 @@ def print_overview(records, band, supplies): for label, _, s in records: t = s["total"] print(f"{label[:20]:<20} {s['cycle_count']:5d} {s['fs']:9.3f} " - f"{calibrate.rate_stability(s['cap']) * 1e6:6.0f}p " + f"{ch.rate_stability(s['cap']) * 1e6:6.0f}p " f"{s['field']:10,.0f} {t['sd']:7.1f} {t['white']:7.2f} " f"{t['asd']:6.2f} {s['allan'][0]:6.2f} | " f"{t['white_ppm']:7.1f} {t['asd_ppm']:7.1f}") @@ -396,17 +380,16 @@ def print_overview(records, band, supplies): f"over {band[0]:g}-{band[1]:.1f} Hz.") print(" Right of the bar, the same two divided by that capture's own |B|: " "ppm and ppm/rtHz.") - print(" Compare conditions on the right. A condition whose gain differs " - "records proportionally\n smaller nanotesla of noise while being no " - "quieter, so the absolute columns would\n credit a low gain as a " - "low noise floor. Use the left against Table 3-1, which is\n " - "quoted in nT at one supply.") - print(" white = sd(diff)/sqrt(2), which rejects drift; sd does not. " + print(" Compare conditions on the right, where each capture is divided by " + "its own |B|, so a\n gain difference stays out of the noise figure. " + "Use the left against Table 3-1,\n which is quoted in nT at one " + "supply.") + print(" white = sd(diff)/sqrt(2), which rejects drift. " "Allan = best sigma by averaging.") - print(" drift is how far the oscillator moved across the run, in ppm -- the " + print(" drift is how far the oscillator moved across the run, in ppm: the " "error bar on that\n capture's rate, and on anything derived from it. " "A tenth of a percent is ordinary\n for an LR relaxation oscillator. " - "It is a timing figure: nothing magnetic reaches it.") + "It is a timing figure.") if supplies: rails = " ".join(f"{c} {v:g} V" for c, v in sorted(supplies.items())) print(f" rails: {rails}") @@ -431,8 +414,7 @@ def print_rate_model(records): f"{fit.axis_overhead_s * 1e6:6.2f} +/- " f"{fit.axis_overhead_s_sd * 1e6:4.2f} us " f"({fit.points} points, " - + ("no residual: two points fit both terms exactly, so a bad " - "point cannot show up)" + + ("exactly determined, so no residual)" if fit.points == 2 else f"residual {fit.residual_ppm:.0f} ppm)")) print(f" {'rm3100.py':<12} C {rm3100.COUNTS_PER_SECOND:10,.0f} counts/s " @@ -468,16 +450,15 @@ def print_locked_lines(records): if distinct: print(" -> SAMPLE-LOCKED. It holds this fraction of fs across " "captures whose rates\n differ, so it moves with the " - "sampling and is not a field at a fixed\n frequency. " - "Decimating a faster capture avoids it; filtering cannot.") + "sampling. Decimating a faster capture\n removes it.") elif len(hits) > 1: - print(" -> seen in several captures, but at rates too close " - "together to tell a\n sample-locked line from a fixed " - "one. Re-run with rates further apart.") + print(" -> seen in several captures, at rates too close together " + "to tell a\n sample-locked line from a fixed one. " + "Re-run with rates further apart.") else: - print(" -> seen in one capture only, so nothing here says " - "whether it is locked to\n the sampling or to a " - "frequency.") + print(" -> seen in one capture only. Deciding whether it is " + "locked to the sampling\n or to a frequency needs a " + "second capture at another rate.") def print_decimation(records, band): @@ -503,27 +484,25 @@ def print_decimation(records, band): raw = rows[0] fir = next(r for r in rows if r["name"].startswith("fir")) native = rows[-1] - # Two ratios, because one is clean and the other is not. Decimating a - # capture cannot change its gain or move the sensor, so the first is a - # measurement of the filter alone. The second brings in a whole separate - # run, and carries every difference between the two runs with it. + # Two ratios. Decimating leaves the gain and the sensor position + # untouched, so the first measures the filter alone; the second brings + # in a separate run, carrying every difference between the two. print(f" -> decimating this capture changes its own floor by " - f"{(fir['asd_ppm'] / raw['asd_ppm'] - 1) * 100:+.1f}% " - f"-- no other difference is involved.") + f"{(fir['asd_ppm'] / raw['asd_ppm'] - 1) * 100:+.1f}%, " + f"this capture against itself.") print(f" Against the separately-recorded slow capture it is " f"{fir['asd_ppm'] / native['asd_ppm']:.3f}x, which also carries " f"whatever\n differs between two runs.") - print("\n Compare the ASD ppm column: each path has a different effective " - "noise bandwidth,\n so per-sample sd figures are not on the same " - "footing, and the two cycle counts do\n not share a gain, so " - "absolute nT are not either.") + print("\n Compare the ASD ppm column, which is quoted over a common band " + "and divided by\n each capture's own |B|. Each path carries its own " + "effective noise bandwidth in\n the per-sample sd, and its own gain " + "in the absolute nT.") # The candidate explanations for a gain that moves with the supply, in the -# order they are reported. `both` is included because it is the obvious thing -# to try and because it is instructive when it fails: the rate change *is* the -# supply change seen through the oscillator, so applying the two on top of each -# other counts the same physics twice. +# order reported. `both` stacks the two, and its residual shows what that costs: +# the rate change *is* the supply change seen through the oscillator, so the two +# together count the same physics twice. CORRECTIONS = [("raw", "raw"), ("rate", "x measured fs"), ("supply", "x rail volts"), @@ -533,11 +512,11 @@ CORRECTIONS = [("raw", "raw"), def correction_names(cal): """Row labels, naming what the rate row is actually referenced to. - Uncalibrated it is `1 + rate_error`, which measures the rate against - rm3100's *nominal* model -- and that model misfits by a cycle-count - dependent amount, so the row means something only inside a matched - cycle-count contrast. With a calibration it is the measured oscillator - against a stated reference, which means the same thing everywhere. + Uncalibrated it is `1 + rate_error`, measured against rm3100's *nominal* + model, whose misfit is cycle-count dependent, so the row carries meaning + inside a matched cycle-count contrast, where that misfit cancels. With a + calibration it is the measured oscillator against a stated reference, which + means the same thing everywhere. """ renamed = {"rate": "x oscillator", "both": "x osc x volts"} return [(key, renamed.get(key, name) if cal else name) @@ -549,16 +528,12 @@ def corrected_fields(s, supplies, reference_volts, cal=None): rate the chip's own clock is measured every run, so a gain that tracks the clock is correctable from the capture alone. With a calibration - this is the measured oscillator against a stated reference. Without - one it falls back to `1 + rate_error`, which references rm3100's - nominal model -- and that model's error is cycle-count dependent - (+0.84% at cc=100, -0.99% at cc=400 on one unit, same supply), so - the fallback is only meaningful inside a matched cycle-count - contrast, where it cancels. - supply a gain that is simply ratiometric with the rail needs the rail - measured, but nothing from the capture. - both the two together, which is only right if they are independent - mechanisms rather than one seen twice. + this is the measured oscillator against a stated reference; without + one it falls back to `1 + rate_error`, which carries meaning inside + a matched cycle-count contrast. + supply a gain simply ratiometric with the rail. Needs the rail measured, + and nothing from the capture. + both the two together, which holds if they are independent mechanisms. """ raw = s["field"] if cal is not None: @@ -574,7 +549,7 @@ def corrected_fields(s, supplies, reference_volts, cal=None): def print_contrasts(records, supplies, cal=None): - """Condition contrasts at matched cycle count -- the designed comparison.""" + """Condition contrasts at matched cycle count: the designed comparison.""" conditions = conditions_in_order(records) if len(conditions) < 2: return @@ -608,10 +583,9 @@ def print_contrasts(records, supplies, cal=None): f"{right[key] - left[key]:+9,.1f} nT " f"{(right[key]/left[key] - 1) * 100:+6.2f}%") contrasts[(ca, cb)].append((cc, math.log(right["raw"] / left["raw"]))) - # Rather than only testing whether gain goes as 1/V, fit the power - # it actually goes as. -1 is exactly ratiometric, 0 is no supply - # dependence at all, and the rate gets the same treatment because - # both are the same oscillator seen from different ends. + # Fit the power gain actually goes as. -1 is exactly ratiometric, + # 0 is no supply dependence. The rate gets the same treatment: same + # oscillator, seen from the other end. va, vb = supplies.get(ca), supplies.get(cb) if va and vb and va != vb: dv = math.log(vb / va) @@ -619,17 +593,25 @@ def print_contrasts(records, supplies, cal=None): (cc, math.log(right["raw"] / left["raw"]) / dv, math.log(b["fs"] / a["fs"]) / dv)) + if cal is not None: + print("\n The x oscillator rows above hold for any " + "gain_offset_counts, including the\n " + f"{cal.gain_offset_counts:+.2f} this calibration carries. The " + "shape term depends only on the cycle\n count, so at a matched " + "cycle count both conditions carry the same one and it\n " + "cancels. A contrast spanning cycle counts is what measures the " + "offset.") + for (ca, cb), points in exponents.items(): va, vb = supplies[ca], supplies[cb] print(f"\n how {ca} -> {cb} ({va:g} V -> {vb:g} V) scales with the rail:") for cc, field_power, rate_power in points: print(f" cycle count {cc}: |B| ~ V^{field_power:+.3f} " f"rate ~ V^{rate_power:+.3f}") - print(" V^-1 on |B| would be exactly ratiometric -- gain simply " - "proportional to the\n rail. V^0 would be no supply " - "dependence. Read the residual of the `x rail\n volts` row " - "above against the reproducibility figure below before " - "believing\n either endpoint.") + print(" V^-1 on |B| is exactly ratiometric, gain proportional to " + "the rail; V^0 is\n no supply dependence. Read the residual " + "of the `x rail volts` row above\n against the " + "reproducibility figure below to see what these support.") for (ca, cb), points in contrasts.items(): if len(points) < 2: @@ -639,11 +621,11 @@ def print_contrasts(records, supplies, cal=None): f"mean {np.mean(values) * 100:+.2f}%, " f"spread {(max(values) - min(values)) * 100:.2f}%") print(" A small spread means the condition acts the same at every " - "cycle count, which\n is what a scale factor does and what an " - "environmental step does not.") + "cycle count, which\n is the signature of a scale factor. An " + "environmental step spreads wider.") # A condition measured twice bounds everything this experiment cannot - # separate -- handling, gradient, drift -- and is the only honest error bar. + # separate, handling and gradient and drift, and is the honest error bar. print("\n reproducibility within one condition (the error bar on all of " "the above):") any_repeat = False @@ -659,20 +641,18 @@ def print_contrasts(records, supplies, cal=None): + f" spread {(max(fields)/min(fields) - 1) * 100:.2f}%") if any_repeat: print(" These repeats differ in cycle count as well, so the spread " - "is an upper bound\n on run-to-run reproducibility, not a " - "clean measurement of it.") + "is an upper bound\n on run-to-run reproducibility.") else: - print(" none -- no condition was captured twice, so there is no " - "error bar at all.") + print(" none: no condition was captured twice, so there is no error " + "bar at all.") def print_pairwise(records): if len(records) < 2: return print("\n=== pairwise: |B|, and how far the sensor moved ===") - print("|B| survives a rotation, so a difference is only movement if the " - "sensor also\ntranslated through a gradient. These size that doubt; " - "they do not veto the result.") + print("|B| survives a rotation, so a difference is movement where the " + "sensor also\ntranslated through a gradient. These size that doubt.") for (la, _, a), (lb, _, b) in itertools.combinations(records, 2): ratio = b["mean"] / a["mean"] spread = axis_ratio_spread(a, b, ratio) @@ -696,8 +676,8 @@ def print_pairwise(records): moved = angle > ROTATION_OK_DEG or residual > RATIO_SPREAD_OK if moved: print(" -> the sensor moved between these two, so some of the " - "|B| difference may be\n a different place in the field " - "rather than a different scale factor.") + "|B| difference may be\n a different place in the " + "field.") else: print(" -> the sensor held still, so the |B| difference is a " "scale change.") @@ -707,34 +687,13 @@ def print_pairwise(records): # Figure # -------------------------------------------------------------------------- -def style(ax): - ax.set_facecolor(ch.SURFACE) - ax.grid(True, which="both", color=ch.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(ch.GRID) - ax.tick_params(colors=ch.TEXT_SECONDARY, labelsize=9, length=0) - ax.xaxis.label.set_color(ch.TEXT_SECONDARY) - ax.yaxis.label.set_color(ch.TEXT_SECONDARY) - - -def title(ax, text): - ax.set_title(text, loc="left", color=ch.TEXT_PRIMARY, fontsize=12, - fontweight="bold", pad=8) - - def binned_asd(freqs, asd, bins=140): """Median ASD in log-spaced frequency bins, for a trace that can be read. A raw periodogram at these lengths draws tens of thousands of points into a - few hundred pixels, which renders as a solid band whose apparent thickness - is chi-square scatter rather than anything about the sensor. Binning by the - median keeps the floor exactly where the quoted numbers put it -- they are - medians too -- while a line stays a line. Narrow lines survive because a - log bin up here is many raw bins wide only at high frequency, where these - captures have no lines. + few hundred pixels and renders as a solid band whose thickness is chi-square + scatter. Binning by the median keeps the floor where the quoted numbers put + it, those being medians too, and keeps a line a line. """ edges = np.geomspace(freqs[0], freqs[-1], bins + 1) index = np.digitize(freqs, edges) - 1 @@ -760,28 +719,28 @@ def make_figure(records, band, supplies, path, cal=None): fig, axs = plt.subplots(2, 2, figsize=(13.5, 9), dpi=150) fig.patch.set_facecolor(ch.SURFACE) for ax in axs.flat: - style(ax) + ch.style(ax) a = axs[0, 0] a.axvspan(band[0], band[1], color=ch.GRID, alpha=0.55, zorder=0, linewidth=0) 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"]))) - # Fractional, because the conditions do not share a gain. In absolute - # nT/rtHz a condition reading 7% smaller plots 7% lower while being no - # quieter, and the chart would say the opposite of the truth. + freqs, asd = binned_asd(*ch.welch_asd( + s["axes"]["total"], s["fs"], nperseg=ch.segment_length(s["fs"]))) + # Fractional, because each condition has its own gain: in absolute + # nT/rtHz a condition reading 7% smaller plots 7% lower while being + # equally quiet. asd = asd / s["field"] * 1e6 a.loglog(freqs, asd, color=color[condition_of(label)], linestyle=dash[variant_of(label)], linewidth=1.4, alpha=0.9, label=label) a.set_xlim(*ch.FRACTIONAL_ASD_HZ); a.set_ylim(*ch.FRACTIONAL_ASD_PPM) ch.decade_ticks(a) - title(a, "Spectral density of |B|, as a fraction of |B|") + ch.title(a, "Spectral density of |B|, as a fraction of |B|") a.set_xlabel("frequency (Hz)"); a.set_ylabel("ppm of |B| per √Hz") a.annotate(f"shaded: {band[0]:g}–{band[1]:.0f} Hz, the compared band." - "\nnormalised by each capture's own |B|, so a gain" - "\ndifference is not read as a noise difference", + "\nnormalised by each capture's own |B|, which keeps" + "\na gain difference out of the noise figure", xy=(0.02, 0.05), xycoords="axes fraction", color=ch.TEXT_SECONDARY, fontsize=9) a.legend(frameon=False, fontsize=9, labelcolor=ch.TEXT_SECONDARY, @@ -789,7 +748,7 @@ def make_figure(records, band, supplies, path, cal=None): a = axs[0, 1] plotted = _plot_decimation(a, records, band, color, dash) - title(a, "Fast capture, filtered and decimated") + ch.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.FRACTIONAL_ASD_HZ); a.set_ylim(*ch.FRACTIONAL_ASD_PPM) @@ -810,7 +769,7 @@ def make_figure(records, band, supplies, path, cal=None): 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|") + ch.title(a, "Allan deviation of |B|, as a fraction of |B|") a.set_xlabel("averaging time τ (s)"); a.set_ylabel("σ (ppm of |B|)") a.annotate("slope −½ = white noise; upturn = drift", xy=(0.02, 0.04), xycoords="axes fraction", @@ -840,18 +799,18 @@ def _plot_decimation(ax, records, band, color, dash): for condition, (_, fast), (_, slow), k in decimation_pairs( [(label, s) for label, _, s in records]): try: - v = ch.decimate(dict(fast["cap"].axes())["total"], k, "fir") + v = ch.decimate(fast["axes"]["total"], k, "fir") except ValueError: continue hue = color[condition] freqs, asd = binned_asd(*ch.welch_asd(v, fast["fs"] / k, nperseg=ch.segment_length(fast["fs"] / k))) - # Fractional again: the two cycle counts do not share a gain either. + # Fractional again: the two cycle counts have their own gains too. asd = asd / fast["field"] * 1e6 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"], + slow["axes"]["total"], slow["fs"], nperseg=ch.segment_length(slow["fs"]))) asd = asd / slow["field"] * 1e6 ax.loglog(freqs, asd, color=hue, linestyle="-", linewidth=1.5, @@ -892,7 +851,7 @@ def _plot_corrections(ax, records, supplies, color, marker, cal=None): ax.set_yticks([-r for r in range(len(rows))], [n for _, n in rows]) ax.set_ylim(-len(rows) + 0.5, 0.85) ax.margins(x=0.24) - title(ax, "Total field under each candidate correction") + ch.title(ax, "Total field under each candidate correction") ax.set_xlabel("|B| (nT)") ax.annotate("a correction that is right collapses the spread", xy=(0.02, 0.06), xycoords="axes fraction", @@ -932,7 +891,7 @@ def main(): "ratiometric correction; repeatable") ap.add_argument("--calibration", default=None, metavar="JSON", help="reference the rate correction to a measured " - "calibration instead of rm3100's nominal timing model") + "calibration in place of rm3100's nominal model") ap.add_argument("-o", "--output", default=None, metavar="PNG", help="write a comparison figure here") args = ap.parse_args() @@ -963,15 +922,14 @@ def main(): file=sys.stderr) # One band for every capture, so no run is scored over more bandwidth than - # another. Without this a 299 Hz capture is judged over 3-135 Hz and a - # 73 Hz one over 3-33, and the difference reads as a noise difference. + # another; otherwise the difference in bandwidth reads as a noise difference. slowest = min(comparable_rates([(l, cap) for l, _, cap, _ in loaded])) band = args.band or ch.band_for(slowest) if band[1] > ch.BAND_NYQUIST_FRACTION * slowest / 2: print(f"warning: --band reaches {band[1]:g} Hz, past " f"{ch.BAND_NYQUIST_FRACTION * slowest / 2:.2f} Hz where the " - f"slowest path compared -- {slowest:.2f} Hz, decimation " - f"included -- rolls off", file=sys.stderr) + f"slowest path compared ({slowest:.2f} Hz, decimation " + f"included) rolls off", file=sys.stderr) records = [(label, path, summarise(cap, band)) for label, path, cap, _ in loaded] diff --git a/diagnose-comms.py b/diagnose-comms.py index 56023ed..419a22e 100644 --- a/diagnose-comms.py +++ b/diagnose-comms.py @@ -2,7 +2,7 @@ """Hardware diagnostics for the RM3100 / CH347 link. Walks the signal chain from USB inwards and stops at the first thing that is -broken, so a failure points at a specific wire rather than "it doesn't work": +broken, so a failure points at a specific wire: adapter -> I2C bus -> sensor identity -> registers -> analog section @@ -23,7 +23,7 @@ BUS_KHZ = ch347.DEFAULT_SPEED_KHZ def check(label, ok, detail=""): - print(f" [{'OK ' if ok else 'FAIL'}] {label}" + (f" -- {detail}" if detail else "")) + print(f" [{'OK ' if ok else 'FAIL'}] {label}" + (f": {detail}" if detail else "")) return ok @@ -76,7 +76,7 @@ def main(): print("\n-> Registers are unreliable. Try a lower bus speed or shorter wires.") return 1 - print("\n5. Analog section (BIST -- drives the sensor coils)") + print("\n5. Analog section (BIST, which drives the sensor coils)") sensor.write_reg(rm3100.REG_BIST, [BIST_RUN]) sensor.write_reg(rm3100.REG_POLL, [rm3100.POLL_ALL_AXES]) time.sleep(0.5) @@ -102,8 +102,8 @@ def main(): print("\n6. Live measurement") sensor.configure() # BIST zeroes the measurement registers and leaves DRDY set. With - # DRC1=1 only a results read clears it, so discard one -- otherwise - # wait_for_data() returns instantly on the stale flag and reports zeros. + # DRC1=1 only a results read clears it, so discard one; wait_for_data() + # would otherwise return instantly on the stale flag and report zeros. sensor.read_raw() sensor.set_rate(0x96) sensor.start_cmm() @@ -129,7 +129,7 @@ def main(): finally: sensor.stop_cmm() - print("\nAll checks passed -- the link is fully working.") + print("\nAll checks passed, the link is fully working.") return 0 diff --git a/logger.py b/logger.py index 18b4d38..cc72120 100755 --- a/logger.py +++ b/logger.py @@ -7,15 +7,13 @@ Run ./setup.sh first to install the udev rule and create the virtualenv, then: ./.venv/bin/python logger.py --scan-only The chip samples on its own internal schedule (manual sections 5.7.2, 5.8.2), so -`sample_index` is a grid coordinate, not a count of reads. A late read returns the -*newer* measurement rather than a delayed one, so an unnoticed miss would silently -compress the time axis. Misses are therefore detected and written as explicit -placeholder rows, keeping the index contiguous and the gap visible; the run -continues and reports the total at the end. +`sample_index` is a grid coordinate. A late read returns the *newer* measurement, +so every miss is detected and written as an explicit placeholder row, which keeps +the index contiguous and the gap visible. The run continues and reports the total +at the end. Only irreducible facts are written: the sample count, the host clock, and the raw -counts. Chip time, tesla and magnitude are all reconstructed on load from the -header -- see capture.py. +counts. Chip time, tesla and magnitude are reconstructed on load by capture.py. """ import argparse @@ -34,35 +32,30 @@ import rm3100 CSV_FIELDS = ["sample_index", "system_time_unix", "x_raw", "y_raw", "z_raw", "warning"] -# Mains fundamental. Anything sampled below twice this folds the line onto -# signal irreversibly, which is the one error post-processing cannot undo. +# Mains fundamental. Below twice this the line folds onto signal permanently, +# which is the one error post-processing cannot undo. MAINS_HZ = 60.0 CONSOLE_REFRESH_S = 0.05 FLUSH_INTERVAL_S = 0.5 # How long to let the writer finish its backlog once sampling has stopped. Far -# longer than a healthy writer needs; exceeding it means rows are lost, which is -# reported rather than waited out indefinitely. +# longer than a healthy writer needs, so reaching it means rows are lost and the +# run says so. WRITER_DRAIN_TIMEOUT_S = 30.0 -# The sampling loop releases the GIL inside each USB transfer, then has to take -# it back. Python's default 5 ms switch interval means the writer thread can -# hold it for that long, stalling the sampler and inflating the DRDY bracket -- -# measured at ~10 ms against a 6.67 ms budget at 150 Hz, while the same loop -# without a writer thread peaked at 3.1 ms. Handing off 10x more often costs -# negligible throughput (the writer is not CPU-bound) and keeps the stall well -# inside the bracket budget. +# The sampling loop releases the GIL inside each USB transfer and has to take it +# back, so Python's default 5 ms switch interval becomes the jitter floor: the +# writer can hold the GIL that long, stalling the sampler and inflating the DRDY +# bracket past its budget. Handing off 10x more often costs little throughput, +# the writer being I/O bound. GIL_SWITCH_INTERVAL_S = 0.0005 # Thread priority. The sampling loop runs on the main thread and is the only -# latency-sensitive part; the writer merely formats a few fields per sample. -# -# Lowering the writer needs no privilege and is always done: only the *relative* -# priority matters, and on Linux nice is per-thread, so it does not touch the -# sampler. Raising the sampler needs CAP_SYS_NICE and is opt-in via -# --high-priority, which fails loudly rather than degrading silently -- asking -# for it and not getting it is worth knowing about. +# latency-sensitive part. Lowering the writer needs no privilege and is always +# done: the *relative* priority is what matters, and on Linux nice is per-thread. +# Raising the sampler needs CAP_SYS_NICE, so it is opt-in via --high-priority, +# which fails loudly when the privilege is missing. SAMPLER_NICE = -10 WRITER_NICE = 10 @@ -70,10 +63,8 @@ WRITER_NICE = 10 def restore_ownership(path): """Hand a file created under sudo back to the invoking user. - --high-priority needs root, and anything root writes stays root-owned -- - which earlier in this project produced captures the normal user could not - rewrite. Undo that here so privilege is needed for scheduling and nothing - else leaks from it. + --high-priority needs root and anything root writes stays root-owned, so + the ownership is handed back here. Privilege buys scheduling and nothing else. """ if os.geteuid() != 0: return @@ -91,9 +82,8 @@ def set_thread_nice(value): """Renice the calling thread. Returns the nice actually in effect, or None. Linux threads are tasks, so PRIO_PROCESS with a native thread id applies to - just this thread. Elsewhere this may be a no-op or affect the process, which - is why failure is tolerated rather than fatal -- priority is an optimisation, - not a correctness requirement. + just this thread. Elsewhere it may be a no-op or reach the whole process. + Failure is tolerated, priority being an optimisation. """ try: tid = threading.get_native_id() @@ -102,59 +92,19 @@ def set_thread_nice(value): except (AttributeError, OSError): return None -# Miss detection uses the DRDY bracket -- the span between the last poll showing -# DRDY clear and the poll showing it set -- rather than the read-to-read -# interval. -# -# The bracket is exact where the interval is only a heuristic. If DRDY reads -# clear at t_c then every earlier measurement has already been read, since a -# results read is what clears it. Measurements complete one period apart, so a -# bracket narrower than one period can contain at most one completion, and DRDY -# going high proves it contained at least one. Exactly one, whatever the host -# was doing beforehand. -# -# The read-to-read interval cannot make that claim: measured on this rig it -# reaches 37 ms against a 28.8 ms period -- a 29% overshoot from host stalls -# alone -- while the bracket stays under 10 ms. Thresholding the interval -# therefore flags healthy captures as lossy. -# -# The bracket is compared against the period measured by calibrate_period(), not -# the table value: the latter is 6-9% out on this unit, which is enough to -# miscount grid points inside a bracket and slip the sample index. - -# Bounded so a stalled writer degrades predictably instead of exhausting memory -# on a long run. Far above the depth a healthy writer ever reaches, so filling it -# means the writer is wedged rather than merely behind -- and the run then ends, -# because a row that never reaches the file would leave a hole in sample_index -# and capture.py rejects a non-contiguous index outright. Stopping keeps -# everything already written valid. +# Bounded, so a stalled writer degrades predictably on a long run. Far above the +# depth a healthy writer reaches, so filling it means the writer is wedged. The +# run then ends, which keeps sample_index contiguous over everything written. QUEUE_MAX = 200_000 _SENTINEL = object() -# Rows carry flags in a `warning` column, space separated, empty when fine. -# Keeping them out of the numeric columns means x/y/z stay parseable as numbers -# (blank for a row with no data) and the column generalises to future flags. -# -# MISSED this row is a placeholder -- the measurement was never read, so it -# has no data. The row exists to keep sample_index contiguous, so it -# remains a valid chip-time grid coordinate and the gap is explicit -# rather than silently compressing the timeline. -# AMBIGUOUS the gap ending at this row was of uncertain length, so the index -# may have slipped from here on. -# -# The two are independent. A gap measuring 1.35 periods rounds to one, so no -# placeholder is written -- yet it sits far enough from an integer to distrust, -# and that case carries AMBIGUOUS on the *real* sample terminating the gap, -# which keeps its data. A placeholder inside an uncertain gap carries both. -# -# The flag is a confidence measure, not a claim that 1.35 and 1.65 are equally -# likely: 1.35 probably is one period and 1.65 probably two. Deciding which -# needs neighbouring timestamps and assumptions, and can still be wrong, which -# is exactly why the rows are flagged rather than silently resolved. -WARN_MISSED = "MISSED" -WARN_AMBIGUOUS = "AMBIGUOUS" +# Rows carry flags in a `warning` column, space separated, empty when fine. See +# the capture format in the README; the two are independent, so a row can carry +# either or both. +WARN_MISSED = "MISSED" # placeholder: the measurement was never read +WARN_AMBIGUOUS = "AMBIGUOUS" # gap of uncertain length ends here # How near a half-period the rounding may fall before the count is a coin toss. AMBIGUITY_MARGIN = 0.25 @@ -176,7 +126,7 @@ def parse_args(): p.add_argument("--tmrc", type=lambda s: int(s, 0), default=None, help="continuous-mode rate register, 0x92 (fastest) to 0x9D " "(slowest). Default is the fastest, which lets the cycle " - "count set the rate -- give one only to sample slower " + "count set the rate; give one only to sample slower " "than the cycle count allows") p.add_argument("--duration", type=float, default=0.0, help="seconds to log, or 0 to run until Ctrl-C (default: %(default)s)") @@ -225,9 +175,9 @@ def print_plan(cfg, bus_speed, requested_rate, sampler_nice=None): f"{rm3100.AXIS_OVERHEAD_S * 1e6:.1f} us))") print(f" sensor ceiling {1 / rm3100.sample_period(cfg.cycle_count):.2f} Hz") print(f" TMRC 0x{cfg.tmrc:02x} = " - f"{rm3100.TMRC_RATES[cfg.tmrc]:g} Hz -- {cfg.governed_by} governs") + f"{rm3100.TMRC_RATES[cfg.tmrc]:g} Hz, {cfg.governed_by} governs") print(f" predicted rate {cfg.predicted_hz:.2f} Hz " - f"+/-{rm3100.RATE_TOLERANCE:.0%} (oscillator tolerance; measured below)") + f"+/-{rm3100.RATE_TOLERANCE:.0%} (nominal-model spread; measured below)") print(f" resolution {lsb:.2f} nT/LSB dither {dither:.2f} LSB " "at spec noise") print(f" duty {cfg.duty:.1%} integration / period") @@ -247,40 +197,46 @@ def print_plan(cfg, bus_speed, requested_rate, sampler_nice=None): # wrong and what would fix it, then run what was requested. warnings = [] - # Below the recommended cycle count the LSB grows faster than the sensor's - # own noise, so the quantiser stops being dithered and averaging stalls. - # The manual's hard floor is 30 (section 5.1), but the margin is already - # slim at 50, so anything under that is worth saying out loud. + # What a low cycle count costs is duty; see rm3100.MIN_CYCLE_COUNT. Dither + # is the more commonly cited worry, so it is reported alongside. if cfg.cycle_count < rm3100.RECOMMENDED_MIN_CYCLE_COUNT: + reference_duty = (rm3100.integration_time(rm3100.DEFAULT_CYCLE_COUNT) + / rm3100.sample_period(rm3100.DEFAULT_CYCLE_COUNT)) + penalty = (reference_duty / cfg.duty) ** 0.5 - 1.0 warnings.append( f"cycle count {cfg.cycle_count} is below the recommended " - f"{rm3100.RECOMMENDED_MIN_CYCLE_COUNT}: dither is {dither:.2f} LSB " - f"at spec noise ({lsb:.1f} nT/LSB) against the ~0.2 where averaging " - "stops recovering sub-LSB resolution. characterize.py prints sd/LSB " + f"{rm3100.RECOMMENDED_MIN_CYCLE_COUNT}: the per-axis overhead holds " + f"duty to {cfg.duty:.0%} against {reference_duty:.0%} at cycle " + f"count {rm3100.DEFAULT_CYCLE_COUNT}, which costs {penalty:.0%} of " + f"post-filter noise floor. Dither is {dither:.2f} LSB at spec noise " + f"({lsb:.1f} nT/LSB); characterize.py prints the measured sd/LSB " "per axis, which settles it from the first capture.") - # The measured worst configuration mistake -- TMRC slower than the - # cycle-count ceiling leaves the sensor idle, and idle time buys nothing. - # Idle is the shortfall against the ceiling, which is what TMRC costs; it is - # not 1 - duty, since the fixed per-axis overhead is active time too. + # The worst configuration mistake: TMRC slower than the cycle-count ceiling + # leaves the sensor idle, and idle time buys nothing. Idle is the shortfall + # against the ceiling; the per-axis overhead is active time too. idle = 1.0 - cfg.predicted_hz * rm3100.sample_period(cfg.cycle_count) if cfg.governed_by == "TMRC" and idle > 0.2: + # Against the duty the *same rate* would reach on the cycle count alone, + # which is not 100%: the per-axis overhead is a fixed cost per sample, so + # even a free-running sensor gives some of it back at high rates. + best_duty = 1.0 - rm3100.AXES * rm3100.AXIS_OVERHEAD_S * cfg.predicted_hz warnings.append( f"TMRC 0x{cfg.tmrc:02x} governs and leaves the sensor idle " - f"{idle:.0%} of each period, holding duty to {cfg.duty:.0%}. " - "Measured cost at 23% duty was 1.43x the noise ASD. Raising the " - "cycle count instead reaches the same rate at ~100% duty.") + f"{idle:.0%} of each period, holding duty to {cfg.duty:.0%} against " + f"the {best_duty:.0%} this rate reaches on the cycle count alone. " + f"Post-filter noise goes as 1/sqrt(duty), so that is " + f"{(best_duty / cfg.duty) ** 0.5:.2f}x the noise floor for nothing.") - # A rate that came back different from the one asked for is easy to miss in - # a note, and every figure above is derived from the rate. + # Every figure above is derived from the rate, so a rate that came back + # different from the one asked for is worth saying out loud. if (requested_rate is not None and abs(cfg.predicted_hz / requested_rate - 1) > rm3100.RATE_TOLERANCE): warnings.append( f"{requested_rate:g} Hz was requested but this configuration runs " f"at {cfg.predicted_hz:.2f} Hz.") - # Mains and its harmonics fold onto signal below 2x mains, and no amount of - # post-filtering undoes an alias. + # Below 2x mains the line folds onto signal, and no filter undoes an alias. if cfg.predicted_hz < 2 * MAINS_HZ: alias = abs(MAINS_HZ - round(MAINS_HZ / cfg.predicted_hz) * cfg.predicted_hz) warnings.append( @@ -289,8 +245,6 @@ def print_plan(cfg, bus_speed, requested_rate, sampler_nice=None): "afterwards. Sampling faster and decimating gives the same noise " "floor with the mains line still visible.") - # Table 3-1 stops at 400, so the noise figure printed above is extrapolation - # past that point rather than a specification. if cfg.cycle_count > rm3100.MAX_SPEC_CYCLE_COUNT: warnings.append( f"cycle count {cfg.cycle_count:,} is past the {rm3100.MAX_SPEC_CYCLE_COUNT} " @@ -342,7 +296,7 @@ def start_sensor(sensor, cycle_count, tmrc): """Identify and configure the sensor, printing each step. Returns REVID.""" revid = sensor.revid() if revid == rm3100.EXPECTED_REVID: - print(f"REVID 0x{revid:02x} -- RM3100 confirmed") + print(f"REVID 0x{revid:02x}, RM3100 confirmed") else: print(f"WARNING: REVID 0x{revid:02x}, expected " f"0x{rm3100.EXPECTED_REVID:02x}", file=sys.stderr) @@ -354,7 +308,7 @@ def start_sensor(sensor, cycle_count, tmrc): f"{(cycle_count,) * 3}", file=sys.stderr) sys.exit(1) lsb_nt = rm3100.tesla_per_count(cycle_count) * rm3100.NT_PER_TESLA - print(f"Cycle counts set to {readback} -- 1 LSB = {lsb_nt:.2f} nT") + print(f"Cycle counts set to {readback}, 1 LSB = {lsb_nt:.2f} nT") # Raises if DRC1 did not take; exactly-once sampling depends on it. sensor.configure() @@ -388,11 +342,10 @@ def writer_thread(q, path, meta, stats): break index, system_time, cx, cy, cz, warning = item - # 6 decimal places is microsecond resolution on an epoch value, - # far finer than the ~1 ms bracket uncertainty on each timestamp. - # Placeholder rows carry zeros rather than blanks: the warning - # column already says the row has no data, so keeping x/y/z strictly - # integer makes the file trivial to parse. + # 6 decimals is microsecond resolution on an epoch value, far finer + # than the ~1 ms bracket uncertainty. Placeholders carry zeros rather + # than blanks: the warning column already says the row has no data, + # so x/y/z stay strictly integer and trivial to parse. out.writerow([index, f"{system_time:.6f}", 0 if cx is None else cx, 0 if cy is None else cy, @@ -405,7 +358,7 @@ def writer_thread(q, path, meta, stats): last_flush = now if now - last_print >= CONSOLE_REFRESH_S and cx is not None: # Placeholders carry no field values, so the display holds the - # last real reading rather than blanking or showing nonsense. + # last real reading. last_print = now x, y, z = cx * lsb_nt, cy * lsb_nt, cz * lsb_nt print(f"\r{stats['rows']:,} samples ({stats['missed']} missed)" @@ -418,14 +371,14 @@ def writer_thread(q, path, meta, stats): def i2c_bus_time(bus_speed_khz): """Seconds of I2C traffic the bus carries per sample. - A combined register read of n payload bytes clocks n+3 bytes -- write - address, register, read address, payload -- at 9 bit-times each (8 data plus - ACK), and spends roughly a bit-time on each START and STOP. One sample costs - a STATUS poll (1 payload byte) plus a results read (9). + A combined register read of n payload bytes clocks n+3 bytes (write address, + register, read address, payload) at 9 bit-times each, 8 data plus ACK, and + spends roughly a bit-time on each START and STOP. One sample costs a STATUS + poll of 1 payload byte plus a results read of 9. - This is the irreducible traffic, not bus occupancy: the loop polls - continuously while waiting for DRDY, so raw occupancy approaches 100% and - says nothing useful. What matters is this against the sample period. + This is the irreducible traffic, and the figure worth reading against the + sample period. The loop polls continuously while waiting for DRDY, so raw + bus occupancy approaches 100%. """ bits = 0 for payload in (1, 9): @@ -436,34 +389,29 @@ def i2c_bus_time(bus_speed_khz): def calibrate_period(sensor, dt_nominal, seconds, attempts=3): """Measure the true measurement period from a clean stretch of samples. - Miss detection needs the real period, not the table value: the RC oscillator - runs up to 7% off (manual section 5.2.1) and this unit is ~6% slow, which is - enough to miscount grid points inside a bracket and slip the sample index. + Miss detection needs the real period. The table value is a nominal model, + out by up to ~6% on this unit, which is enough to miscount grid points + inside a bracket and slip the sample index. - Least-squares slope of read time against index. Over one second at 35 Hz - (~35 points, ~0.5 ms jitter) that pins the period to roughly 0.03%, against - the 6-9% error of the nominal value -- and 0.03% is far finer than the - resolution needed to tell k grid points from k+1. + Least-squares slope of read time against index, which over one second pins + the period to roughly 0.03%, far finer than the resolution needed to tell k + grid points from k+1. - Requires a stretch with no miss, since an interval spanning a miss is a - multiple of the period, not one period. Returns None if no clean stretch is - available, which is itself the answer: the rate is not sustainable. + Requires a stretch with no miss, since an interval spanning one is a multiple + of the period. Returns None when no clean stretch is available, which is + itself the answer: the rate is not sustainable. - Also returns the mean host cost of one sample's traffic -- a DRDY poll plus - a results read -- timed here at no extra cost, since those are the very - transactions being performed. + Also returns the mean host cost of one sample's traffic, a DRDY poll plus a + results read, timed here for free from the transactions being performed. """ poll_ready, read_raw = sensor.poll_ready, sensor.read_raw monotonic = time.monotonic - # Bootstrap the period from observation rather than trusting the TMRC table. - # When the cycle count governs instead of TMRC -- which it does whenever the - # cycle-count ceiling falls below the requested rate (manual section 5.2.1) - # -- the table value is not merely 7% out but wrong by a large factor. At - # cc=400 with TMRC 0x92 the table says 1.67 ms and the truth is 13.5 ms, so - # a threshold built on it rejects every interval and calibration can never - # start. A handful of raw intervals settles it: the median is robust to the - # occasional stall, and only needs to be close enough to seed the real fit. + # Bootstrap the period from observation rather than the TMRC table. When the + # cycle count governs, the table value is out by a large factor: at cc=400 + # with TMRC 0x92 it says 1.67 ms against a real 13.5 ms, and a threshold + # built on it rejects every interval. The median of a handful of raw + # intervals is robust to the occasional stall and seeds the real fit. probe = [] clear = monotonic() deadline = clear + max(seconds, 0.25) @@ -496,9 +444,8 @@ def calibrate_period(sensor, dt_nominal, seconds, attempts=3): if ready - clear >= dt_nominal: # conservative: nominal is short clean = False break - # Timed from here, not from `ready`: `ready` is the instant the - # STATUS reply landed on the wire, so measuring against it would - # charge the tail of the poll transaction to the read as well. + # Timed from here. `ready` is the instant the STATUS reply landed + # on the wire, so it carries the tail of the poll transaction. before = monotonic() read_raw() read_cost += monotonic() - before @@ -519,23 +466,27 @@ def calibrate_period(sensor, dt_nominal, seconds, attempts=3): def sample_loop(sensor, q, duration, dt_nominal, stats): """Read every measurement, recording placeholders for any that are lost. - The timestamp recorded per sample is the midpoint of the bracket between the - last poll showing DRDY clear and the poll showing it set: the measurement - completed somewhere in that window, and the midpoint is the best estimate - available without a hardware DRDY line. + Miss detection uses the **DRDY bracket**: the span between the last poll + showing DRDY clear and the poll showing it set. The bracket is exact. A + results read is what clears DRDY, so everything earlier has been read; + measurements complete one period apart, so a bracket narrower than a period + holds at most one completion, and DRDY going high proves it held at least + one. Host stalls push the read-to-read interval well past a period, and this + is unaffected by them. - Ends the run rather than dropping a row if the queue fills: sample_index has - to stay contiguous to remain a chip-time grid coordinate, so a truncated - capture is worth more than a longer one with a hole in it. + The timestamp recorded per sample is the midpoint of that bracket, the best + estimate available without a hardware DRDY line. + + Ends the run if the queue fills, which keeps sample_index contiguous and so + keeps every row written a valid chip-time grid coordinate. """ read_raw = sensor.read_raw poll_ready = sensor.poll_ready monotonic = time.monotonic put = q.put_nowait - # DRDY may already be set from whatever ran before -- BIST in particular - # zeroes the result registers and leaves it high. Discard one reading so - # sample 0 is a genuinely fresh measurement. + # DRDY may already be set from whatever ran before; BIST zeroes the result + # registers and leaves it high. Discard one reading so sample 0 is fresh. if poll_ready()[0]: read_raw() @@ -546,7 +497,7 @@ def sample_loop(sensor, q, duration, dt_nominal, stats): # dt_nominal here is the *calibrated* period from calibrate_period(), not # the table value. Keep refining it from clean intervals so the estimate - # follows the oscillator's thermal drift (~2500 ppm measured over 13 h). + # follows the oscillator's thermal drift. dt_est = dt_nominal previous_mono = previous_wall = None DT_SMOOTHING = 0.02 @@ -555,7 +506,7 @@ def sample_loop(sensor, q, duration, dt_nominal, stats): if monotonic() >= deadline: return # The stamps come back with the STATUS reply itself, taken the instant - # it landed rather than several frames later up here. + # it landed. ready, ready_mono, ready_wall = poll_ready() if not ready: clear_mono, clear_wall = ready_mono, ready_wall @@ -566,34 +517,20 @@ def sample_loop(sensor, q, duration, dt_nominal, stats): if bracket > stats["max_bracket"]: stats["max_bracket"] = bracket - # Counting lost measurements in two stages, because the two questions - # have different best answers. - # - # "Did we lose any?" is settled by the bracket, and rigorously: a - # bracket shorter than one period cannot contain two completions, so - # nothing was lost regardless of timing precision. - # - # "How many?" is better answered by the interval since the previous - # accepted sample. Each accepted read sits on a grid point, so that - # interval is very close to an exact multiple of the period, and - # rounding recovers the multiple. The bracket alone cannot -- it only - # spans back to the last poll that saw DRDY clear, which discards where - # the grid actually is. - # - # Precision: the period is calibrated to ~0.07% and each completion is - # located to about half a poll interval (~0.2 ms, under 1% of a period). - # Over a gap of k periods the total error is well under a tenth of a - # period for any small k, so the rounding is unambiguous -- which is - # the case that matters, since a long gap means the run is unusable - # anyway. + # Two stages, because the two questions have different best answers. + # "Did we lose any?" is settled rigorously by the bracket, above. "How + # many?" is better answered by the interval since the previous accepted + # sample: each accepted read sits on a grid point, so that interval is + # very close to an exact multiple of the period and rounding recovers + # the multiple. The bracket reaches back only to the last poll that saw + # DRDY clear, which loses track of where the grid sits. lost = 0 uncertain = False if bracket >= dt_est and previous_mono is not None: periods = (ready_mono - previous_mono) / dt_est lost = max(0, round(periods) - 1) - # Uncertainty is a property of the gap, not of the placeholders: - # a gap rounding to zero losses can still be a coin toss, and then - # there is no placeholder to carry the flag. + # Uncertainty is a property of the gap: one rounding to zero losses + # can still be a coin toss, with no placeholder to carry the flag. uncertain = abs(periods - round(periods)) > AMBIGUITY_MARGIN if uncertain: stats["ambiguous"] += 1 @@ -603,11 +540,9 @@ def sample_loop(sensor, q, duration, dt_nominal, stats): dt_est += (ready_mono - previous_mono - dt_est) * DT_SMOOTHING if lost: stats["missed"] += lost - # Placeholders sit on the grid -- one period after the previous - # accepted sample, and so on -- because that is where the lost - # measurements actually completed. Spreading them across the - # bracket instead would bunch them at the end of the gap and skew - # the rate fit. + # Placeholders sit on the grid, one period apart from the previous + # accepted sample, which is where the lost measurements completed + # and what keeps the rate fit straight. flags = (f"{WARN_MISSED} {WARN_AMBIGUOUS}" if uncertain else WARN_MISSED) for k in range(lost): @@ -648,7 +583,7 @@ def main(): "Re-run under sudo, using the venv interpreter by absolute " "path:\n" f" sudo {sys.executable} {' '.join(sys.argv)}\n" - "Or drop the flag -- the writer thread already steps aside, " + "Or drop the flag: the writer thread already steps aside, " f"which is the half that needs no privilege.") else: try: @@ -694,9 +629,8 @@ def main(): print() nominal_rate = cfg.predicted_hz - # Establish the real period before recording. Everything downstream -- - # miss detection, placeholder counting, the sample-index grid -- depends - # on it, and the table value is 6-9% out on this unit. + # Miss detection, placeholder counting and the sample-index grid all + # depend on the real period, so measure it first. dt, host_cost = calibrate_period(sensor, 1.0 / nominal_rate, args.calibrate) if dt is None: @@ -708,16 +642,15 @@ def main(): "misses cannot be\ncounted reliably and the sample index would " "not track chip time.\n" "Lower the rate (higher --tmrc) or raise --bus-speed.") - # Compared against the predicted rate, not TMRC's nominal: when the - # cycle count governs, the TMRC figure is not what was aimed for and the - # error against it is meaningless. + # Against the predicted rate, which is what the configuration aimed + # for once the cycle count governs. print(f"Calibrated period {dt * 1e3:.4f} ms = {1 / dt:.4f} Hz " f"({(1 / dt) / cfg.predicted_hz - 1:+.2%} vs predicted)") - # Whether the host is the constraint is a measurement, not a given: it - # costs a near-constant ~0.9 ms per sample regardless of configuration, - # so it dominates at short periods and disappears at long ones. + # Whether the host is the constraint is a measurement. It costs a + # near-constant ~0.9 ms per sample, so it dominates at short periods + # and disappears at long ones. share = host_cost / dt - margin = ("little margin -- expect misses under load" if share > 0.5 + margin = ("little margin, expect misses under load" if share > 0.5 else "modest margin" if share > 0.25 else "ample margin") print(f"Host cost {host_cost * 1e3:.3f} ms/sample " @@ -725,9 +658,7 @@ def main(): meta = { "rm3100_capture": 1, - # The rate this configuration is predicted to produce, whichever of - # the two ceilings governs -- not TMRC's table value, which is not - # what was aimed for when the cycle count wins. + # Whichever of the two ceilings governs. "nominal_rate_hz": nominal_rate, "tmrc_nominal_hz": rm3100.TMRC_RATES[args.tmrc], "tmrc": f"0x{args.tmrc:02x}", @@ -737,9 +668,8 @@ def main(): "i2c_address": f"0x{address:02x}", "bus_speed_khz": args.bus_speed, "revid": f"0x{revid:02x}", - # What the logger actually used for miss detection. Not recoverable - # from the data, since capture.py fits the whole run rather than the - # first second. + # What miss detection actually used. capture.py fits the whole run, + # so this figure lives only here. "calibrated_period_s": repr(dt), } if args.note: @@ -747,7 +677,7 @@ def main(): meta["note"] = " ".join(args.note.split()) path = args.output or datetime.now().strftime("rm3100_%Y%m%d_%H%M%S.csv") - print(f"Logging to {path} -- Ctrl-C to stop\n") + print(f"Logging to {path}, Ctrl-C to stop\n") q = queue.Queue(maxsize=QUEUE_MAX) stats = {"rows": 0, "missed": 0, "ambiguous": 0, "max_bracket": 0.0, @@ -772,54 +702,53 @@ def main(): print(f"WARNING: could not stop CMM: {exc}", file=sys.stderr) print(f"Wrote {stats['rows']} samples to {path}") - # How close the run came to losing a measurement: the useful number for - # judging whether a rate is sustainable before committing to a long run. - # Quoted against the calibrated period, which is the threshold actually - # enforced -- the nominal one is 6-9% out. + # How close the run came to losing a measurement: the number to watch + # before committing to a long run. Quoted against the calibrated period, + # which is the threshold enforced. print(f"Worst DRDY bracket {stats['max_bracket'] * 1e3:.2f} ms of " f"{dt * 1e3:.2f} ms allowed " f"({stats['max_bracket'] / dt * 100:.0f}% of margin used)") - # Every one of these means the capture is usable but compromised in a - # way the file alone does not show, so each is said out loud and each - # makes the exit status non-zero. + # Each of these leaves the capture usable but compromised in a way the + # file alone conceals, so each is said out loud and each makes the exit + # status non-zero. failed = False if stats["truncated"]: failed = True - print("\nWARNING: recording stopped early -- the writer could not " - "keep up and the queue\nfilled. Continuing would have left a " - "hole in sample_index, which capture.py rejects,\nso the run " - "ended instead. Everything written is valid.", file=sys.stderr) + print("\nWARNING: recording stopped early. The writer could not " + "keep up and the\nqueue filled, so the run ended to keep " + "sample_index contiguous.\nEverything written is valid.", + file=sys.stderr) if not drained: failed = True - # The writer is a daemon thread, so it dies with the process and - # anything still queued never reaches the file. Say so: a short tail - # is not visible from the file itself. + # The writer is a daemon thread, so anything still queued dies with + # the process. Say so, since a short tail is invisible in the file. print(f"\nWARNING: the writer did not finish within " - f"{WRITER_DRAIN_TIMEOUT_S:g} s, so the last {q.qsize():,} " - "queued row(s) were\nnot written. The file is valid, just " - "short.", file=sys.stderr) + f"{WRITER_DRAIN_TIMEOUT_S:g} s, so the last\n" + f"{q.qsize():,} queued row(s) were dropped. The file is valid " + "up to where it stops.", file=sys.stderr) if stats["ambiguous"]: failed = True print(f"\nWARNING: {stats['ambiguous']:,} gap(s) could not be " - "counted confidently -- the interval fell near a half-period, " - "so\nthe number of lost measurements is a guess and " - f"sample_index may have slipped.\nThose rows carry " - f"{WARN_AMBIGUOUS} in the warning column. Do not trust this " - "capture's\ntime axis for spectral work.", file=sys.stderr) + "counted confidently. The interval fell near a\nhalf-period, " + "so the number of lost measurements is a guess and " + "sample_index\nmay have slipped from there on. Those rows " + f"carry {WARN_AMBIGUOUS} in the\nwarning column. Re-record " + "before using this capture's time axis.", + file=sys.stderr) if stats["missed"]: failed = True pct = stats["missed"] / max(stats["rows"], 1) * 100 print(f"\nWARNING: {stats['missed']:,} measurement(s) were lost " f"({pct:.3f}% of rows), flagged {WARN_MISSED}.\n" - "The sample index still tracks chip time -- the gaps are " - "explicit, not compressed --\nbut those rows carry no field " - "data. Lower the rate (higher --tmrc) or raise --bus-speed\n" - "to remove them.", file=sys.stderr) + "The gaps are explicit, so the sample index still tracks chip " + "time, and those\nrows carry no field data. Lower the rate " + "(higher --tmrc) or raise --bus-speed\nto remove them.", + file=sys.stderr) return 1 if failed else 0 return 0 diff --git a/plot.py b/plot.py index 242e7ee..d7ad160 100644 --- a/plot.py +++ b/plot.py @@ -4,10 +4,9 @@ ./.venv/bin/python plot.py capture_60s.csv ./.venv/bin/python plot.py capture_60s.csv -o out.png -Small multiples rather than one shared axis: the three axes sit at very -different DC offsets, so a single scale would flatten the variation that -matters. Each panel therefore has its own y-scale -- read the panels -independently, and note the per-panel mean/sd annotation for context. +Small multiples, since the three axes sit at very different DC offsets and a +shared scale would flatten the variation that matters. Each panel has its own +y-scale, with a mean/sd annotation for context. """ import argparse @@ -20,24 +19,15 @@ import numpy as np import calibrate import capture +import characterize +from characterize import GRID, SURFACE, TEXT_PRIMARY, TEXT_SECONDARY -# Light-mode design tokens. -SURFACE = "#fcfcfb" -TEXT_PRIMARY = "#0b0b0b" -TEXT_SECONDARY = "#52514e" -GRID = "#e3e2df" - -# Categorical slots 1-3 for the three peer axes. Validated all-pairs in light -# mode (worst CVD dE 9.2, normal-vision 24.0). The norm is a derived quantity, -# not 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 for small -# multiples. -SERIES = [ - ("x", "X axis", "#2a78d6"), - ("y", "Y axis", "#eb6834"), - ("z", "Z axis", "#1baf7a"), - ("total", "Norm |B| = sqrt(X^2 + Y^2 + Z^2)", TEXT_PRIMARY), -] +# The three peer axes take characterize's categorical slots, so a capture keeps +# the same colours in every figure. The norm is derived from them, so it takes +# neutral ink. +SERIES = [(key, f"{name} axis", color) + for key, name, color in characterize.AXES] +SERIES.append(("total", "Norm |B| = sqrt(X^2 + Y^2 + Z^2)", TEXT_PRIMARY)) def rolling_mean(v, window): @@ -45,7 +35,7 @@ def rolling_mean(v, window): A plain convolution tapers toward zero at the edges. Dividing by the number of samples that actually contributed gives a true partial-window mean - instead, so the ends carry no artefact. + instead, so the ends carry no artifact. """ if window < 2: return v @@ -66,15 +56,13 @@ def main(): args = ap.parse_args() try: - # Either format: a raw capture, or one already converted to nT by - # calibrate.py. Which it is, is a property of the file. + # Either format: a raw capture, or one already converted to nT. cap = calibrate.load_any(args.csv) except (OSError, capture.CaptureError, calibrate.CalibrationError) as exc: sys.exit(str(exc)) - # Drift-corrected time: the nominal grid would be ~6% off real seconds. - t = cap.elapsed + t = cap.elapsed # drift-corrected; the nominal grid is ~6% off series = dict(cap.axes()) duration = cap.duration rate = cap.true_rate_hz @@ -87,8 +75,7 @@ def main(): v = series[key] ax.set_facecolor(SURFACE) - # Raw trace kept thin and translucent: at ~376 Hz there are far more - # samples than pixels, so a full-weight line would read as a solid band. + # Thin and translucent, there being far more samples than pixels. ax.plot(t, v, color=color, linewidth=0.4, alpha=0.30, solid_capstyle="round", rasterized=True) if window > 1: @@ -96,8 +83,7 @@ def main(): solid_capstyle="round") ax.set_ylabel("nT", color=TEXT_SECONDARY, fontsize=10) - # Direct label instead of a legend: one series per panel, so the title - # names it. This is also the relief the palette's contrast WARN requires. + # Direct label: one series per panel. ax.set_title(label, color=TEXT_PRIMARY, fontsize=12, loc="left", pad=8, fontweight="bold") ax.annotate(f"mean {v.mean():,.0f} nT sd {v.std():.0f} nT " @@ -123,8 +109,8 @@ def main(): if window > 1 else "") fig.suptitle("RM3100 magnetometer capture", color=TEXT_PRIMARY, fontsize=15, fontweight="bold", x=0.5, y=0.985) - # A calibrated file has no nominal rate to be wrong against -- it has a - # gain factor instead, which is the thing worth stating about it. + # A calibrated file carries a gain factor, which is what is worth stating + # about it; a raw one carries a nominal rate to compare against. provenance = (cap.describe() if hasattr(cap, "describe") else f"nominal {cap.nominal_rate_hz:g} Hz, " f"{cap.rate_error * 100:+.1f}%") diff --git a/rm3100.py b/rm3100.py index e4b0432..52b8cd7 100644 --- a/rm3100.py +++ b/rm3100.py @@ -46,22 +46,21 @@ STATUS_DRDY = 0x80 EXPECTED_REVID = 0x22 -# 100 cycle count gives a good balance of sensor duty cycle, spectrum, and post-filter noise floor. -# This should be raised if data rates are higher than desired or bus speed / host speed can't support. -# Lower than 50 is not recommended under any circumstances. See the README for more info. +# Balances duty cycle, spectrum and post-filter noise floor; see the README for +# the trade against 50 and 200. DEFAULT_CYCLE_COUNT = 100 -# Not a hard floor -- the manual's is MIN_CYCLE_COUNT below -- but the point -# under which the dither margin is too thin to defend, so logger.py warns. +# The point under which the duty cycle starts costing real noise floor, so +# logger.py warns. MIN_CYCLE_COUNT below is the hard floor. RECOMMENDED_MIN_CYCLE_COUNT = 50 # Table 3-1 tabulates gain and noise up to a cycle count of 400. Past that both # _gain_lsb_per_tesla and expected_noise_nt are extrapolating their fits. MAX_SPEC_CYCLE_COUNT = 400 -# The fastest rate register. Used by default so the cycle count, not TMRC, sets -# the rate -- TMRC offers only factor-of-two steps and has no effect at all once -# the cycle count governs. +# The fastest rate register, used by default so the cycle count sets the rate. +# TMRC offers factor-of-two steps, and once the cycle count governs it has no +# effect at all. TMRC_FASTEST = 0x92 @@ -86,8 +85,7 @@ def gain_lsb_per_tesla(cycle_count): def tesla_per_count(cycle_count): """Calibration constant: multiply a raw count by this to get tesla. - Stored in capture headers, so a reader needs no knowledge of the gain fit -- - and it is a multiply at the point of use rather than a divide. + Stored in capture headers, so a reader needs no knowledge of the gain fit. """ return 1.0 / gain_lsb_per_tesla(cycle_count) @@ -96,41 +94,27 @@ def tesla_per_count(cycle_count): # # per-axis time = cycle_count / COUNTS_PER_SECOND + AXIS_OVERHEAD_S # -# COUNTS_PER_SECOND comes from the specification, not from a fit: Table 3-1 -# gives a 180 kHz circuit oscillation frequency, and section 4.1 measures each -# cycle count in *both* bias directions, so one count costs two oscillations. -# Rates measured here across cycle counts 229..29769 agree with that figure to -# within 0.8%, which confirms the spec value rather than improving on it -- the -# residual is this particular unit's oscillator sitting inside ordinary -# component tolerance, and another part would sit somewhere else. -# -# The overhead has no specified value and must be measured; 68.7 us reproduces -# the observed rates to better than 1% when the divisor is held at spec. It is -# why the naive "rate x cycle_count" constant is not constant, drifting from -# 84,429 at cycle count 100 to 89,191 at 1200. -# # **This pair is a nominal model, not a description of any particular part.** -# Solving it freely against one unit's two cycle counts gave 88,546 counts/s -# with a 40.6 us overhead at 3.006 V, and 92,889 with 38.1 us at 3.353 V -- so -# both terms are per-unit *and* move with the supply, because the oscillator -# does. Holding the divisor at spec and solving for the overhead alone cannot -# absorb it: that gives 58.9 us at cycle count 100 and 113.6 us at 400 on the -# same data, fitting neither. +# Both terms are per-unit and move with the supply, because the oscillator does, +# and the resulting error is cycle-count dependent, which makes +# `Capture.rate_error` a configuration diagnostic. calibrate.py holds a measured +# per-unit pair; `sweep.py --calibration-out` produces one. # -# The consequence to know about is that the error is *cycle-count dependent*. -# On that unit the model was +0.84% fast at cycle count 100 and -0.99% slow at -# 400 -- same supply, opposite signs -- so `Capture.rate_error` measured against -# it is not a gain indicator. calibrate.py holds the per-unit pair, and -# `sweep.py --calibration-out` measures one. +# COUNTS_PER_SECOND is specified rather than fitted: Table 3-1 gives a 180 kHz +# circuit oscillation and section 4.1 measures each cycle count in *both* bias +# directions, so one count costs two oscillations. The overhead has no specified +# value; 68.7 us is what Table 3-1's own maximum single-axis rates imply. It is +# why the naive "rate x cycle_count" constant is not constant. COUNTS_PER_SECOND = 90000.0 # 180 kHz (Table 3-1) / 2 bias directions AXIS_OVERHEAD_S = 68.7e-6 # nominal; measured parts run 38-41 us AXES = 3 -# Rates predicted from the model are good to roughly this much on a given unit, -# and no better across units, since the manual quotes no tolerance on the -# oscillator. Anything needing the real number measures it: logger.py calibrates -# the period against the host clock before recording. -RATE_TOLERANCE = 0.02 +# How far a rate predicted from the nominal model can sit from the rate actually +# achieved: the *model's* error envelope. The ~7% of section 5.2.1 is a spread +# on the TMRC update rates of Table 5-4, which is a different quantity from the +# cycle-count ceiling that governs here. Anything needing the real number +# measures it: logger.py calibrates the period against the host clock. +RATE_TOLERANCE = 0.06 # Table 3-1 quotes 30/20/15 nT at cycle counts 50/100/200, which fits # K/sqrt(cycle_count). Extrapolation past ~400 is unverified: the manual calls @@ -157,20 +141,18 @@ def cycle_count_for_rate(rate_hz, axes=AXES): """Cycle count that makes the sensor free-run at about rate_hz. Inverse of sample_period(). Use with a TMRC faster than the target so the - cycle count governs -- then this is a continuous rate knob, where TMRC only - offers factor-of-two steps, and the duty cycle is ~100% by construction. - - The achieved rate will sit within about RATE_TOLERANCE of the target, set by - oscillator tolerance rather than by this calculation. Measure it if it - matters. + cycle count governs. It is then a continuous rate knob, where TMRC offers + factor-of-two steps, and the duty cycle is ~100% by construction. + The achieved rate lands within about RATE_TOLERANCE of the target. """ cc = round((1.0 / (rate_hz * axes) - AXIS_OVERHEAD_S) * COUNTS_PER_SECOND) return max(1, min(0xFFFF, cc)) # Section 5.1: "quantization issues generally dictate working above a cycle -# count value of ~30". Below that the LSB grows faster than the sensor's own -# noise, so the quantiser stops being dithered and averaging stalls. +# count value of ~30". Kept as the manual's floor, though what actually gives +# way below 50 is the duty cycle: the per-axis overhead is a fixed cost against +# a shrinking integration time, and post-filter noise goes as 1/sqrt(duty). MIN_CYCLE_COUNT = 30 MAX_CYCLE_COUNT = 0xFFFF MIN_RATE_BY_CYCLE_COUNT = 1.0 / (AXES * (MAX_CYCLE_COUNT / COUNTS_PER_SECOND @@ -186,10 +168,8 @@ def plan(rate_hz=None, cycle_count=None, tmrc=None): cycle count sets how long a measurement takes, TMRC sets how often one is started. Leaving TMRC faster than the cycle-count ceiling means the sensor free-runs at ~100% duty; setting it slower makes the sensor idle, which - costs sensitivity for nothing unless low power is the goal. - - So the default in every branch is TMRC_FASTEST, and the cycle count is the - rate knob -- it is continuous where TMRC offers only factor-of-two steps. + costs sensitivity for nothing unless low power is the goal. So the default + in every branch is TMRC_FASTEST and the cycle count is the rate knob. """ notes = [] if cycle_count is None: @@ -269,10 +249,9 @@ class RM3100: Returns (data, mono, wall): the bus stamps the clocks when the reply lands, and the stamp travels with its own data. - Prefers a single combined transaction (repeated START) where the bus - offers one, since the host round trip is what limits the sample rate. - Falls back to the manual's STOP-then-START form (sections 4.5.2, 5.8.4) - for a bus that cannot do it. + Prefers a combined transaction (repeated START) where the bus offers + one, since the host round trip is what limits the sample rate. Falls + back to the manual's STOP-then-START form (sections 4.5.2, 5.8.4). """ combined = getattr(self.bus, "write_read", None) if combined is not None: @@ -311,9 +290,8 @@ class RM3100: def configure(self): """Put DRDY into a state where polling STATUS actually works. - Verified rather than assumed: exactly-once sampling relies on DRC1=1 - clearing DRDY when the measurement registers are read, which is what - makes a second read of the same measurement impossible. + The readback is verified here, since exactly-once sampling relies on + DRC1=1 clearing DRDY only when the measurement registers are read. """ self.write_reg(REG_HSHAKE, [HSHAKE_DRDY_ON_READ_ONLY]) # Bits 4-6 are read-only NACK status, so compare only the writable ones. @@ -355,10 +333,10 @@ class RM3100: time.sleep(interval) def read_raw(self): - """Return ((x, y, z) counts, mono, wall) -- the fast path. + """Return ((x, y, z) counts, mono, wall): the fast path. - Callers logging at high rates should use this and defer the microtesla - conversion, so the sampling loop does I2C and nothing else. + Callers logging at high rates use this and defer the unit conversion, so + the sampling loop does I2C and nothing else. """ data, mono, wall = self.read_reg(REG_MX, 9) return decode_measurements(data), mono, wall diff --git a/sweep.py b/sweep.py index b4afa63..f86104f 100644 --- a/sweep.py +++ b/sweep.py @@ -5,11 +5,10 @@ ./.venv/bin/python sweep.py --from 8 --to 128 ./.venv/bin/python sweep.py --rates 10,25,50 -Cycle count is the rate knob, not TMRC. TMRC only offers factor-of-two steps -and, once the cycle count governs, has no effect at all -- measured 73.85 Hz at -TMRC 0x92 against 73.86 Hz at 0x94 for the same cycle count. So each point runs -with TMRC set fast and the cycle count chosen from rm3100.cycle_count_for_rate(), -which makes the rate continuous and the duty cycle ~100% by construction. +Cycle count is the rate knob, not TMRC: TMRC offers only factor-of-two steps and +has no effect at all once the cycle count governs. So each point runs with TMRC +set fast and the cycle count from rm3100.cycle_count_for_rate(), which makes the +rate continuous and the duty cycle ~100% by construction. Each row is measured, not predicted: the logger calibrates the true period against the host clock before recording, and that figure is what is reported. @@ -33,7 +32,7 @@ MISSED = re.compile(r"WARNING: ([\d,]+) measurement\(s\) were lost") AMBIGUOUS = re.compile(r"WARNING: ([\d,]+) gap\(s\) could not be counted") # Each sweep point is a short capture, so the oscillator barely drifts within -# one -- but the points are minutes apart and the chip warms through the sweep. +# one, and the points are minutes apart while the chip warms through the sweep. # This is the fractional period uncertainty that implies, and it is what the # timing fit's error bars are built on. PERIOD_SD = 500e-6 @@ -65,10 +64,8 @@ def measure(rate, bus_speed, tmrc, output, record_seconds): out = proc.stdout + proc.stderr found = CALIBRATED.search(out) if found: - # Both warnings matter and they are independent. Reporting only losses - # would hide the worse case: an ambiguous gap means the sample index - # itself may have slipped, so the timeline is suspect even where no - # measurement was lost. + # Both warnings matter and are independent: an ambiguous gap means the + # sample index itself may have slipped, even where nothing was lost. notes = [] for pattern, label in ((MISSED, "lost"), (AMBIGUOUS, "ambiguous")): hit = pattern.search(out) @@ -83,11 +80,10 @@ def measure(rate, bus_speed, tmrc, output, record_seconds): def write_calibration(measured, path, note): """Fit the timing model across the sweep and save it as a calibration. - The sweep already measures what a calibration is made of -- a rate at each - of several cycle counts -- so fitting it here is the difference between a - per-unit reference that was measured and one that was assumed. The reference - oscillator is set to the fitted value, which makes this sweep the definition - of "uncorrected" and every later capture correct relative to it. + The sweep already measures what a calibration is made of: a rate at each of + several cycle counts. The reference oscillator is set to the fitted value, + which makes this sweep the definition of "uncorrected" and every later + capture correct relative to it. """ if len({cc for cc, _ in measured}) < 2: print(f"\nnot writing {path}: the sweep produced " @@ -126,8 +122,8 @@ def write_calibration(measured, path, note): print(f" fit residual {fit.residual_ppm:.0f} ppm over {fit.points} points") print(" gain_exponent is 1.0, which the oscillator-clocked measurement " "interval predicts,\n and gain_offset_counts is 0.0. Neither is " - "measured by a rate sweep -- both need\n a field that does not move " - "between cycle counts. See NOTES.md.") + "measured by a rate sweep. Both need\n a field that holds still " + "between cycle counts. See Noise_Floor_Testing/NOTES.md.") return path @@ -183,9 +179,8 @@ def main(): noise = rm3100.expected_noise_nt(cc) # Fraction of the period spent integrating: what actually reduces noise. duty = rm3100.integration_time(cc) * hz if hz else None - # A note always trails the row rather than sitting in a cell: "host - # cannot sustain" is wider than any sensible column and would shove the - # rest of the line out of alignment. + # A note trails the row: "host cannot sustain" is wider than any + # sensible column. print(row([ f"{rate:g} Hz", f"{hz:.2f} Hz" if hz else "-", @@ -201,12 +196,13 @@ def main(): "measurement: 208/sqrt(cc) nT, fitted to its 30/20/15 nT at 50/100/200.\n" "Past cycle count ~400 the manual gives no data, so those are\n" "extrapolation. Measure the real floor with characterize.py.\n" - "duty is integration time against the period -- what reduces noise.\n" + "duty is integration time against the period, which is what reduces\n" + "noise.\n" "It falls at high rates as the fixed per-axis overhead grows relative\n" "to the integration, and would fall further if TMRC rather than the\n" "cycle count governed, leaving the sensor idle between measurements.\n" - "bus use is irreducible I2C traffic against the period -- not occupancy,\n" - "which approaches 100% because the loop polls continuously for DRDY.") + "bus use is irreducible I2C traffic against the period. Raw occupancy\n" + "approaches 100%, the loop polling continuously for DRDY.") if args.calibration_out: write_calibration(measured, args.calibration_out, args.note) diff --git a/tests/test_analysis.py b/tests/test_analysis.py index efb339c..752da98 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -5,6 +5,7 @@ so a normalisation slip -- the easy mistake in Welch and Allan code, and an invisible one on real data -- shows up as a factor rather than a wobble. """ +import math from pathlib import Path import numpy as np @@ -388,3 +389,99 @@ def test_make_time_base_declines_a_constant_field(tmp_path): cap = _drifted(tmp_path / "c.csv", noise_nt=0.0) with pytest.raises(ValueError, match="constant|no lines"): characterize.make_time_base(cap, str(tmp_path / "c_timebase.png")) + + +# -------------------------------------------------------------------------- +# Measuring the oscillator's drift +# -------------------------------------------------------------------------- + +def test_window_rates_recovers_a_planted_ramp(tmp_path): + """The rate falls at the planted rate, in fractional terms per second.""" + ppm_per_second = 8.0 + cap = _drifted(tmp_path / "d.csv", ppm_per_second=ppm_per_second) + centres, rates = characterize.window_rates(cap, 10.0) + assert len(rates) == 10 + # A period ramping up as (1 + r*t) is a rate falling as (1 - r*t). + slope = np.polyfit(centres, rates / rates.mean(), 1)[0] + assert slope == pytest.approx(-ppm_per_second * 1e-6, rel=0.05) + + +def test_window_rates_is_flat_without_drift(tmp_path): + cap = _drifted(tmp_path / "f.csv") + _, rates = characterize.window_rates(cap, 10.0) + assert np.ptp(rates) / rates.mean() < 1e-9 + + +def test_window_rates_centres_span_the_capture(tmp_path): + cap = _drifted(tmp_path / "c.csv") + centres, rates = characterize.window_rates(cap, 10.0) + assert len(centres) == len(rates) + assert centres[0] == pytest.approx(5.0) + assert centres[-1] == pytest.approx(95.0) + + +@pytest.mark.parametrize("window", [0.0, -1.0, 60.0, 1e6]) +def test_window_rates_declines_impossible_windows(tmp_path, window): + """Non-positive, or too long to fit two windows in the capture.""" + cap = _drifted(tmp_path / "s.csv") + centres, rates = characterize.window_rates(cap, window) + assert len(centres) == 0 and len(rates) == 0 + + +def test_rate_stability_recovers_a_planted_ramp(tmp_path): + """Peak to peak across eight windows, against the ramp that produced it.""" + cap = _drifted(tmp_path / "r.csv", ppm_per_second=8.0) + # Eight windows of 12.5 s: centres 6.25 s and 93.75 s apart, so the spread + # is the ramp over 87.5 s, not over the full 100 s. + expected = 8.0e-6 * cap.duration * (1 - 1 / 8) + assert characterize.rate_stability(cap) == pytest.approx(expected, rel=0.05) + + +def test_rate_stability_is_zero_on_a_uniform_grid(tmp_path): + assert characterize.rate_stability(_drifted(tmp_path / "u.csv")) < 1e-9 + + +def test_rate_stability_declines_a_capture_too_short_to_window(tmp_path): + cap = _drifted(tmp_path / "t.csv", rows=200) + assert characterize.rate_stability(cap) == 0.0 + + +def test_rate_allan_of_a_linear_drift_rises_with_tau(tmp_path): + """sigma_y(tau) = |D| * tau / sqrt(2) for a deterministic frequency ramp. + + This is the shape the real captures show above a few seconds, and it is what + separates a warming oscillator from a random walk. + """ + ppm_per_second = 8.0 + cap = _drifted(tmp_path / "a.csv", ppm_per_second=ppm_per_second) + taus, devs = characterize.rate_allan(cap, [4.0, 8.0, 16.0]) + assert list(taus) == [4.0, 8.0, 16.0] + for tau, dev in zip(taus, devs): + assert dev == pytest.approx(ppm_per_second * 1e-6 * tau / math.sqrt(2), + rel=0.05) + + +def test_rate_allan_falls_with_tau_when_only_jitter_is_present(tmp_path): + """Independent read jitter averages down; drift does not. Opposite slopes.""" + cap = _drifted(tmp_path / "j.csv", jitter_s=2e-3) + taus, devs = characterize.rate_allan(cap, [4.0, 16.0]) + assert devs[0] > devs[1] * 2 + + +def test_rate_allan_is_negligible_on_a_uniform_grid(tmp_path): + cap = _drifted(tmp_path / "z.csv") + _, devs = characterize.rate_allan(cap, [4.0, 8.0]) + assert np.all(devs < 1e-9) + + +def test_rate_allan_drops_taus_that_yield_too_few_windows(tmp_path): + """Two windows give one difference, which is not an estimate of anything.""" + cap = _drifted(tmp_path / "d2.csv") + taus, devs = characterize.rate_allan(cap, [10.0, 40.0, 200.0]) + assert list(taus) == [10.0] # 40 s gives two windows, 200 s none + assert len(devs) == 1 + + +def test_rate_allan_accepts_an_empty_request(tmp_path): + taus, devs = characterize.rate_allan(_drifted(tmp_path / "e.csv"), []) + assert len(taus) == 0 and len(devs) == 0 diff --git a/tests/test_calibrate.py b/tests/test_calibrate.py index 4febbe5..8abb219 100644 --- a/tests/test_calibrate.py +++ b/tests/test_calibrate.py @@ -20,16 +20,47 @@ from synthetic import drifting_times, field_counts, write_capture FIELD = (-40000.0, 16000.0, 12000.0) +REFERENCE_HZ = 88546.0 +REFERENCE_OVERHEAD_S = 40.61e-6 + + def planted_period(cycle_count, counts_per_second, overhead): + """Period for a capture whose per-axis overhead is a fixed *duration*. + + Within one condition this is indistinguishable from the clocked form below, + which is why `fit_timing` -- which only ever sees one condition at a time -- + is tested with it. + """ return rm3100.AXES * (cycle_count / counts_per_second + overhead) -def make(path, cycle_count=100, counts_per_second=88546.0, overhead=40.61e-6, +def planted_period_clocked(cycle_count, counts_per_second, + overhead=REFERENCE_OVERHEAD_S, + reference_hz=REFERENCE_HZ): + """Period for a capture whose overhead is a fixed *count* of cycles. + + The two forms only diverge once a capture runs at a frequency other than the + one the overhead was measured at, and that is exactly the case a calibration + is applied to. This is the model `calibrate.oscillator_hz` inverts: the + overhead's duration goes as 1/f, so a capture at 92,889 counts/s carries + 38.71 us of what was 40.61 us at 88,546. + """ + return rm3100.AXES * (cycle_count + overhead * reference_hz) / counts_per_second + + +def make(path, cycle_count=100, counts_per_second=REFERENCE_HZ, + overhead=REFERENCE_OVERHEAD_S, clocked=False, rows=2000, mean=FIELD, noise_nt=5.0, seed=0): + """A synthetic capture with a planted timing model. + + `clocked` chooses which hypothesis about the overhead the capture obeys; it + makes no difference unless `counts_per_second` differs from REFERENCE_HZ. + """ + period = (planted_period_clocked if clocked else planted_period) counts = field_counts(rows, mean, cycle_count, noise_nt, seed) return capture.load(write_capture( path, rows=rows, cycle_count=cycle_count, counts=counts, - dt=planted_period(cycle_count, counts_per_second, overhead))) + dt=period(cycle_count, counts_per_second, overhead))) def written(tmp_path, **fields): @@ -147,13 +178,56 @@ def test_the_nominal_overhead_costs_far_more_and_varies_with_cycle_count(tmp_pat def test_oscillator_refuses_a_calibration_that_cannot_describe_the_capture(tmp_path): - """An overhead longer than the whole per-axis time is not merely wrong.""" + """An overhead that rivals the cycle count is not merely wrong.""" cap = make(tmp_path / "c.csv", cycle_count=100) absurd = calibrate.calibration(88546.0, axis_overhead_s=1.0) with pytest.raises(calibrate.CalibrationError, match="does not describe"): calibrate.oscillator_hz(cap, absurd) +@pytest.mark.parametrize("cycle_count", [50, 100, 200, 400]) +def test_the_overhead_is_carried_as_counts_not_microseconds(tmp_path, + cycle_count): + """A capture at another frequency is inverted with its own overhead. + + The overhead is a fixed count of oscillator cycles, so a condition running + 4.9% faster spends 4.9% less time on it. Subtracting the reference + condition's *duration* instead would leave a residual, and one that grows as + the overhead's share of the per-axis time grows -- which is why the low + cycle counts are parametrised here. + """ + cap = make(tmp_path / f"{cycle_count}.csv", cycle_count=cycle_count, + counts_per_second=92889.0, clocked=True) + cal = calibrate.calibration(REFERENCE_HZ, REFERENCE_OVERHEAD_S) + assert calibrate.oscillator_hz(cap, cal).value == pytest.approx( + 92889.0, rel=1e-6) + + +def test_a_second_condition_is_recovered_as_its_own_timing_fit_would(tmp_path): + """The single-capture inversion agrees with the two-point fit it replaces. + + `fit_timing` needs two cycle counts and recovers both terms for a condition + outright; `oscillator_hz` claims to get the frequency from one capture and a + calibration taken elsewhere. Those two routes have to land in the same + place, or a calibration cannot be carried across supplies at all -- which is + the entire premise of the file. + """ + periods = {cc: planted_period_clocked(cc, 92889.0) for cc in (100, 400)} + fit = calibrate.fit_timing(list(periods), list(periods.values()), + period_sd=1e-6) + assert fit.counts_per_second == pytest.approx(92889.0, rel=1e-9) + # The overhead scaled by exactly the frequency ratio, which is the claim. + assert fit.axis_overhead_s == pytest.approx( + REFERENCE_OVERHEAD_S * REFERENCE_HZ / 92889.0, rel=1e-9) + + cal = calibrate.calibration(REFERENCE_HZ, REFERENCE_OVERHEAD_S) + for cycle_count in periods: + cap = make(tmp_path / f"{cycle_count}.csv", cycle_count=cycle_count, + counts_per_second=92889.0, clocked=True) + assert calibrate.oscillator_hz(cap, cal).value == pytest.approx( + fit.counts_per_second, rel=1e-6) + + # -------------------------------------------------------------------------- # gain_factor # -------------------------------------------------------------------------- @@ -167,17 +241,20 @@ def test_a_capture_at_its_own_reference_is_not_corrected(tmp_path): def test_gain_factor_is_the_oscillator_ratio_at_exponent_one(tmp_path): - overhead = 40.61e-6 - cap = make(tmp_path / "c.csv", counts_per_second=92889.0, overhead=overhead) - cal = calibrate.calibration(88546.0, overhead, reference_oscillator_hz=88546.0) + overhead = REFERENCE_OVERHEAD_S + cap = make(tmp_path / "c.csv", counts_per_second=92889.0, overhead=overhead, + clocked=True) + cal = calibrate.calibration(REFERENCE_HZ, overhead, + reference_oscillator_hz=REFERENCE_HZ) assert calibrate.gain_factor(cap, cal).value == pytest.approx( - 92889.0 / 88546.0, rel=1e-4) + 92889.0 / REFERENCE_HZ, rel=1e-4) def test_a_zero_exponent_disables_the_oscillator_term(tmp_path): - overhead = 40.61e-6 - cap = make(tmp_path / "c.csv", counts_per_second=92889.0, overhead=overhead) - cal = calibrate.calibration(88546.0, overhead, gain_exponent=0.0) + overhead = REFERENCE_OVERHEAD_S + cap = make(tmp_path / "c.csv", counts_per_second=92889.0, overhead=overhead, + clocked=True) + cal = calibrate.calibration(REFERENCE_HZ, overhead, gain_exponent=0.0) assert calibrate.gain_factor(cap, cal).value == 1.0 @@ -277,35 +354,6 @@ def test_the_oscillator_term_scales_with_the_exponent(tmp_path): assert b.relative == pytest.approx(2 * a.relative, rel=1e-9) -def test_rate_stability_is_zero_for_a_constant_period(tmp_path): - """Not exactly zero: the capture format stores host time to a microsecond, - which is ~0.05 ppm of scatter here. Real drifts run 186 to 1241 ppm, so a - 1 ppm bar separates "nothing" from anything worth reporting.""" - cap = make(tmp_path / "c.csv", rows=4000) - assert calibrate.rate_stability(cap) < 1e-6 - - -def test_rate_stability_recovers_a_planted_ramp(tmp_path): - """The oscillator warms up and slows; this is what measures that.""" - rows, dt, ramp = 8000, 1 / 300.0, 0.002 # 2000 ppm across the run - index = np.arange(rows) - # Period growing linearly, so the instantaneous rate falls by `ramp`. - times = np.cumsum(dt * (1 + ramp * index / rows)) - path = tmp_path / "c.csv" - write_capture(path, rows=rows, dt=dt) - text = path.read_text().splitlines() - header = [l for l in text if l.startswith("#")] + [text[len( - [l for l in text if l.startswith("#")])]] - body = text[len(header):] - rebuilt = header + [ - ",".join([row.split(",")[0], f"{1_700_000_000.0 + t:.6f}"] - + row.split(",")[2:]) - for row, t in zip(body, times)] - path.write_text("\n".join(rebuilt) + "\n") - cap = capture.load(path) - assert calibrate.rate_stability(cap) == pytest.approx(ramp, rel=0.15) - - def test_a_two_point_fit_refuses_to_invent_a_confidence(): """It has no residual, so it cannot estimate its own uncertainty.""" ccs = [100, 400] @@ -527,8 +575,8 @@ def test_a_calibrated_file_is_not_a_capture(tmp_path): def test_the_header_records_the_calibration_that_was_applied(tmp_path): source = tmp_path / "raw.csv" - cap = make(source, counts_per_second=92889.0, overhead=40.61e-6) - cal = calibrate.calibration(88546.0, 40.61e-6) + cap = make(source, counts_per_second=92889.0, clocked=True) + cal = calibrate.calibration(REFERENCE_HZ, REFERENCE_OVERHEAD_S) cal_path = tmp_path / "bench.json" calibrate.save_calibration(cal_path, cal) out = tmp_path / "out.csv" @@ -595,114 +643,3 @@ def test_cli_reports_a_bad_calibration_rather_than_writing_output(tmp_path): with pytest.raises(SystemExit, match="rm3100_calibration"): calibrate.main() assert not out.exists() - - -# -------------------------------------------------------------------------- -# Measuring the oscillator's drift -# -------------------------------------------------------------------------- - -DRIFT_DT = 1.0 / 250.0 -# 102 s: ten 10 s windows with room to spare. Landing exactly on a multiple of -# the window would make the window count depend on whether the planted drift -# happens to push the duration over it. -DRIFT_ROWS = 25_500 - - -def drifting(path, ppm_per_second=0.0, jitter_s=0.0, rows=DRIFT_ROWS, - dt=DRIFT_DT, seed=0): - """A capture whose sample period ramps, with optional host-side jitter.""" - return capture.load(write_capture( - path, rows=rows, dt=dt, counts=field_counts(rows, FIELD, 100, 5.0, seed), - times=drifting_times(rows, dt, ppm_per_second, jitter_s, seed))) - - -def test_window_rates_recovers_a_planted_ramp(tmp_path): - """The rate falls at the planted rate, in fractional terms per second.""" - ppm_per_second = 8.0 - cap = drifting(tmp_path / "d.csv", ppm_per_second=ppm_per_second) - centres, rates = calibrate.window_rates(cap, 10.0) - assert len(rates) == 10 - # A period ramping up as (1 + r*t) is a rate falling as (1 - r*t). - slope = np.polyfit(centres, rates / rates.mean(), 1)[0] - assert slope == pytest.approx(-ppm_per_second * 1e-6, rel=0.05) - - -def test_window_rates_is_flat_without_drift(tmp_path): - cap = drifting(tmp_path / "f.csv") - _, rates = calibrate.window_rates(cap, 10.0) - assert np.ptp(rates) / rates.mean() < 1e-9 - - -def test_window_rates_centres_span_the_capture(tmp_path): - cap = drifting(tmp_path / "c.csv") - centres, rates = calibrate.window_rates(cap, 10.0) - assert len(centres) == len(rates) - assert centres[0] == pytest.approx(5.0) - assert centres[-1] == pytest.approx(95.0) - - -@pytest.mark.parametrize("window", [0.0, -1.0, 60.0, 1e6]) -def test_window_rates_declines_impossible_windows(tmp_path, window): - """Non-positive, or too long to fit two windows in the capture.""" - cap = drifting(tmp_path / "s.csv") - centres, rates = calibrate.window_rates(cap, window) - assert len(centres) == 0 and len(rates) == 0 - - -def test_rate_stability_recovers_a_planted_ramp(tmp_path): - """Peak to peak across eight windows, against the ramp that produced it.""" - cap = drifting(tmp_path / "r.csv", ppm_per_second=8.0) - # Eight windows of 12.5 s: centres 6.25 s and 93.75 s apart, so the spread - # is the ramp over 87.5 s, not over the full 100 s. - expected = 8.0e-6 * cap.duration * (1 - 1 / 8) - assert calibrate.rate_stability(cap) == pytest.approx(expected, rel=0.05) - - -def test_rate_stability_is_zero_on_a_uniform_grid(tmp_path): - assert calibrate.rate_stability(drifting(tmp_path / "u.csv")) < 1e-9 - - -def test_rate_stability_declines_a_capture_too_short_to_window(tmp_path): - cap = drifting(tmp_path / "t.csv", rows=200) - assert calibrate.rate_stability(cap) == 0.0 - - -def test_rate_allan_of_a_linear_drift_rises_with_tau(tmp_path): - """sigma_y(tau) = |D| * tau / sqrt(2) for a deterministic frequency ramp. - - This is the shape the real captures show above a few seconds, and it is what - separates a warming oscillator from a random walk. - """ - ppm_per_second = 8.0 - cap = drifting(tmp_path / "a.csv", ppm_per_second=ppm_per_second) - taus, devs = calibrate.rate_allan(cap, [4.0, 8.0, 16.0]) - assert list(taus) == [4.0, 8.0, 16.0] - for tau, dev in zip(taus, devs): - assert dev == pytest.approx(ppm_per_second * 1e-6 * tau / math.sqrt(2), - rel=0.05) - - -def test_rate_allan_falls_with_tau_when_only_jitter_is_present(tmp_path): - """Independent read jitter averages down; drift does not. Opposite slopes.""" - cap = drifting(tmp_path / "j.csv", jitter_s=2e-3) - taus, devs = calibrate.rate_allan(cap, [4.0, 16.0]) - assert devs[0] > devs[1] * 2 - - -def test_rate_allan_is_negligible_on_a_uniform_grid(tmp_path): - cap = drifting(tmp_path / "z.csv") - _, devs = calibrate.rate_allan(cap, [4.0, 8.0]) - assert np.all(devs < 1e-9) - - -def test_rate_allan_drops_taus_that_yield_too_few_windows(tmp_path): - """Two windows give one difference, which is not an estimate of anything.""" - cap = drifting(tmp_path / "d2.csv") - taus, devs = calibrate.rate_allan(cap, [10.0, 40.0, 200.0]) - assert list(taus) == [10.0] # 40 s gives two windows, 200 s none - assert len(devs) == 1 - - -def test_rate_allan_accepts_an_empty_request(tmp_path): - taus, devs = calibrate.rate_allan(drifting(tmp_path / "e.csv"), []) - assert len(taus) == 0 and len(devs) == 0 diff --git a/tests/test_capture.py b/tests/test_capture.py index 82d8371..79e66a2 100644 --- a/tests/test_capture.py +++ b/tests/test_capture.py @@ -75,11 +75,11 @@ def test_rate_error_is_signed_like_the_rate_not_the_period(tmp_path): assert slow.rate_error == pytest.approx(250.0 / 282.0 - 1, rel=1e-6) -def test_elapsed_bases_differ_by_exactly_the_rate_error(tmp_path): +def test_elapsed_is_the_grid_scaled_by_the_measured_period(tmp_path): cap = capture.load(write_capture(tmp_path / "c.csv")) assert cap.elapsed[0] == 0.0 - ratio = cap.elapsed[-1] / cap.elapsed_nominal[-1] - assert ratio == pytest.approx(cap.dt_true / cap.dt_nominal) + assert cap.elapsed[-1] == pytest.approx(cap.sample_index[-1] * cap.dt_true) + assert cap.dt_true / cap.dt_nominal == pytest.approx(1 / (1 + cap.rate_error)) def test_duration_spans_the_capture(tmp_path):