diff --git a/.gitignore b/.gitignore index 3a4cfd4..158db13 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,6 @@ __pycache__ .venv *.png -*.csv \ No newline at end of file +*.csv + +/Noise_Floor_Testing \ No newline at end of file diff --git a/NOTES.md b/NOTES.md new file mode 100644 index 0000000..8eee567 --- /dev/null +++ b/NOTES.md @@ -0,0 +1,448 @@ +# 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 +corrected rather than deleted when a later measurement supersedes it. + +Every capture cited here predates the current header format and has been +deleted (`*.csv` is gitignored). The numbers stand as measurements; the files +cannot be re-loaded. Re-record before relying on any of it. + +**Rigs.** Two hosts appear below. *Bench* is bare metal. *VM* is a +resource-limited virtual machine with virtualised USB passthrough, shared cores +and background load; every host-side ceiling measured there is pessimistic and +labelled as such. + +--- + +## Superseded conclusions + +Read this first — these appear as firm recommendations in earlier entries and +are no longer right. + +| Earlier conclusion | Superseded by | Now | +|---|---|---| +| "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 | + +--- + +## First working capture + +10 s, cycle count 200, TMRC 0x96, sensor stationary on a desk beside a PC. + +``` +349 samples over 10.02 s -> 34.8 Hz + X: mean +6.473 uT sd 177 nT + Y: mean +50.914 uT sd 106 nT + Z: mean +18.944 uT sd 396 nT + |B|: mean 54.709 uT sd 245 nT + sample interval: 28.78 ms, sd 1.52 ms +``` + +- **Rate.** 34.8 Hz against a nominal ~37 Hz is in spec: §5.2.1 quotes roughly + 7% one-sigma tolerance on the update rate. +- **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. +- **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 + warm-up before trusting absolute values. + +Nothing here is calibrated: these are raw sensor outputs with no hard- or +soft-iron correction, so individual axis values reflect local distortion as much +as Earth's field. + +--- + +## Supply A/B: 3.3 V adapter vs 3.0 V LDO — indistinguishable + +Two 60 s captures on the bench, cycle count 200, TMRC 0x96, 400 kHz. + +| | 3V3 adapter | 3.0 V LDO | +|---|---|---| +| \|B\| | 40,345 nT | 42,731 nT (**+5.9%**) | +| sd \|B\| absolute | 35.4 nT | 38.5 nT (+8.8%) | +| sd \|B\| **fractional** | **878 ppm** | **900 ppm** (+2.5%) | + +**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. + +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 +adapter rail most likely. Supply quality was never the binding constraint; the +magnetic environment is. + +**The 5.9% scale change cannot be attributed to the supply from this data.** The +sensor also moved: the field direction rotated 6.16°, and the per-axis ratios +(X 1.136, Y 1.051, Z 0.870) spread by 26.6%, where a pure gain change would +scale all three identically. |B| is preserved under rotation but not under +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 +3.0 V, noting that performance "will differ" elsewhere. At 3.3 V the gain is an +extrapolation off the calibration point. + +This experiment is what `compare.py` was written for; it refuses to attribute a +magnitude difference to gain unless both movement checks pass. + +### Designing an A/B supply test + +- **Clamp the sensor and never touch it.** Swapping supplies must not disturb + it. This is the single thing that invalidated the first attempt. +- **Interleave A/B/A/B**, never all of A then all of B. Drift over a session is + large enough to swamp the effect, and only interleaving lets adjacent pairs be + differenced. +- **Repeat 3–4 times per configuration.** With one run each, a 2.5% difference + is indistinguishable from run-to-run variation. +- **Let it thermally settle** before each run. Coil resistance moves 0.4 %/°C + (Table 3-3), so an unsettled sensor drifts through the measurement. +- Record the configuration with `--note`, so the comparison does not depend on + filenames surviving: + +```bash +./.venv/bin/python logger.py --duration 300 --note "3V3" --output a1.csv +./.venv/bin/python logger.py --duration 300 --note "3V0 LDO" --output b1.csv +./.venv/bin/python compare.py --group note *.csv +``` + +--- + +## 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**. + +It is 60 Hz aliasing, proven by making it move: + +| capture | fs | observed peak | predicted 60 Hz alias | +|---|---|---|---| +| 3V3 | 34.7565 Hz | 9.5063 Hz | 9.5129 Hz | +| 3.0 V LDO | 34.5557 Hz | 9.0725 Hz | 9.1114 Hz | +| 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 +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. + +**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 +75 Hz it folds to 9.304 Hz against a 9.233 Hz prediction, with only 4 nT/√Hz +left at 60 Hz because that frequency is no longer sampled. This is why the +default configuration runs at ~282 Hz: at that rate mains *and* its second +harmonic are inside the band and can be notched afterwards, which no amount of +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 +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. +Hourly |B| sd tells the story: + +| hours | \|B\| sd | usable | +|---|---|---| +| 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 | + +A 5-minute slice from the middle of the quiet window gave |B| = 42,132 nT with +sd 32.5 nT (771 ppm). + +**2. A single rate does not describe a 13-hour capture.** One slope across the +whole run leaves a residual sd of 5.3 s and a 22.8 s span. Hour-by-hour fits +show the oscillator speeding up monotonically: + +| hour | fitted period | vs run mean | +|---|---|---| +| 2 | 29.0187 ms | +1254 ppm | +| 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 +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 +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 +bottoms at ~13 nT near τ = 0.12 s and climbs steadily after. + +--- + +## Host-side rate ceilings + +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 +START**, so `ch347.write_read()` does it in one. Measured at 400 kHz: + +| | separate | combined | +|---|---|---| +| `data_ready()` | 0.652 ms | **0.350 ms** | +| `read_raw()` | 0.838 ms | **0.548 ms** | + +### 2. The GIL, not the OS scheduler + +Even after that, `logger.py` stalled ~10 ms at 150 Hz while the *same sampling +loop with no writer thread* peaked at 3.1 ms (p99.9 = 2.9 ms, zero threshold +crossings in 2079 samples). The sampler releases the GIL inside each USB +transfer and must re-acquire it, waiting up to Python's default **5 ms** switch +interval while the writer holds it. Setting `sys.setswitchinterval(0.0005)` cut +the worst bracket from 9.96 ms to 5.65 ms and made 150 Hz pass. + +Worth remembering generally: with a latency-sensitive loop and a helper thread, +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%**. + +### Measured ceilings (VM, cycle count 50, 400 kHz) + +| TMRC | Nominal | Result | Worst DRDY bracket | +|---|---|---|---| +| 0x92 | 600 Hz | fail | 1.86 ms of 1.67 ms (112%) | +| 0x93 | 300 Hz | fail | 3.45 ms of 3.33 ms (104%) | +| 0x94 | 150 Hz | borderline | 5.65 / 5.95 / 6.53 / 7.19 ms of 6.67 ms | +| 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, +> 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 +> run the capture. The "% of margin used" line `logger.py` prints after every +> run is the number to watch. + +### The bus is never the constraint + +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 | +|---|---|---| +| 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 | — | + +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 +is always achievable. Every miss observed has been a host stall. + +### Bus speed, three runs each at cycle count 228 + +| bus | bus time | host cost | implied USB overhead | +|---|---|---|---| +| 100 kHz | 1.500 ms | 2.03–2.19 ms | ~0.63 ms | +| 400 kHz | 0.375 ms | 0.88–0.95 ms | ~0.55 ms | +| **750 kHz** | 0.200 ms | **0.78–0.80 ms** | ~0.59 ms | + +Host cost is `bus time + ~0.6 ms` of fixed USB round-trip latency, so a faster +bus helps but cannot go below that floor. An earlier note claimed 400 kHz +benchmarked faster than 750; that rested on a single short run whose 4% +difference was inside the run-to-run spread and does not survive repetition. +400 kHz remains a sane fallback if 750 proves marginal on long wiring. + +--- + +## 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 +the DRDY bracket stays under 10 ms. The first implementation thresholded the +interval and aborted a perfectly healthy capture at sample 88. + +**Period calibration is worth 100×.** Counting how many grid points passed +unseen inside a bracket needs the real period, and the TMRC table value is 6–9% +out on this unit: + +| | period | error | +|---|---|---| +| nominal (TMRC table) | 27.03 ms | −6.4% | +| **calibrated (1 s)** | **28.8849 ms** | **0.068%** | +| whole-run fit (25 s) | 28.8653 ms | — | + +**Counting needs the interval, detection needs the bracket.** Fit residual after +adopting the two-signal scheme: + +| capture | lost | residual | +|---|---|---| +| 37 Hz clean | 0 | **0.019 periods** | +| 150 Hz clean | 0 | 0.082 periods | +| 150 Hz under CPU load | 0 | 0.444 periods | +| 300 Hz | 10 | 1.241 periods | + +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)` +estimator would have inserted a spurious placeholder and slipped the index from +there on. At 300 Hz the ambiguous counter fired 19 times. + +**Run-to-run stability.** Fitted periods land 6.0–6.6% below nominal, and two +consecutive 60 s runs agreed to **0.017%**, so a single-slope drift correction +is sound over at least a minute. + +--- + +## Rate sweep + +`sweep.py` end to end. Every target from 1 to 128 Hz landed within 0.1% +(columns as of the run; the current tool also prints duty): + +``` + target measured err cc LSB/uT nT/count noise bus use + 1H 1.00H +0.0% 29769 10929.7 0.091 1.21n 0.0% + 8H 8.01H +0.1% 3718 1366.4 0.732 3.41n 0.3% + 32H 32.00H -0.0% 927 341.8 2.926 6.83n 1.2% + 128H 127.99H -0.0% 229 85.6 11.687 13.75n 4.8% + 256H 249.23H -2.6% 113 43.0 23.265 19.57n 9.3% 2 lost +``` + +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. + +**Cycle count sets the ceiling, not TMRC.** Measured at 400 kHz: + +| Cycle count | Measured | Table 3-1 predicts | +|---|---|---| +| 200 | 152 Hz | 147 Hz | +| 100 | 297 Hz | 283 Hz | +| 50 | 538 Hz | 533 Hz | + +**The manual contradicts itself here and Table 3-1 is the correct one.** The +§5.2.1 note claims cycle count 200 allows a "~430 Hz" three-axis update rate; +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 +TMRC 0x92 and 73.86 Hz at 0x94. Varying only cc at TMRC 0x94: + +| cc | measured | governed by | duty | +|---|---|---|---| +| 50 | 138.53 Hz | TMRC | 25% | +| 200 | 138.84 Hz | boundary | 99% | +| 400 | 73.83 Hz | cycle count | 100% | +| 800 | 37.09 Hz | cycle count | 100% | + +**Idling costs sensitivity.** Reaching ~130 Hz two ways: + +| | rate | cc | duty | ASD | +|---|---|---|---|---| +| 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, +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 +different rates: + +``` +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 +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 +the time and quantising 3.8× more coarsely. Its extra sensor speed is capability +the host cannot collect. + +--- + +## 60 s capture at maximum rate (office rig, historical) + +28,864 samples, 481 Hz, cycle count 50, 400 kHz. The sensor was nudged by hand +at t ≈ 41 s, so statistics are split around that: + +| Axis | sd, 0–40 s | sd, 47–60 s | step across the event | +|---|---|---|---| +| X | 84 nT | 85 nT | +280 nT | +| Y | 388 nT | 376 nT | +217 nT | +| Z | 187 nT | 198 nT | +232 nT | +| \|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 + 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 + 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 + Table 3-1 quotes for cycle count 50; Y is ~13×. +- **Timing was jittery**: intervals averaged 2.08 ms (481 Hz) with sd 1.09 ms + and a 19.6 ms worst case. This capture predates the chip-grid time base, so it + is unsuitable for spectral work without resampling. A capture taken today + carries an exact grid regardless of read jitter. + +--- + +## Supply ripple: what the manual does and does not say + +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. + +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 +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 +noise near the ~180 kHz LR oscillation. + +Unverified on this rig: whether the breakout carries local decoupling, and what +the CH347's 3V3 rail ripple actually measures against the 50 mV limit. + +--- + +## Open questions + +- Re-run the LDO comparison with the sensor **clamped**, to separate gain from + movement. This is the one experiment the A/B entry above failed to perform. +- Deal with 60 Hz coupling at source — standoff distance, then a repeat + characterisation — rather than only sampling fast enough to see it. +- Check the dither assumption at cycle count 100 against a genuinely quiet site. + `characterize.py` prints `sd/LSB` per axis; the concern is that a sensor + quieter than spec has less dither than the 0.79 LSB the spec figure implies. +- Reliability over extended runs on bare metal, where the VM ceilings above do + not apply. diff --git a/README.md b/README.md index 9fb857f..ff64bc5 100644 --- a/README.md +++ b/README.md @@ -1,23 +1,36 @@ # RM3100 logging over a CH347 USB-I2C adapter -Logs a PNI RM3100 geomagnetic sensor connected to the I2C pins of a Waveshare -USB to UART/I2C/SPI/JTAG adapter (CH347, USB `1a86:55db`). +Logs a PNI RM3100 geomagnetic sensor connected to the I²C pins of a Waveshare +USB to UART/I2C/SPI/JTAG adapter (CH347, USB `1a86:55db`), with a time base good +enough for spectral work: the sample index is a chip-clock grid coordinate, not +a count of host reads. -Reference: *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. +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. + +## Quickstart + +```bash +./setup.sh # udev rule + venv, sudo for the rule only +./.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 +``` ## Hardware Waveshare adapter in **Mode 1** (UART1 + I2C + SPI), voltage selector at **3V3**. -The RM3100 is on a breakout board. - -> 3V3 is required, though 3.0v would be ideal for the analog side. - -### Wiring - -All connections below are **confirmed working** — this exact wiring produces -valid data. +The RM3100 is on a breakout board. All connections below are confirmed working. | Adapter | RM3100 | Notes | |---|---|---| @@ -31,92 +44,57 @@ valid data. 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). -The two supply pins and I2CEN are the non-obvious ones; both are explained below, -because omitting either produces a failure that is easy to misdiagnose. +Four pins are worth knowing about; two of them cause failures that are easy to +misdiagnose. -#### I2CEN (pin 22) must be tied HIGH +**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. -`HIGH = I2C, LOW = SPI` (§4.3.1). The RM3100's SDA/SCL pins are shared with the -SPI MOSI/SCLK pins, and I2CEN is what selects between the two interfaces. 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. +**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. -At DVDD = 3.3 V it needs at least `VIH = 0.7 × DVDD` = 2.31 V, so tie it -directly to 3V3. +**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. -#### AVDD/AVSS (pins 4, 5) must be powered to take measurements +**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). -§4.3.1: *"AVDD can be turned off when not making a measurement to conserve -power, since all other operations are supported with DVDD."* - -So with only DVDD connected, the part is **half alive** in a way that is easy to -misread as working: the I2C address responds, REVID reads back `0x22`, and -registers write and read back correctly — but the analog section that drives the -sensor coils is dead, so DRDY never asserts and the measurement registers stay -at zero. `BIST` reports `XOK=YOK=ZOK=0` (see Diagnostics below), which is the -unambiguous tell. - -**Both rails are needed — AVDD is not a substitute for DVDD.** Two constraints -(§4.3.1, Table 3-5): -- DVDD must come up **before or at the same time as** AVDD, never after. Moving - the supply from DVDD to AVDD instead of adding it violates this. -- AVDD must stay within **0.1 V** of DVDD while on — sharing one 3V3 rail - satisfies this for free. - -#### Power quality - -The manual specifies very little here, but what it does specify is firm. +### Power | Parameter | Limit | Source | |---|---|---| -| Ripple on AVDD or DVDD | **50 mV peak-to-peak** | Table 3-5, `VDD_ripple` | -| DVDD − AVDD while running | **±0.1 V** | Table 3-5, `ΔVDD_OP` | +| Ripple on AVDD or DVDD | 50 mV peak-to-peak | Table 3-5, `VDD_ripple` | +| DVDD − AVDD while running | ±0.1 V | Table 3-5, `ΔVDD_OP` | | Supply range | 2.0–3.6 V (typ 3.0) | Table 3-5 | | Absolute maximum | 3.7 V | Table 3-4 | -The ripple limit carries **no frequency qualifier** — it is stated flat, so -switching noise and mains hum are not distinguished. +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 +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. -Reference decoupling (Figures 4-1 and 4-2) is a **10 µF bulk capacitor in -parallel with 0.1 µF ceramic**, rail to ground. One rail feeds AVDD, both DVDD -pins, and I2CEN. That is the entire filtering specification — no ferrite, no -split analog/digital 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 (ferromagnetic packaging), not an electrical one. - -**The manual gives no transfer function from ripple to field error** — no -µT-per-mV, no PSRR, no ripple-vs-noise curve. Only the 50 mV limit. - -Slow drift should matter less than fast ripple. §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 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 is rejected far better than noise -near the ~180 kHz LR oscillation. (That last step is inference from the -described architecture, not an explicit claim in the manual.) - -**This rig runs at 3.3 V, but every number in Table 3-1 is quoted at 3.0 V.** -Footnote 1: *"Other bias resistors, external timing resistors and operating -voltages may be used, but performance will differ from the values listed."* -3.3 V is comfortably in range, but the 75 LSB/µT gain this driver uses is -strictly a 3.0 V figure — a candidate scale-factor error the manual does not -quantify. Unverified here: whether the breakout carries local decoupling, and -what the CH347's 3V3 rail ripple actually measures against the 50 mV limit. - -#### SA0 / SA1 set the I2C address - -The top 5 address bits are fixed at `0b01000`; SA0 (pin 3) and SA1 (pin 28) set -the low two, giving 0x20–0x23 (§4.5). On this breakout both are strapped high, -so the sensor answers at **0x23**. These 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 (§5.4.1), which the -manual explicitly offers as an alternative. +**This rig runs at 3.3 V, but every number in Table 3-1 is quoted at 3.0 V** +(footnote 1: *"performance will differ"*). 3.3 V is comfortably in range, but the +gain fit behind `tesla_per_count` is strictly a 3.0 V figure, so the scale factor +is an extrapolation off the calibration point. A 3.0 V LDO was built and compared +— see [NOTES.md](NOTES.md), which also covers what the manual does and does not +say about ripple. ## Setup @@ -124,47 +102,60 @@ manual explicitly offers as an alternative. ./setup.sh ``` -Idempotent. It installs a udev rule, creates `.venv`, and installs pyusb. It -needs `sudo` for the udev rule only. +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. -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: +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: ``` -SUBSYSTEM=="usb", ATTRS{idVendor}=="1a86", ATTRS{idProduct}=="55db", TAG+="uaccess", GROUP="plugdev", MODE="0660" +# /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" ``` -- **`TAG+="uaccess"`** makes systemd-logind put an ACL on the node for whoever - is logged in at the local seat. No group, no logout, and it is the only part - that works unaided on atomic Fedora. `setup.sh` triggers with - `--action=add`, not the default `change`, because the `uaccess` builtin only - runs on add — a change event installs the rule but grants nothing until the - next replug. -- **`GROUP=`/`MODE=`** is the fallback for ssh sessions and seatless systems, - where there is no local seat for `uaccess` to grant to. The group is chosen - from those that exist — `plugdev` is a Debian convention and absent on - Fedora, and naming a group that does not exist makes udev log an error and - leave the node owned by root. If neither `plugdev` nor `dialout` exists the - clause is omitted entirely. +`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. -#### Atomic Fedora (Bazzite, Silverblue, Kinoite) +**Two details in that file are load-bearing**, and getting either wrong produces +a rule that installs cleanly and grants nothing: -`usermod -aG dialout $USER` fails there with *"group 'dialout' does not -exist"* even though `getent group dialout` finds it. Fedora ships -`nss_altfiles`, so `getent` reads both `/etc/group` and `/usr/lib/group`, while -`usermod` writes `/etc/group` alone — and on an rpm-ostree system the group -frequently lives only in the latter. Copy the line across first: +- **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. -```bash -grep -E '^dialout:' /usr/lib/group | sudo tee -a /etc/group -sudo usermod -aG dialout $USER # then log out and back in -``` +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. -None of that is needed for a local desktop session: `uaccess` covers it, which -is why the rule carries both. `setup.sh` tests the node itself at the end and -prints whichever fix applies rather than guessing up front. +Two more atomic-Fedora traps `setup.sh` prints fixes for: -Python 3.14 here is PEP 668 externally-managed, hence the venv rather than a +- `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`. ## Usage @@ -178,34 +169,71 @@ 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 the rate and the LSB; see Recommended configuration | -| `--tmrc` | fastest | rate register. Default lets the cycle count set the rate; give one only to sample *slower* | +| `--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 scan | -| `--bus-speed` | 750 | I²C kHz; 100 would spend 42% of the period on the bus at cc=100 | +| `--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 | +| `--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 | -Exit status is non-zero if an interval was missed. The partial capture is kept -and is valid up to that point. +**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. -### Capture format +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. + +### 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 --start 60 --end 360 +./.venv/bin/python compare.py --group note *.csv # A/B two conditions +./.venv/bin/python sweep.py --rates 10,50,150 # measure what each rate delivers +``` + +- **`plot.py`** — four stacked panels with independent y-scales, since the three + axes sit at very different DC offsets. +- **`characterize.py`** — amplitude spectral density against the 1.2 nT/√Hz of + Table 3-1, Allan deviation, residual distribution, and host read latency. The + `sd/LSB` column is what settles whether the quantiser is properly dithered. +- **`compare.py`** — reports noise **fractionally** (ppm of |B|), because a gain + change carries the noise with it and an absolute comparison reads a pure scale + change as a noise difference. It also refuses to attribute a magnitude + difference to gain unless the sensor demonstrably held still. +- **`sweep.py`** — runs `logger.py` once per target rate and tabulates what each + configuration actually achieved. + +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. + +## 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. Everything else -(chip time, elapsed time, tesla, magnitude) is reconstructed on load by -`capture.py`, so there is exactly one source of truth for each quantity. +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. ``` # rm3100_capture: 1 -# nominal_rate_hz: 37.0 -# tmrc: 0x96 -# cycle_count: 200 -# tesla_per_count: 1.3347570742124934e-08 -# calibrated_period_s: 0.028891621... +# nominal_rate_hz: 282.53110196547465 +# tmrc_nominal_hz: 600.0 +# tmrc: 0x92 +# cycle_count: 100 +# tesla_per_count: 2.617115938236064e-08 # i2c_address: 0x23 -# bus_speed_khz: 400 +# bus_speed_khz: 750 # revid: 0x22 +# calibrated_period_s: 0.0035391156462585034 +# note: bench, 3V0 LDO sample_index,system_time_unix,x_raw,y_raw,z_raw,warning 0,1755930856.722866,-2765,767,378, 5,1755930856.867436,0,0,0,MISSED @@ -213,10 +241,15 @@ sample_index,system_time_unix,x_raw,y_raw,z_raw,warning 21,1755930857.272108,-2761,769,377,AMBIGUOUS ``` +`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. + The `warning` column carries space-separated flags, empty when the row is fine, -and the column 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. +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. | flag | meaning | |---|---| @@ -233,185 +266,214 @@ likely — 1.35 probably is one period and 1.65 probably two. Deciding which nee neighbouring timestamps and assumptions and can still be wrong, which is exactly why these rows are flagged rather than silently resolved. -`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`. Rows are ~33 bytes, -against ~90 for the old redundant schema. +`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. -All times are unix epoch seconds; no ISO-8601 anywhere in the data. +## How it works -```bash -./.venv/bin/python plot.py capture_60s.csv # -> capture_60s.png -./.venv/bin/python plot.py capture.csv --smooth 0 # raw only, no moving average -``` +**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: -### One-step capture and analysis +- **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. -```bash -./noise_floor_test.sh --duration 60 # capture, then plot + characterize -./noise_floor_test.sh # until Ctrl-C, then analyse -``` - -Activates the venv and forwards every flag to `logger.py`, then runs `plot.py` -and `characterize.py` on the result. It defaults `--bus-speed 400` (measured -best) and picks the output name so the analysis can find it; anything you pass -explicitly wins. - -Analysis runs whether the capture ended on `--duration`, on **Ctrl-C**, or on an -abort — a partial capture is still valid up to the abort. Exit status is -non-zero if any stage failed, so it is usable from a scheduler. - -### Timing model - -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 is uniform even when - host reads are late, so `sample_index` is a grid coordinate, not a count of - reads. -- **A late read returns the *newer* measurement**, not a delayed one. Missing an - interval therefore silently skips a grid point and corrupts the time axis. - -So `logger.py` **records a lost measurement rather than dropping it**: placeholder -rows carrying `MISSED` in the count fields keep `sample_index` contiguous, so the -gap is explicit instead of silently compressing the timeline. The run continues -and reports the total at the end, exiting non-zero. Double-counting is -structurally impossible: `HSHAKE` DRC1=1 makes DRDY clear when the results are -read, and `configure()` verifies the readback. +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 inside a bracket 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 against real time, which is exactly what -placeholders exist to prevent. `calibrate_period()` takes a least-squares slope -over one second of loss-free samples: - -| | period | error | -|---|---|---| -| nominal (TMRC table) | 27.03 ms | −6.4% | -| **calibrated (1 s)** | **28.8849 ms** | **0.068%** | -| whole-run fit (25 s) | 28.8653 ms | — | - -A 100× improvement, from ~35 points at ~0.5 ms jitter. It is then refined from -clean intervals during the run so it follows thermal drift. - -**If no clean second can be found, the run aborts before recording anything** — -without a clean stretch the period cannot be measured, so misses cannot be -counted and the index would not track chip time. That is also the honest signal -that the requested rate is not sustainable. +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 **bracket**, rigorously: a bracket shorter than one - period cannot contain two completions, whatever the timing precision. +- *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 alone cannot say, since it only reaches - back to the last poll that saw DRDY clear and so discards where the grid 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. -With the period calibrated to 0.07% and each completion located to under 1% of -a period, the rounding is unambiguous for small gaps. Measured fit residual: - -| capture | lost | residual | -|---|---|---| -| 37 Hz clean | 0 | **0.019 periods** | -| 150 Hz clean | 0 | 0.082 periods | -| 150 Hz under CPU load | 0 | 0.444 periods | -| 300 Hz | 10 | 1.241 periods | - -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)` -estimator would have inserted a spurious placeholder and slipped the index. - -`logger.py` reports an **ambiguous** count when a rounding lands near a -half-period. That is the signal the count itself is a guess, and it is distinct -from merely having missed something — at 300 Hz it fired 19 times. A capture -with ambiguous gaps should be re-recorded before spectral work. - -Detection uses the **DRDY bracket** — the span between the last poll showing DRDY -clear and the poll showing it set — not the read-to-read interval. The bracket 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. - -The interval cannot make that claim, and measurement settled it: on this rig the -read-to-read interval reaches **37 ms against a 28.8 ms period** — a 29 % overshoot -from host stalls alone — while the bracket stays under 10 ms. An interval -threshold aborts on perfectly healthy captures; the first attempt here did -exactly that at sample 88. +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 -chip's RC oscillator is *regular* but only accurate to ±7% (§5.2.1) — this unit -runs ~6% slow, 34.7 Hz against a 37 Hz nominal. `capture.py` regresses the host -clock on `sample_index` to recover the true period: +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 ~6% wrongly scaled +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 otherwise displace every spectral feature by 6%. +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. -### Sample rate +**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. -**Cycle count, not TMRC, sets the ceiling.** Measured with `--bus-speed 400`: +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. -| Cycle count | Measured (polled) | Table 3-1 predicts | Noise (Table 3-1) | -|---|---|---|---| -| 200 | 152 Hz | 147 Hz | 15 nT | -| 100 | 297 Hz | 283 Hz | 20 nT | -| 50 | 538 Hz | 533 Hz | 30 nT | +## Configuration -**The manual contradicts itself here and Table 3-1 is the correct one.** The -§5.2.1 note claims cycle count 200 allows a "~430 Hz" 3-axis update rate; -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. +**Cycle count is the rate knob, not TMRC.** Two ceilings compete and the slower +one wins (§5.2.1): the cycle count sets how *long* a measurement takes, TMRC sets +how *often* one is started. -Bus speed matters less than expected: 400 kHz measured best (538 Hz), 100 kHz -gave 444 Hz, and 750 kHz was slightly *worse* than 400 kHz at 515 Hz. Use 400. +| | sets | granularity | +|---|---|---| +| cycle count | `3 × (cc/90,000 + 68.7 µs)` per measurement | continuous | +| TMRC | how often a measurement starts | factor-of-two steps | -Skipping the DRDY poll and free-running reads reaches 1100–1460 Hz, but 62–87% -of those samples are duplicates, so it buys nothing real. +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. -### Threading +**Cycle count also sets resolution**, which is the real trade: -The sampling loop does I2C and nothing else — it pushes raw counts onto a queue, -and a writer thread does unit conversion, CSV formatting, flushing and the -console display. Wall-clock timestamps are reconstructed in the writer from one -start time plus each sample's monotonic offset, rather than calling -`datetime.now()` per sample; that is both cheaper and immune to an NTP step -mid-capture. +| cc | rate | Nyquist | nT/LSB | dither | duty | | +|---|---|---|---|---|---|---| +| 50 | 534 Hz | 267 Hz | 50.37 | 0.58 | 89.0% | spectrum, thin dither | +| **100** | **283 Hz** | **141 Hz** | **26.17** | **0.79** | **94.2%** | **default** | +| 200 | 145 Hz | 73 Hz | 13.35 | 1.10 | 97.0% | resolution, 60 Hz only | +| 400 | 74 Hz | 37 Hz | 6.74 | 1.54 | 98.5% | resolution, mains aliases | -This is worth roughly **28%**: the same 60 s capture yielded 22,558 samples -(376 Hz) with inline writing and 28,864 (481 Hz) with the writer split off. -The queue is bounded, and any drop is counted and reported rather than -silently losing data. +*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. + +**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 one assumption here not yet + checked against a quiet-site capture. + +Fall back to cycle count 200 if `characterize.py` reports `sd/LSB` below ~0.35. + +**Bus speed is independent of the rate** — it appears in neither mechanism. What +it sets is latency: how long a read takes, hence how tightly DRDY can be +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. + +| bus | traffic/sample | share of a cc=100 period | +|---|---|---| +| 100 kHz | 1.500 ms | **42%** | +| 400 kHz | 0.375 ms | 11% | +| **750 kHz** | **0.200 ms** | **6%** | + +Below ~0.46 Hz the 16-bit cycle-count register runs out and TMRC must set the +cadence; below cycle count 30 the manual warns of quantisation (§5.1). Both +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. +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 | +| 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 | +| cycle count past 400 | `cc > MAX_SPEC_CYCLE_COUNT` | Table 3-1 ends there, so the printed gain and noise are extrapolated | +| bus over half the period | traffic / period > 50% | names the speed that would fit | + +The aliasing one is the easiest to walk into: `--rate 32` yields a clean-looking +2.9 nT/LSB capture with 60 Hz mains sitting at 4.0 Hz, indistinguishable from +signal. Sampling fast and decimating afterwards gives the same noise floor with +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. +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 | +| `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` | +| `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 | +| `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 | + +**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. + +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 +than merely "data arrived". ## Files | File | Role | |---|---| -| `setup.sh` | udev rule, venv, pyusb | -| `noise_floor_test.sh` | capture then analyse in one step — wraps `logger.py`, then runs `plot.py` and `characterize.py` | -| `ch347.py` | CH347 I2C transport — USB only, knows nothing about the sensor | -| `rm3100.py` | RM3100 driver — registers and decoding, knows nothing about USB | -| `logger.py` | CLI: scan → identify → configure → log, aborting on a missed interval | +| `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 | +| `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 | | `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 | -| `sweep.py` | sweep target rates, reporting measured rate, LSB, noise and bus use per point | +| `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 | | `diagnose-comms.py` | walks USB → bus → identity → registers → BIST → live read, stopping at the first failure | +| `tests/` | pytest suite; needs no hardware | -The adapter/sensor split is deliberate: `ch347.py` is a general I2C master usable +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)`. +`write(addr, data)` and `read(addr, count)` — optionally `write_read()`, which it +prefers when available. ## Implementation notes @@ -419,20 +481,18 @@ 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 and carries I2C/SPI/JTAG with no kernel driver bound, so 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. +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 in `ch347.py` follows the -`aystarik/ch347-i2c-spi-gpio` Linux driver, whose id-table entry +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: @@ -442,592 +502,76 @@ Wire format — bulk OUT `0x06`, bulk IN `0x86`, max 63 bytes per transfer: | 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, byte 0 -is the address ACK and the rest is payload. +**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 read is STOP-then-START, not a repeated START** (§4.5.2, §5.8.4), which is - exactly what two separate CH347 transactions produce. No special handling. +- **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. + `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. -## Diagnostics - -Symptoms map cleanly onto causes, so work down this list. - -| Symptom | Cause | -|---|---| -| `Cannot claim CH347 interface 2: Access denied` | udev rule missing or not yet applied to this node — run `./setup.sh`, which diagnoses which of `uaccess` or group membership fell through. Do **not** reach for `sudo`: it works as your normal user, and masking a permissions problem with root only defers it. `--high-priority` is the one case that legitimately needs root, and it hands the capture file back to `$SUDO_UID` afterwards | -| `No CH347 adapter found` | not plugged in, or not in Mode 1 | -| `ABORTED: missed N interval(s)` | the host could not keep up — raise `--tmrc` (slower) or `--bus-speed`. Expected at cycle count 50 near the sensor's maximum | -| `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; both of these were hit during bring-up | -| 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, which points at AVDD or the -REXT timing resistor rather than anything on the I2C side. - -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 gain is wrong rather -than merely "data arrived". - -## Observed performance - -First working capture, 10 s at TMRC 0x96, cycle count 200, sensor stationary on -a desk next to a PC: - -``` -349 samples over 10.02 s -> 34.8 Hz - X: mean +6.473 uT sd 177 nT - Y: mean +50.914 uT sd 106 nT - Z: mean +18.944 uT sd 396 nT - |B|: mean 54.709 uT sd 245 nT - sample interval: 28.78 ms, sd 1.52 ms -``` - -- **Rate.** 34.8 Hz against a nominal ~37 Hz is in spec: §5.2.1 quotes roughly - 7% one-standard-deviation tolerance on the update rate. -- **Magnitude.** 54.7 µT sits in the expected 25–65 µT band for Earth's field, - which is the real confirmation that gain and the int24 decode are right. -- **Noise** is well above the 15 nT that Table 3-1 quotes for cycle count 200, - and is not a driver problem. Most likely proximity to a PC and the USB - adapter; a secondary candidate is ripple on the CH347's 3V3 rail against the - 50 mV limit (see Power quality above), which has not been measured. Expect far - better readings away from mains wiring and switching supplies. -- **Drift.** All three axes drifted monotonically down over the 10 s (Z most, by - about 1.3 µT), so |B| fell from 55.5 to 54.3. Consistent with thermal settling - after AVDD is first powered — coil DC resistance moves 0.4 %/°C (Table 3-3). - Worth allowing a warm-up period before trusting absolute values. - -None of this is calibrated: the figures are raw sensor output with no hard- or -soft-iron correction, so the individual axis values reflect local distortion as -much as Earth's field. - -### Bench results (desk, cc=200, TMRC 0x96, 400 kHz) - -Timing held up well throughout: fitted periods land 6.0–6.6% below nominal, and -two consecutive 60 s runs agreed to **0.017%**, so a single-slope drift -correction is sound over at least a minute. - -#### Supply: 3.3 V adapter vs 3.0 V LDO — indistinguishable - -Two 60 s captures, `..._DESK_3v3.csv` and `..._DESK_3vLDO.csv`: - -| | 3V3 adapter | 3.0 V LDO | -|---|---|---| -| \|B\| | 40,345 nT | 42,731 nT (**+5.9%**) | -| sd \|B\| absolute | 35.4 nT | 38.5 nT (+8.8%) | -| sd \|B\| **fractional** | **878 ppm** | **900 ppm** (+2.5%) | - -**Compare noise fractionally, not absolutely.** The measured scale moved 5.9% -between these runs, and any scale change carries the noise with it — so the -absolute sd difference is mostly an artifact. Normalised, the two supplies differ -by 2.5%, which is 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 adapter rail, most likely. Supply quality was never the binding constraint; -the magnetic environment is. - -**The 5.9% scale change is not attributable to supply from this data.** The -sensor also moved between captures — the direction rotated **6.16°**, and the -per-axis ratios (X 1.136, Y 1.051, Z 0.870) spread by **26.6%**, whereas a pure -gain change would scale all three identically. Translation through a field -gradient changes |B| just as a gain change does, so the two are confounded here. - -There is still good 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 **3.0 V**, noting performance "will differ" elsewhere. -At 3.3 V the gain is an extrapolation off the calibration point. - -To settle it, clamp the sensor so it cannot move and switch only the supply. -Two diagnostics then make the answer unambiguous: a pure gain change gives -**equal per-axis ratios** (spread ≈ 0) and **≈ 0° of rotation**. - -`compare.py` runs both checks and refuses to attribute a magnitude difference to -gain unless they pass: +## Tests ```bash -./.venv/bin/python logger.py --duration 300 --note "3V3" --output a1.csv -./.venv/bin/python logger.py --duration 300 --note "3V0 LDO" --output b1.csv -./.venv/bin/python compare.py --group note *.csv +./.venv/bin/python -m pytest # whole suite, no hardware needed +./.venv/bin/python -m pytest -k logger # one module ``` -Use `--note` to record the configuration in the header, so the comparison does -not depend on filenames surviving. - -### Designing an A/B supply test - -- **Clamp the sensor and never touch it.** Swapping supplies must not disturb - it. This is the single thing that invalidated the first attempt. -- **Interleave A/B/A/B**, never all of A then all of B. Drift over a session is - large enough to swamp the effect, and only interleaving lets you difference - adjacent pairs. -- **Repeat at least 3–4 times per configuration.** With one run each, a 2.5% - difference is indistinguishable from run-to-run variation. -- **Pick cycle count to match a TMRC step** — see below. -- **Let it thermally settle** before each run. Coil resistance moves 0.4 %/°C - (Table 3-3), so an unsettled sensor drifts through the measurement. - -#### 60 Hz mains aliases into the band — the dominant feature - -The largest spectral feature by far is a line near 9 Hz reaching ~100–125 nT/√Hz, -against a ~4–5 nT/√Hz broadband floor: **25× the noise floor**. - -It is 60 Hz folding in, proven by making it move: - -| capture | fs | peak | predicted 60 Hz alias | -|---|---|---|---| -| 3V3 | 34.7565 Hz | 9.5063 Hz | 9.5129 Hz | -| 3.0 V LDO | 34.5557 Hz | 9.0725 Hz | 9.1114 Hz | -| 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 a -second-order fold, exactly as predicted. - -The RM3100 has **no anti-alias filter** — its cycle-count integration is a boxcar -with poor stopband — so out-of-band interference folds in freely. Until this is -dealt with (notch in post, a sample rate that folds mains somewhere harmless, or -removing the coupling) it will dominate every sd figure and mask any supply-level -improvement. - -#### Overnight run: 13.2 h, 1.63 M samples - -`rm3100_20260823_014804.csv`. **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, and no missed intervals (worst -DRDY bracket 11.0 ms of 27.03 ms, 41% of margin). Two problems in the data, both -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. -Hourly |B| sd tells the story plainly: - -| hours | \|B\| sd | usable | -|---|---|---| -| 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 | - -The quiet window is extracted as `overnight_quiet.csv` (hours 2.0–10.5, -1,055,585 samples). A 5-minute slice from its middle gives |B| = 42,132 nT with -sd 32.5 nT (771 ppm). - -**2. A single rate does not describe a 13-hour capture.** Fitting one slope -across the whole run leaves a residual sd of **5.3 s** and a **22.8 s span**. -Hour-by-hour fits show the oscillator speeding up monotonically: - -| hour | period | vs run mean | -|---|---|---| -| 2 | 29.0187 ms | +1254 ppm | -| 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 -expected from an RC oscillator. The consequence is that **long captures cannot -be given a single frequency axis**; spectra must be computed on windows short -enough that the rate is constant. `capture.py` now warns when the residual sd -exceeds one sample period, and reports read jitter separately from fit residual -— 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 have read as a -20,000× worse host. - -#### Long run: drift-limited beyond a second - -`rm3100_20260823_014804.csv`, 26,604 samples over 769 s: sd rises to 214–284 nT -per axis 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 -bottoms at ~13 nT near τ = 0.12 s and climbs steadily after. Over 13 minutes, -thermal and environmental drift dominate everything else. - -### Achievable rate - -Two fixes roughly doubled the achievable rate. Measured with -`--cycle-count 50 --bus-speed 400`: - -| TMRC | Nominal | Result here | Worst DRDY bracket | -|---|---|---|---| -| 0x92 | 600 Hz | abort | 1.86 ms of 1.67 ms (112%) | -| 0x93 | 300 Hz | abort | 3.45 ms of 3.33 ms (104%) | -| 0x94 | 150 Hz | borderline | 5.65 / 5.95 / 6.53 / 7.19 ms of 6.67 ms | -| 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.** They were measured in a -> resource-limited VM with virtualized USB, shared cores and background load -> (2.76 at the time of the one 150 Hz abort). Bracket tails there include -> hypervisor scheduling and USB passthrough latency that bare metal does not -> pay. The *shape* of the result should carry over — the limit is host round-trip -> latency, not I²C bandwidth or the sensor — but the numbers should be -> re-measured on the machine that will actually run the capture. The -> "% of margin used" line that `logger.py` prints after every run is the number -> to watch. - -**The I²C bus is never the constraint.** Per sample it carries a DRDY poll plus -a results read — 150 bit-times — against a host cost 2–3× larger. `logger.py` -prints both at the start of every run: - -| | 400 kHz | 100 kHz | -|---|---|---| -| 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 | — | - -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 -is always achievable. Every miss observed has been a host stall. - -The ceiling was therefore two host-side costs: - -**1. Two USB round trips per register read.** 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 in fact accepts a **repeated START**, so `write_read()` -in [ch347.py](ch347.py) does it in one. Measured at 400 kHz: - -| | separate | combined | -|---|---|---| -| `data_ready()` | 0.652 ms | **0.350 ms** | -| `read_raw()` | 0.838 ms | **0.548 ms** | - -**2. The GIL, not the OS scheduler.** Even after that, `logger.py` stalled ~10 ms -at 150 Hz while the *same sampling loop with no writer thread* peaked at 3.1 ms -(p99.9 = 2.9 ms, zero threshold crossings in 2079 samples). The sampler releases -the GIL inside each USB transfer and must re-acquire it, waiting up to Python's -default **5 ms** switch interval while the writer holds it. Setting -`sys.setswitchinterval(0.0005)` cut the worst bracket from 9.96 ms to 5.65 ms -and made 150 Hz pass. - -That is worth remembering generally: with a latency-sensitive loop and a helper -thread, the default GIL switch interval *is* the jitter floor. - -### Rate, resolution and the bus - -Four settings interact. Three of them are worth understanding together, and the -fourth is independent of the rest. - -**Two mechanisms set the rate, and the slower one wins** (§5.2.1): - -| | sets | granularity | -|---|---|---| -| cycle count | how *long* a measurement takes: `3 x (cc/90,000 + 68.7 us)` | continuous | -| TMRC | how *often* one is started | factor-of-two steps | - -If TMRC asks for something the cycle count cannot deliver, the cycle count wins -and TMRC has no effect at all — measured 73.85 Hz at TMRC 0x92 against 73.86 Hz -at 0x94 for the same cycle count. If TMRC asks for something slower, TMRC wins -and **the sensor idles**, which is the case to avoid. - -**Idling costs sensitivity.** Duty is integration time over period, and noise -after filtering scales as `1/sqrt(duty)`. Reaching ~130 Hz two ways: - -| | rate | cc | duty | ASD | -|---|---|---|---|---| -| 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. -So: **hold TMRC fast and let the cycle count set the rate.** That is what the -defaults do, and `--rate` derives both for you. - -**Cycle count also sets resolution**, which is the real trade. It buys rate and -LSB in opposite directions, and the dither margin follows the LSB: - -| cc | rate | nT/LSB | dither | Nyquist | -|---|---|---|---|---| -| 50 | 534 Hz | 50.37 | 0.58 | 267 Hz | -| 228 | 128 Hz | 11.74 | 1.17 | 64 Hz | -| 931 | 32 Hz | 2.91 | 2.34 | 16 Hz | - -Below ~0.46 Hz the 16-bit cycle-count register runs out and TMRC has to set the -cadence; below cycle count 30 the manual warns of quantisation (§5.1). Both -bounds are enforced by `rm3100.plan()`. - -**Bus speed is independent of the rate.** It does not appear in either rate -mechanism. What it sets is *latency* — how long a read takes, hence how tightly -a DRDY event 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 and 750 kHz is the default everywhere. - -### Recommended configuration - -**Cycle count 100 at 750 kHz**, which is what the defaults do. That runs the -sensor at its own ceiling of ~282 Hz, with TMRC held fast so the cycle count -governs and the duty cycle stays at 94%. - -```bash -./.venv/bin/python logger.py --duration 60 # cc=100, 750 kHz, ~282 Hz -``` - -It sits deliberately between the two things that pull in opposite directions: - -- **Against cycle count 200** it costs **1.5%** in post-filter noise — 94% duty - against 97%, and nothing else — 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**, so seeing it beats 1.5% of noise. At - 282 Hz both mains and its second harmonic sit in band and can be notched; at - 145 Hz only the fundamental does. -- **Against cycle count 50** it gives up half the spectrum (141 Hz of Nyquist - against 267 Hz) and buys **36% more dither margin**, 0.79 LSB against 0.58, - plus a 1.9× coarser rate demand on the host. cc=50 is the right choice when - something above 141 Hz needs identifying; it is not the right default, - because its dither margin is the one assumption in the whole configuration - that has not yet been checked against a real quiet-site capture. - -750 kHz is not optional at this cycle count. The period is 3.539 ms, and the -irreducible I2C traffic per sample is: - -| bus | traffic | share of period | -|---|---|---| -| 100 kHz | 1.500 ms | **42%** | -| 400 kHz | 0.375 ms | 11% | -| 750 kHz | 0.200 ms | **6%** | - -750 kHz is also the fastest measured, not merely the least bus traffic. Host -cost per sample, three runs each at cycle count 228: - -| bus | bus time | host cost | implied USB overhead | -|---|---|---|---| -| 100 kHz | 1.500 ms | 2.03–2.19 ms | ~0.63 ms | -| 400 kHz | 0.375 ms | 0.88–0.95 ms | ~0.55 ms | -| **750 kHz** | 0.200 ms | **0.78–0.80 ms** | ~0.59 ms | - -Host cost is `bus time + ~0.6 ms` of fixed USB round-trip latency, so a faster -bus helps but cannot go below that floor. An earlier note here claimed 400 kHz -benchmarked faster than 750; that rested on a single short run whose 4% -difference was inside the run-to-run spread, and does not survive repetition. -400 kHz remains a sane fallback if 750 proves marginal on long wiring. - -**Fall back to cycle count 200 if dither is still thin.** The assumption to check -is that cc=100 stays dithered: its intrinsic noise is 0.79 LSB against a -simulated failure point below ~0.2 LSB, so there is roughly 4× of margin at the -spec noise figure. That margin shrinks if the real sensor is quieter than spec, -which is the entire point of the quiet-site experiment. -`characterize.py` prints `sd/LSB` per axis, which settles it from the first -capture — below ~0.35 and averaging stops recovering sub-LSB resolution, at -which point 200 (1.10 LSB) buys the margin back at the cost of bandwidth. - -| cycle count | rate | Nyquist | nT/LSB | dither | duty | | -|---|---|---|---|---|---|---| -| 50 | 534 Hz | 267 Hz | 50.37 | 0.58 | 89.0% | spectrum, thin dither | -| **100** | **282 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 | - -`logger.py` warns on stderr for anything below cycle count 50 -(`rm3100.RECOMMENDED_MIN_CYCLE_COUNT`), along with five other configuration -traps — see [Configuration warnings](#configuration-warnings). - -### Configuration warnings - -`rm3100.plan()` resolves a configuration and `logger.py` prints its derivation -before every run, then checks it against six known traps. Each one is a *silent* -failure: the capture completes, the numbers look plausible, and the defect is -only visible 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 | -| TMRC governs and the sensor idles | TMRC-governed and 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 figure above it is derived from the rate | -| Nyquist below 60 Hz | rate < 120 Hz | mains folds onto signal and no later filter undoes it | -| cycle count past 400 | `cc > MAX_SPEC_CYCLE_COUNT` | Table 3-1 ends there, so the printed gain and noise are extrapolated, not specified | -| bus over half the period | traffic / period > 50% | names the speed that would fit | - -The aliasing one is the easiest to walk into: `--rate 32` yields a clean-looking -2.91 nT/LSB capture with 60 Hz mains sitting at 4.02 Hz, indistinguishable from -signal. Sampling fast and decimating afterwards gives the same noise floor with -the line still visible. - -### Rate model, and how TMRC is chosen - -Two things set the rate, and **the slower one wins** (§5.2.1). TMRC requests a -rate; the cycle count imposes a ceiling. Measured at TMRC 0x94 (150 Hz nominal, -138.8 Hz actual), varying only cc: - -| cc | measured | governed by | duty | -|---|---|---|---| -| 50 | 138.53 Hz | TMRC | 25% | -| 200 | 138.84 Hz | boundary | 99% | -| 400 | 73.83 Hz | cycle count | 100% | -| 800 | 37.09 Hz | cycle count | 100% | - -**When the cycle count governs, TMRC is irrelevant** — cc=400 gives 73.85 Hz at -TMRC 0x92 and 73.86 Hz at 0x94. So the simplest operating model is to set TMRC -fast (`0x92`) and use **cycle count as the single rate knob**: it is a continuous -integer where TMRC is coarse factor-of-two steps, and a cc-governed rate is -always ~100% duty, which sidesteps the quantisation trap entirely. - -TMRC is then only worth choosing when you want a rate *slower* than the cycle -count allows — to save power by letting the chip idle, or to hit a specific -cadence. - -The rate follows from a fitted model rather than Table 3-1's three data points -(`rm3100.sample_period` / `cycle_count_for_rate`): - -``` -per-axis time = cycle_count / 90,000 + 68.7 us -``` - -The divisor is the **specified** value, not a fitted one: Table 3-1 gives a -180 kHz circuit oscillation and §4.1 measures each count in *both* bias -directions, so one count costs two oscillations. Rates measured here across -cycle counts 229–29,769 agree with that to within **0.8%**, which confirms the -spec figure rather than improving on it — the residual is this unit's oscillator -inside ordinary component tolerance, and another part would sit elsewhere. - -The overhead has no specified value and must be measured. It is why the naive -`rate × cc` "constant" is not constant, drifting from 84,429 at cc=100 to -89,191 at cc=1200. - -**Treat predicted rates as ±2%.** The manual quotes no oscillator tolerance, so -anything needing the real number measures it — `logger.py` calibrates the period -against the host clock before recording, and that is what every capture reports. - -`sweep.py` measures this end to end. Every target from 1 to 128 Hz lands within -**0.1%**: - -``` - target measured err cc LSB/uT nT/count noise bus use - 1H 1.00H +0.0% 29769 10929.7 0.091 1.21n 0.0% - 8H 8.01H +0.1% 3718 1366.4 0.732 3.41n 0.3% - 32H 32.00H -0.0% 927 341.8 2.926 6.83n 1.2% - 128H 127.99H -0.0% 229 85.6 11.687 13.75n 4.8% - 256H 249.23H -2.6% 113 43.0 23.265 19.57n 9.3% 2 lost -``` - -256 Hz is where this host gives out; everything below it is exact. Note the -noise column past cycle count ~400 is extrapolation — the manual calls that its -useful upper limit and gives no data beyond. - -### Choosing cycle count - -Run each cycle count at *its own* fastest clean rate, not a fixed one. What then -matters is **duty** — the fraction of wall time the sensor is actually -integrating — because after decimating to a common bandwidth, noise scales as -1/√duty. Two ceilings compete: the sensor's (`~84,333/cc/3` Hz) and the host's -(~145 Hz here). - -| cc | sensor max | TMRC | actual | duty | noise | LSB | Nyquist | -|---|---|---|---|---|---|---|---| -| 50 | 562 Hz | 0x94 | 138.8 Hz | **25%** | **2.01×** | 50.37 nT | 69.4 Hz | -| 100 | 281 Hz | 0x94 | 138.8 Hz | 49% | 1.42× | 26.17 nT | 69.4 Hz | -| **200** | 141 Hz | **0x94** | **138.8 Hz** | **99%** | 1.01× | 13.35 nT | 69.4 Hz | -| 250 | 112 Hz | 0x95 | 69.4 Hz | 62% | 1.27× | 10.72 nT | 34.7 Hz | -| **400** | 70 Hz | **0x95** | **69.4 Hz** | **99%** | 1.01× | 6.74 nT | 34.7 Hz | -| 600 | 47 Hz | 0x96 | 34.2 Hz | 73% | 1.17× | 4.51 nT | 17.1 Hz | -| **800** | 35 Hz | **0x96** | **34.2 Hz** | **97%** | 1.01× | 3.39 nT | 17.1 Hz | - -**Because TMRC is quantised, cycle count should be chosen so the sensor ceiling -sits just *above* a TMRC step, not just below it.** Land just below and the rate -halves while integration time does not, wasting ~40% of the duty — cc=250 and -cc=600 are exactly that trap. - -The natural operating points are `cc ≈ 203 / 405 / 821`, i.e. **200 / 400 / 800** -paired with TMRC `0x94 / 0x95 / 0x96`. Both extremes verified on hardware: - -``` -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 -choice between them is purely **bandwidth versus resolution**, not noise. Use -cc=200/0x94 when you need 60 Hz inside the band; cc=800/0x96 when resolution -matters more than bandwidth. - -**Low cycle counts are dominated, not fast.** cc=50 reaches the same 138.8 Hz as -cc=200 — both are host-limited, not sensor-limited — but idles 75% of the time -and quantises 3.8× more coarsely. Its extra sensor speed is capability the host -cannot collect. - -### 60 Hz: resolved rather than aliased - -At 150 Hz, Nyquist is 69.3 Hz and mains lands **in band** — the strongest line -sits at **59.923 Hz**, observed directly. At 75 Hz it folds to 9.304 Hz, matching -the 9.233 Hz prediction, with only 4 nT/√Hz left at 60 Hz because that frequency -is no longer sampled. - -This settles the mains hypothesis by observation rather than inference, and gives -a way to *measure* interference before deciding how to reject it. - -### 60 s capture at maximum rate (office rig, historical) - -`capture_60s.csv` / `capture_60s.png` — 28,864 samples, 481 Hz, cycle count 50, -400 kHz bus. The sensor was nudged by hand at t ≈ 41 s, so statistics are split -around that to keep the noise figures honest: - -| Axis | sd, 0–40 s | sd, 47–60 s | step across the event | -|---|---|---|---| -| X | 84 nT | 85 nT | +280 nT | -| Y | 388 nT | 376 nT | +217 nT | -| Z | 187 nT | 198 nT | +232 nT | -| \|B\| | 369 nT | 358 nT | +313 nT | - -Three things worth noting: - -- **The t ≈ 41 s step was the sensor being physically moved**, confirmed at the - time — not an electrical artefact. Note that 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 field gradient samples a - different local field and changes the magnitude. Near a PC the gradients are - steep (§4.2.2: field falls off as 1/distance³), so a few centimetres is ample - to produce the +313 nT seen here. Treat |B| changes as "not a pure rotation", - nothing more. -- **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 asymmetry 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 Table 3-1 quotes for cycle count 50; Y is ~13×. -- **Sample timing is jittery**: intervals average 2.08 ms (481 Hz) with sd - 1.09 ms and a 19.6 ms worst case. Fine for logging, but the irregular spacing - makes this data unsuitable for spectral analysis without resampling. +The suite fakes the hardware at three seams: `FakeUsbDevice` (libusb, so +`ch347.py`'s framing is itself under test), `FakeBus` (an I²C master with a +register map), and `FakeSensor` (a scripted DRDY timeline against a fake clock, +so a host stall of an exact size can be injected). See `tests/conftest.py`. + +The miss-counting logic gets the most attention, since it is the part whose +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 and venv setup +- [x] udev rule, venv and dependency setup - [x] CH347 vendor protocol — verified byte-for-byte against the kernel driver -- [x] I2C bus scan — sensor found at 0x23 -- [x] Sensor identified — REVID `0x22` -- [x] Register write/read verified — cycle counts written and read back +- [x] I²C bus scan, sensor identified at 0x23, REVID `0x22` - [x] Analog section verified — BIST reports `XOK=YOK=ZOK=1` -- [x] **Live measurements logging to CSV at ~35 Hz** -- [x] Exactly-once sampling with abort-on-miss, and a chip-grid time base +- [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 - [ ] Re-test the LDO with the sensor **clamped**, to separate gain from movement -- [ ] Reliability testing over extended runs -- [ ] **Deal with 60 Hz mains aliasing into the band** — the dominant term -- [ ] Standoff distance, then a repeat characterisation +- [ ] Deal with 60 Hz coupling at source — the dominant term +- [ ] Reliability testing over extended runs on bare metal -> **Captures written before the header format are unreadable** and must be -> re-recorded — `capture.py` says so explicitly rather than guessing. The -> figures quoted under *Observed performance* came from those older files and -> stand as historical measurements only. +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. A proper Python API is planned; these modules are the initial-communication milestone, deliberately kept simple. diff --git a/capture.py b/capture.py index c7591a1..ce29ec4 100644 --- a/capture.py +++ b/capture.py @@ -57,14 +57,22 @@ class Capture: self.x, self.y, self.z = (counts[a] * self.lsb_nt for a in "xyz") self.total = np.sqrt(self.x**2 + self.y**2 + self.z**2) - # Least-squares fit of host clock against grid coordinate. numpy's - # polyfit is centred internally, so the large epoch offset is harmless. - slope, intercept = np.polyfit(sample_index, system_time, 1) + # 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. + 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. self.rate_error = self.dt_nominal / self.dt_true - 1.0 - self.residuals = system_time - (intercept + slope * sample_index) + 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 @@ -213,8 +221,10 @@ 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 aborts on a missed interval, so a gap here means the file was - # damaged or hand-edited rather than merely cut short. + # 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. gaps = np.diff(sample_index) if np.any(gaps != 1): bad = int(sample_index[np.argmax(gaps != 1)]) diff --git a/ch347.py b/ch347.py index 3604435..67a223a 100644 --- a/ch347.py +++ b/ch347.py @@ -46,6 +46,17 @@ 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. +SPEEDS = { + 20: SPEED_20KHZ, + 100: SPEED_100KHZ, + 400: SPEED_400KHZ, + 750: SPEED_750KHZ, +} +DEFAULT_SPEED_KHZ = 750 + # The count field is 6 bits, and a write also spends one byte on the address. MAX_XFER = 0x3F @@ -60,7 +71,7 @@ class CH347I2C: followed by the payload. """ - def __init__(self, speed=SPEED_750KHZ): + def __init__(self, speed=SPEEDS[DEFAULT_SPEED_KHZ]): self._dev = usb.core.find(idVendor=VENDOR_ID, idProduct=PRODUCT_ID) if self._dev is None: raise IOError( @@ -125,7 +136,7 @@ class CH347I2C: def set_speed(self, speed): """Select the I2C clock rate (one of the SPEED_* constants).""" - if speed not in (SPEED_20KHZ, SPEED_100KHZ, SPEED_400KHZ, SPEED_750KHZ): + if speed not in SPEEDS.values(): raise ValueError(f"Invalid I2C speed {speed}, expected 0-3") self._xfer([CMD_STREAM, CMD_SET | speed, CMD_END], 0) diff --git a/characterize.py b/characterize.py index 1388a3e..4211b4f 100644 --- a/characterize.py +++ b/characterize.py @@ -15,22 +15,41 @@ Produces a four-panel figure and a text summary: Read latency host-side diagnostics only -- the measurement grid itself is uniform regardless. +Two noise figures are printed, and the difference between them matters. `sd` is +the plain standard deviation, which on a drifting capture measures the drift: +one run here reads 167.9 nT of sd on an axis whose actual noise is 12.1 nT. +`white sd` is sd(diff)/sqrt(2), a first difference that rejects anything slower +than the sample rate, and is what Table 3-1's 208/sqrt(cycle_count) should be +compared against. + +`--trim` drops equal time from both ends. A capture usually starts while the +sensor is still settling and ends while it is being handled for the next run, +and neither belongs in a noise figure. Nothing is dropped unless asked. + The sample grid is exact, not assumed: the sensor samples on its own internal -schedule (manual sections 5.7.2, 5.8.2), logger.py aborts rather than skip an -interval, 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. +schedule (manual sections 5.7.2, 5.8.2), logger.py writes a placeholder rather +than skip a grid point, and capture.py refuses a non-contiguous sample index. +Frequencies use the rate measured against the host clock, not the nominal table +value, which is ~6% out on this unit. + +Rows lost during capture carry no data and are linearly interpolated on load, so +a capture with many of them will read as artificially smooth at high frequency. +capture.py's summary reports how many, and characterize.py prints it above. """ import argparse +import math import sys +from collections import namedtuple import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt +import matplotlib.ticker as mticker import numpy as np import capture +import rm3100 SURFACE = "#fcfcfb" TEXT_PRIMARY = "#0b0b0b" @@ -50,6 +69,48 @@ SERIES = AXES + [TOTAL] # Table 3-1: "Noise Density @ Max. Single-Axis Sample Rate". SPEC_ASD_NT = 1.2 +# Every panel is drawn on a fixed scale so two runs can be laid side by side and +# compared by eye. Autoscaling defeats that completely: a quieter capture simply +# redraws its own axis and looks identical to a noisy one. The ranges are wide +# enough for any configuration this driver offers -- cycle counts 50 to 800, so +# rates from ~2 to ~530 Hz and LSBs from 3 to 50 nT -- and are deliberately whole +# decades so the tick labels are powers of ten. +ASD_HZ = (1e-2, 1e3) # 0.05 Hz is the lowest bin a 20 s segment gives +ASD_NT = (1e-1, 1e2) # floor ~1 nT/rtHz, lines to ~15 +ALLAN_TAU_S = (1e-3, 1e3) # one sample period to a 1000 s average +ALLAN_NT = (1e-1, 1e2) +RESIDUAL_NT = (-150.0, 150.0) # +/-6 LSB at cycle count 100 +RESIDUAL_DENSITY = (0.0, 0.045) +LATENCY_MS = (0.0, 20.0) # covers 1.9 ms at cc=50 to 13.7 ms at cc=400 +LATENCY_COUNT = (0.8, 1e5) + +# Lower edge of the band every broadband figure is quoted over. Below a few Hz +# the spectrum is drift, not noise floor, and it varies far more between runs +# than the floor does -- so including it would compare environments rather than +# sensors. +BAND_LO_HZ = 3.0 +# Upper edge, as a fraction of Nyquist. This has to stay clear of the corner of +# any filter a compared series has been through, or the comparison scores one +# path partway down a rolloff and reads it as a quieter sensor. The binding one +# is decimation: DECIMATE_CUTOFF_FRACTION is where that filter turns over, so +# BAND_NYQUIST_FRACTION / 2 is held below it with margin to spare. +BAND_NYQUIST_FRACTION = 0.8 + +# Welch segment duration. Frequency resolution is 1/this, so 20 s of segment +# resolves 0.05 Hz -- fine against a band starting at 3 Hz -- and a 5-minute +# capture still yields ~30 overlapping segments to average. +SEGMENT_SECONDS = 20.0 + +# A sample-locked tone is reported at this many standard errors. Under the null +# the statistic is Rayleigh(1), so 5 is a one-in-270,000 bin and stays clear of +# the ~20 fractions tested. Real artefacts here reach 30. +LINE_SIGMA = 5.0 + + +def band_for(fs): + """The (lo, hi) analysis band for a capture sampled at fs.""" + return BAND_LO_HZ, BAND_NYQUIST_FRACTION * fs / 2 + def welch_asd(v, fs, nperseg=4096): """Amplitude spectral density in units/sqrt(Hz) via Welch's method.""" @@ -94,6 +155,267 @@ def allan_deviation(v, fs, points=40): return np.array(taus), np.array(devs) +def white_sd(v): + """Per-sample noise with drift removed, in the units of v. + + The first difference of white noise has twice its variance, so dividing by + sqrt(2) recovers the original sd -- while any component slower than the + sample rate is differenced away. On a capture that drifts, this and the + plain sd measure entirely different things, and only this one is the noise. + + It is not free: differencing has a +6 dB/octave response, so a narrowband + line near Nyquist is weighted about twice as heavily as the broadband floor. + Use band_stats for a figure that is flat across the band. + """ + if len(v) < 2: + return 0.0 + return float(np.diff(v).std() / np.sqrt(2)) + + +def segment_length(fs, seconds=SEGMENT_SECONDS): + """Welch segment length in samples for a fixed segment *duration*. + + Holding the duration fixed rather than the sample count is what makes two + captures at different rates comparable. A fixed nperseg gives them different + frequency resolutions and, worse, different numbers of segments to average: + the median of a Welch estimate is biased low by an amount that depends on + that count, so a slower capture would read as a quieter one by a few percent + for no physical reason at all. + """ + return max(64, int(round(fs * seconds))) + + +def band_stats(v, fs, band=None, nperseg=None): + """Broadband figures over a stated band: (median ASD, RMS, peak, peak Hz). + + The median is the floor -- robust to the handful of bins a line occupies -- + while the RMS integrates everything in the band including those lines, so + the two together say how much of the band is line rather than floor. + """ + lo, hi = band if band is not None else band_for(fs) + freqs, asd = welch_asd(v, fs, nperseg=nperseg or segment_length(fs)) + keep = (freqs >= lo) & (freqs <= hi) + if not keep.any(): + nan = float("nan") + return nan, nan, nan, nan + f, a = freqs[keep], asd[keep] + return (float(np.median(a)), float(np.sqrt(np.trapezoid(a ** 2, f))), + float(a.max()), float(f[np.argmax(a)])) + + +def fir_lowpass(cutoff, ntaps): + """Linear-phase low-pass FIR: a sinc truncated by a Blackman window. + + `cutoff` is in cycles per sample (0 to 0.5). Written out rather than taken + from scipy, which is not a dependency here and would be a large one for one + filter. Blackman buys a -74 dB stopband for a transition band of about + 5.5/ntaps, which is the trade that matters: the whole point of the filter is + that what it stops does not fold back. + + Normalised to unit gain at DC, so decimating does not rescale the field. + """ + if not 0 < cutoff < 0.5: + raise ValueError(f"cutoff {cutoff} must be in (0, 0.5) cycles/sample") + ntaps = int(ntaps) | 1 # odd, so the delay is a whole sample + if ntaps < 3: + raise ValueError(f"ntaps {ntaps} is too short to filter anything") + n = np.arange(ntaps) - (ntaps - 1) / 2 + h = 2 * cutoff * np.sinc(2 * cutoff * n) * np.blackman(ntaps) + return h / h.sum() + + +# Where the anti-alias filter turns over, as a fraction of the decimated rate. +# It must sit above BAND_NYQUIST_FRACTION / 2 so that nothing compared is +# scored on the rolloff, and below 0.5 so that nothing folds back. The +# transition band left over is (0.5 - this)/k, and Blackman gives about +# 5.5/ntaps, so the tap count has to grow with k: 128*k+1 leaves roughly 2x +# margin on that and is odd by construction. +DECIMATE_CUTOFF_FRACTION = 0.45 +assert BAND_NYQUIST_FRACTION / 2 < DECIMATE_CUTOFF_FRACTION < 0.5 + + +def decimate(v, k, method="fir"): + """Downsample by an integer k, low-passing first so nothing folds in. + + Two filters, because they answer different questions: + + boxcar average k consecutive samples. This is what the sensor itself + does over a cycle count, so it is the right comparison for "would + a slower cycle count have given me this?" -- and it shares the + boxcar's poor stopband, which is why the sensor aliases. + fir a proper anti-alias filter. This is what you would actually do in + post, and it is the one that removes out-of-band interference + instead of folding it. + + Returns the decimated series; the new rate is fs/k. + """ + k = int(k) + if k < 1: + raise ValueError(f"decimation factor {k} must be at least 1") + if k == 1: + return np.asarray(v, dtype=float) + v = np.asarray(v, dtype=float) + if method == "boxcar": + n = len(v) // k * k + if n < k: + raise ValueError(f"{len(v)} samples cannot be decimated by {k}") + return v[:n].reshape(-1, k).mean(axis=1) + if method != "fir": + raise ValueError(f"unknown decimation method {method!r}") + h = fir_lowpass(DECIMATE_CUTOFF_FRACTION / k, 128 * k + 1) + if len(v) <= 2 * len(h): + raise ValueError( + f"{len(v)} samples is too few for a {len(h)}-tap filter; " + f"decimating by {k} needs at least {2 * len(h) + 1}") + # Discard a full filter length at each end rather than half: 'same' pads + # with zeros, so the taper reaches ntaps//2 in and a half-length trim would + # leave the tail of it in the data. + return np.convolve(v, h, mode="same")[len(h):-len(h)][::k] + + +SampleLine = namedtuple("SampleLine", "numerator period amplitude sigma") + + +def sample_locked_lines(v, max_period=8, threshold=LINE_SIGMA): + """Tones sitting at a simple rational fraction of the sample rate. + + Returns SampleLine(numerator, period, amplitude, sigma) for each surviving + fraction numerator/period of fs, strongest first. `amplitude` is the peak + amplitude of the tone in the units of v. + + What this finds that a spectrum cannot distinguish is *what a line is locked + to*. A line at a fixed frequency lands on a different fraction of fs when + the rate changes; one that stays at fs/4 across captures whose rates differ + is locked to the sampling, so it is an artefact of the measurement rather + than a field. Comparing two captures settles it -- see compare.py. + + Evaluating the transform at exactly j/period rather than reading a spectrum + puts the whole capture behind one number, which is what makes a 1 nT tone + detectable at 30 sigma under 20 nT of broadband noise. + + Only fractions in lowest terms are tested, so each frequency is reported + once. Without that, a period-4 pattern would also be reported at 6 and 8, + which are multiples of it and carry no additional information. + + A cubic trend is removed first, so drift cannot leak into a low fraction. + """ + v = np.asarray(v, dtype=float) + if len(v) < 4 * max_period: + return [] + index = np.arange(len(v)) + v = v - np.polyval(np.polyfit(index, v, 3), index) + sd = v.std() + if sd == 0: + return [] + found = [] + for period in range(2, max_period + 1): + for numerator in range(1, period // 2 + 1): + if math.gcd(numerator, period) != 1: + continue + if 2 * numerator == period: + # Nyquist has no phase to fit -- it is a real alternation, so + # one degree of freedom rather than two, and a tighter null. + amplitude = abs(float((v * (-1.0) ** index).mean())) + sigma = amplitude / (sd / np.sqrt(len(v))) + else: + phase = np.exp(-2j * np.pi * numerator / period * index) + amplitude = 2 * abs(complex((v * phase).mean())) + # Real and imaginary parts each have variance 2 sd^2 / len, so + # the magnitude is Rayleigh(1) in these units under the null: + # 5 sigma is a one-in-270,000 fluctuation. + sigma = amplitude / (sd * np.sqrt(2 / len(v))) + if sigma >= threshold: + found.append(SampleLine(numerator, period, float(amplitude), + float(sigma))) + return sorted(found, key=lambda line: -line.sigma) + + +def dither_check(v, lsb, longest=1024): + """Does averaging still recover resolution below one quantiser step? + + Returns [(n, sd, sd_in_lsb, ratio_to_ideal), ...] for block averages of n + samples. Ideal is sd(1)/sqrt(n), so a ratio near 1 means averaging is + buying everything it should. + + This is the check the whole cycle-count choice rests on. Section 5.1 warns + that quantisation "generally dictates working above a cycle count of ~30", + and the failure it warns about is specific: if the sensor's own noise is + small against the LSB, samples stop straddling the boundary, the quantiser + stops being dithered, and averaging stalls at a fixed fraction of an LSB no + matter how long you average. What that looks like here is a ratio that + climbs while the LSB column stops falling. + + A slow trend is removed first. Otherwise drift -- which averaging cannot + remove and is not supposed to -- is what stops the average shrinking, and + the answer would be about the site rather than the quantiser. + """ + v = np.asarray(v, dtype=float) + index = np.arange(len(v)) + v = v - np.polyval(np.polyfit(index, v, 5), index) + base = white_sd(v) + if base == 0: + return [] + out = [] + n = 1 + while n <= longest and len(v) // n >= 16: + usable = len(v) // n * n + means = v[:usable].reshape(-1, n).mean(axis=1) + sd = float(means.std()) + out.append((n, sd, sd / lsb, sd / (base / np.sqrt(n)))) + n *= 4 + return out + + +def decade_ticks(ax, which="both", minimum_decades=1.0): + """Label a log axis at powers of ten, with unlabelled minors between. + + Matplotlib falls back to labelling minor ticks when a log axis spans less + than a decade, which produces a scale reading 2x10^0, 3x10^0, 4x10^0 -- + dense, and not what a log axis is read as. Decades are the convention, so + the range is widened to hold at least one rather than the labelling being + changed to suit a narrow range. + """ + axes = {"x": [ax.xaxis], "y": [ax.yaxis], + "both": [ax.xaxis, ax.yaxis]}[which] + for axis in axes: + setter = ax.set_xlim if axis is ax.xaxis else ax.set_ylim + lo, hi = ax.get_xlim() if axis is ax.xaxis else ax.get_ylim() + if lo > 0 and hi > lo: + short = minimum_decades - math.log10(hi / lo) + if short > 0: + pad = 10 ** (short / 2) + setter(lo / pad, hi * pad) + axis.set_major_locator(mticker.LogLocator(base=10.0)) + axis.set_minor_locator( + mticker.LogLocator(base=10.0, subs=tuple(range(2, 10)))) + axis.set_minor_formatter(mticker.NullFormatter()) + + +def trimmed(cap, seconds): + """Drop `seconds` from both ends of a capture. Returns (capture, note). + + A capture typically opens while the sensor is still settling after power-up + and closes while it is being handled for whatever comes next, and neither + end is a noise measurement. Trimming is never silent -- the note says what + went -- and never automatic, because discarding data the caller did not ask + to discard is how a figure ends up quietly describing a subset. + + Refuses rather than trims when there would be too little left to mean + anything, since a short deliberate capture is a legitimate thing to analyse. + """ + if not seconds: + return cap, "" + if seconds < 0: + raise capture.CaptureError(f"--trim {seconds} must not be negative") + if cap.duration < 4 * seconds: + return cap, (f"not trimming: {cap.duration:.1f} s is under 4x the " + f"{seconds:g} s requested from each end") + kept = cap.restrict(seconds, cap.duration - seconds) + return kept, (f"trimmed {seconds:g} s from each end: " + f"{len(cap.sample_index):,} -> {len(kept.sample_index):,} " + f"samples, {cap.duration:.1f} -> {kept.duration:.1f} s") + + def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) @@ -103,12 +425,24 @@ def main(): help="ignore samples before this elapsed time (s)") ap.add_argument("--end", type=float, default=None, help="ignore samples after this elapsed time (s)") + ap.add_argument("--trim", type=float, default=0.0, metavar="SECONDS", + help="drop this many seconds from BOTH ends -- settling at " + "the start, handling at the end (default: %(default)s)") args = ap.parse_args() + # --trim and --start/--end both choose a window, and silently letting one + # win would make the reported figures depend on argument order. + if args.trim and (args.start or args.end is not None): + ap.error("--trim sets the window from both ends; it cannot be combined " + "with --start/--end, which set it explicitly") + try: cap = capture.load(args.csv) if args.start or args.end is not None: cap = cap.restrict(args.start or None, args.end) + whole = cap.duration + cap, trim_note = trimmed(cap, args.trim) + trim_applied = cap.duration < whole except (OSError, capture.CaptureError) as exc: sys.exit(str(exc)) @@ -128,9 +462,17 @@ def main(): # differs -- a gain change carries the noise with it, so absolute sd alone # will read as a noise difference when only the scale moved. field = data["total"].mean() + band = band_for(fs) + # Table 3-1's noise figure for this cycle count -- what `white sd` is the + # measurement of, so the two belong in the same table. + spec_nt = rm3100.expected_noise_nt(cap.cycle_count) - print(cap.summary() + "\n") - print("series sd sd/|B| p2p sd/LSB white-noise ASD median ASD") + print(cap.summary()) + if trim_note: + print(f" {trim_note}") + print() + print(f"{'series':6s} {'sd':>9s} {'white sd':>9s} {'dither':>7s} " + f"{'sd/|B|':>10s} {'p2p':>9s} {'vs spec':>8s} {'ASD':>7s}") fig, axs = plt.subplots(2, 2, figsize=(13.5, 9), dpi=150) fig.patch.set_facecolor(SURFACE) @@ -140,8 +482,9 @@ def main(): for key, label, color in SERIES: v = data[key] sd = v.std() - # A flat (white) spectrum of this sd would sit at this level. - implied = sd / np.sqrt(fs / 2) + # Drift-free: what Table 3-1 quotes, and what the dither margin is + # really made of. The plain sd above it is drift on a drifting capture. + wsd = white_sd(v) freqs, asd = welch_asd(v, fs) axs[0, 0].loglog(freqs, asd, color=color, linewidth=1.2, @@ -151,28 +494,85 @@ def main(): 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. + # a fat tail. Bin edges are locked to quarter-LSB steps: the residual + # of an integer-valued signal is a comb at one LSB by construction, and + # arbitrary bin edges beat against it into ragged spikes that look like + # structure. Aligned bins render the comb as what it is. resid = v - np.polyval(np.polyfit(t, v, 3), t) - axs[1, 0].hist(resid, bins=120, histtype="step", linewidth=1.4, + step = lsb_nt / 4 + edge = np.ceil(max(abs(resid.min()), abs(resid.max())) / step) * step + axs[1, 0].hist(resid, bins=np.arange(-edge, edge + step, step), + histtype="step", linewidth=1.4, color=color, label=label, density=True) - print(f"{label:6s} {sd:8.1f} {sd/field*1e6:7.1f}ppm {v.max()-v.min():9.1f} nT " - f"{sd/lsb_nt:7.2f} {implied:9.2f} nT/rtHz " - f"{np.median(asd):9.2f} nT/rtHz") + median_asd, _, _, _ = band_stats(v, fs, band) + print(f"{label:6s} {sd:9.1f} {wsd:9.2f} {wsd/lsb_nt:7.2f} " + f"{sd/field*1e6:9.1f}p {v.max()-v.min():9.1f} " + f"{(wsd/spec_nt - 1) * 100:+7.1f}% {median_asd:7.2f}") + + print(f" sd, white sd, p2p in nT; dither = white sd in LSB " + f"(1 LSB = {lsb_nt:.2f} nT);") + print(f" sd/|B| in ppm of {field:,.0f} nT; vs spec against Table 3-1's " + f"{spec_nt:.1f} nT at cycle count {cap.cycle_count};") + print(f" ASD = median nT/rtHz over {band[0]:g}-{band[1]:.1f} Hz.") + + # The assumption the cycle-count choice rests on, checked rather than + # asserted: with dither under one LSB, does averaging still get below it? + print("\ndither: averaging against the quantiser") + for key, label, _ in SERIES: + steps = dither_check(data[key], lsb_nt) + if not steps: + continue + longest = steps[-1] + print(f" {label:6s} {steps[0][2]:.2f} LSB per sample -> " + f"{longest[2]:.3f} LSB after {longest[0]:,} " + f"({longest[0] / fs:.1f} s), {longest[3]:.1f}x the ideal 1/sqrt(n)") + print(" Sub-LSB resolution is being recovered wherever the LSB column " + "keeps falling.\n A ratio climbing while it stalls is the " + "quantiser losing its dither (section 5.1);\n a ratio climbing " + "while it falls is just drift, which averaging cannot remove.") + + # Structure locked to the sample index rather than to a frequency. Reported + # per axis because it is not isotropic -- one axis here carries 5x another. + lines = [(label, sample_locked_lines(data[key])) + for key, label, _ in SERIES] + if any(found for _, found in lines): + print("\nlines at simple fractions of the sample rate:") + for label, found in lines: + for line in found: + fraction = line.numerator / line.period + print(f" {label:6s} {line.numerator}/{line.period} of fs " + f"= {fraction * fs:8.3f} Hz {line.amplitude:6.2f} nT " + f"amplitude {line.sigma:5.1f} sigma") + print(" These may be sample-locked -- moving with the rate rather than " + "sitting at a\n fixed frequency -- which would make them " + "artefacts of measuring, not field.\n compare.py decides it, by " + "checking whether they hold the same fraction of fs\n in a " + "capture taken at a different rate.") a = axs[0, 0] + # The band the quoted median comes from, so the number and the picture + # cannot drift apart. Drawn under the traces, hence the low zorder. + a.axvspan(band[0], band[1], color=GRID, alpha=0.55, zorder=0, linewidth=0) a.axhline(SPEC_ASD_NT, color=REFERENCE, linestyle="--", linewidth=1.2) a.annotate(f"Table 3-1 spec {SPEC_ASD_NT} nT/√Hz", xy=(freqs[1], SPEC_ASD_NT), xytext=(0, 5), textcoords="offset points", color=REFERENCE, fontsize=9) a.set_title("Amplitude spectral density", loc="left", color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) + a.annotate(f"shaded: {band[0]:g}–{band[1]:.0f} Hz, the quoted band", + xy=(0.02, 0.04), xycoords="axes fraction", + color=TEXT_SECONDARY, fontsize=9) a.set_xlabel("frequency (Hz)"); a.set_ylabel("nT/√Hz") + a.set_xlim(*ASD_HZ); a.set_ylim(*ASD_NT) + decade_ticks(a) a = axs[0, 1] a.set_title("Allan deviation", loc="left", color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) a.set_xlabel("averaging time τ (s)"); a.set_ylabel("σ (nT)") + a.set_xlim(*ALLAN_TAU_S); a.set_ylim(*ALLAN_NT) + decade_ticks(a) a.annotate("slope −½ = white noise; upturn = drift", xy=(0.02, 0.04), xycoords="axes fraction", color=TEXT_SECONDARY, fontsize=9) @@ -180,10 +580,20 @@ def main(): a = axs[1, 0] a.set_title("Residual distribution (cubic trend removed)", loc="left", color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) + a.annotate(f"bins are {lsb_nt / 4:.2f} nT, a quarter LSB.\n" + f"Integer counts make this a comb at {lsb_nt:.1f} nT however\n" + f"well dithered; how sharp it looks only tracks how far\n" + f"the axis drifted. Whether averaging still beats the\n" + f"quantiser is the dither figure printed above.", + xy=(0.02, 0.97), xycoords="axes fraction", va="top", + color=TEXT_SECONDARY, fontsize=9) a.set_xlabel("nT"); a.set_ylabel("density") + a.set_xlim(*RESIDUAL_NT); a.set_ylim(*RESIDUAL_DENSITY) a = axs[1, 1] - a.hist(latency * 1000, bins=120, color=REFERENCE) + # Fixed bin edges as well as fixed limits: comparing two histograms whose + # bins are different widths compares the binning as much as the data. + a.hist(latency * 1000, bins=np.linspace(*LATENCY_MS, 121), color=REFERENCE) a.set_yscale("log") a.axvline(cap.dt_true * 1000, color=TEXT_PRIMARY, linestyle="--", linewidth=1.2) a.annotate(f"grid {cap.dt_true * 1000:.2f} ms", xy=(cap.dt_true * 1000, 1), @@ -195,6 +605,8 @@ def main(): xy=(0.98, 0.94), xycoords="axes fraction", ha="right", va="top", color=TEXT_SECONDARY, fontsize=9) a.set_xlabel("interval between reads (ms)"); a.set_ylabel("count") + a.set_xlim(*LATENCY_MS); a.set_ylim(*LATENCY_COUNT) + decade_ticks(a, which="y") for ax in axs.flat: ax.grid(True, which="both", color=GRID, linewidth=0.7) @@ -218,7 +630,9 @@ 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"Sample grid is exact.", + f"Sample grid is exact." + + (f" {args.trim:g} s trimmed from each end." if trim_applied + else ""), color=TEXT_SECONDARY, fontsize=10, ha="center") fig.tight_layout(rect=[0, 0, 1, 0.935]) diff --git a/compare.py b/compare.py index 8e11892..9db0cac 100644 --- a/compare.py +++ b/compare.py @@ -2,141 +2,966 @@ """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 # group by header note + ./.venv/bin/python compare.py --group note *.csv # label from header + ./.venv/bin/python compare.py LDO/cc100=a.csv 3V3/cc100=b.csv ... -Built around one hard lesson: a naive comparison of absolute noise is wrong when -the measured scale differs between runs. A gain change carries the noise with it, -so a run that reads 6% larger also reads ~6% noisier while being physically -identical. Everything here is therefore reported **fractionally**, in ppm of the -field magnitude. +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 +are treated as the same configuration. `--supply LDO=3.0` then attaches a rail +voltage to every capture in that condition. -The second trap is attributing a scale change to the variable under test when -the sensor simply moved. Two diagnostics separate them: +Three lessons are built in, each learned by getting it wrong first. - per-axis ratio spread a pure gain change scales X, Y and Z identically, - so the spread is ~0. Anything larger means the sensor - moved. - rotation angle the angle between mean field directions. ~0 deg means - the sensor held still. +**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. -Both must be small before a magnitude difference can be blamed on gain. Note -that |B| is preserved under rotation but *not* under translation through a field -gradient, so a moved sensor can change magnitude on its own. +**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. + +**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: + + 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. """ import argparse import itertools +import math +import os import sys from collections import defaultdict +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt import numpy as np import capture import characterize as ch +import rm3100 # A pure gain change scales every axis by the same factor. Allow a little for # noise on the means before calling it movement. 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. +RATIO_AXIS_MIN_FRACTION = 0.05 # of |B| -def summarise(cap): +# Two captures whose rates differ by less than this are treated as the same +# rate, so a line holding a fixed fraction of fs across them proves nothing. +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. +CONDITION_COLORS = ["#2a78d6", "#eb6834", "#1baf7a"] +# Cycle count is the secondary encoding, carried by line style and marker, so +# identity never rests on colour alone. +VARIANT_STYLES = ["-", "--", ":", "-."] +VARIANT_MARKERS = ["o", "s", "^", "D"] + + +# -------------------------------------------------------------------------- +# Arguments +# -------------------------------------------------------------------------- + +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 + 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. + """ + label, separator, path = argument.partition("=") + if not separator: + return None, argument + if not label: + raise ValueError(f"{argument!r} has an empty label before '='") + if not os.path.exists(path) and os.path.exists(argument): + return None, argument + return label, path + + +def parse_supply(argument): + """Parse `CONDITION=VOLTS` for --supply.""" + condition, separator, volts = argument.partition("=") + if not separator or not condition: + raise argparse.ArgumentTypeError( + f"--supply {argument!r} must be CONDITION=VOLTS, e.g. LDO=3.0") + try: + value = float(volts) + except ValueError: + raise argparse.ArgumentTypeError( + f"--supply {argument!r}: {volts!r} is not a number") + if not 0 < value < 100: + raise argparse.ArgumentTypeError( + f"--supply {argument!r}: {value} V is not a plausible rail") + return condition, value + + +def parse_band(text): + """Parse `LO,HI` for --band.""" + parts = text.split(",") + if len(parts) != 2: + raise argparse.ArgumentTypeError(f"--band {text!r} must be LO,HI in Hz") + try: + lo, hi = (float(p) for p in parts) + except ValueError: + raise argparse.ArgumentTypeError(f"--band {text!r} must be two numbers") + if not 0 < lo < hi: + raise argparse.ArgumentTypeError( + f"--band {text!r} needs 0 < LO < HI") + return lo, hi + + +def label_for(cap, path, group_by, explicit): + if explicit: + return explicit + if group_by == "note": + return cap.meta.get("note", "(no note)") + return os.path.basename(path) + + +def condition_of(label): + """The part of a label before the first '/': the thing under test.""" + return label.split(CONDITION_SEP, 1)[0] + + +def variant_of(label): + """The part after the first '/', or '' when the label carries none.""" + _, separator, variant = label.partition(CONDITION_SEP) + return variant if separator else "" + + +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. + """ + seen = [] + for label, _, _ in records: + condition = condition_of(label) + if condition not in seen: + seen.append(condition) + return seen + + +# -------------------------------------------------------------------------- +# Per-capture summary +# -------------------------------------------------------------------------- + +def summarise(cap, band): d = dict(cap.axes()) mean = np.array([d["x"].mean(), d["y"].mean(), d["z"].mean()]) field = float(np.linalg.norm(mean)) fs = cap.true_rate_hz - out = {"cap": cap, "mean": mean, "field": field, "axes": d, "fs": fs} + out = {"cap": cap, "mean": mean, "field": field, "axes": d, "fs": fs, + "cycle_count": cap.cycle_count} for key, _, _ in ch.SERIES: v = d[key] - freqs, asd = ch.welch_asd(v, fs, nperseg=min(4096, len(v) // 4 * 2)) - band = freqs > min(3.0, fs / 8) + median_asd, band_rms, peak, peak_hz = ch.band_stats(v, fs, band) out[key] = { "sd": v.std(), + "white": ch.white_sd(v), "ppm": v.std() / field * 1e6, - "asd": float(np.median(asd[band])) if band.any() else float("nan"), - "peak": float(asd[band].max()) if band.any() else float("nan"), - "peak_hz": float(freqs[band][np.argmax(asd[band])]) if band.any() else float("nan"), + "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, + "band_rms_ppm": band_rms / field * 1e6, + "peak": peak, + "peak_hz": peak_hz, } taus, devs = ch.allan_deviation(d["total"], fs) i = int(np.argmin(devs)) out["allan"] = (float(devs[i]), float(taus[i])) + out["allan_curve"] = (taus, devs) + out["lines"] = ch.sample_locked_lines(d["total"]) return out -def label_for(cap, path, group_by): - if group_by == "note": - return cap.meta.get("note", "(no note)") +# -------------------------------------------------------------------------- +# Cross-capture analysis +# -------------------------------------------------------------------------- + +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. + """ + usable = [i for i in range(3) + if abs(a["mean"][i]) >= RATIO_AXIS_MIN_FRACTION * a["field"] + and abs(b["mean"][i]) >= RATIO_AXIS_MIN_FRACTION * b["field"]] + if len(usable) < 2: + return None + kept = ratio[usable] + return float(np.ptp(kept) / np.abs(np.mean(kept))) + + +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. + """ + ma, mb = a["mean"], b["mean"] + scale = float(ma @ mb / (ma @ ma)) + residual = float(np.linalg.norm(mb - scale * ma) / np.linalg.norm(mb)) + ua, ub = ma / a["field"], mb / b["field"] + angle = float(np.degrees(np.arccos(np.clip(ua @ ub, -1, 1)))) + return scale, residual, angle + + +def fit_rate_model(caps): + """Solve rm3100's timing model for a set of captures: (counts/s, overhead). + + The model is `period = AXES * (cycle_count / C + h)`, so period/AXES is + 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 + h together by one factor, leaving the *shape* of the model alone. + """ + counts = np.array([float(c.cycle_count) for c in caps]) + if len(np.unique(counts)) < 2: + raise ValueError("need at least two distinct cycle counts to fit " + "both the count rate and the per-axis overhead") + periods = np.array([c.dt_true for c in caps]) / rm3100.AXES + slope, intercept = np.polyfit(counts, periods, 1) + if slope <= 0: + raise ValueError("fitted count rate is not positive; the captures do " + "not follow the timing model") + return 1.0 / float(slope), float(intercept) + + +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. + + 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. + """ + v_fast = dict(fast["cap"].axes())["total"] + v_slow = dict(slow["cap"].axes())["total"] + fs_fast, fs_slow = fast["fs"], slow["fs"] + rows = [("raw fast", v_fast, fs_fast, fast["field"])] + for method in ("boxcar", "fir"): + rows.append((f"{method} /{k}", ch.decimate(v_fast, k, method), + fs_fast / k, fast["field"])) + rows.append(("native slow", v_slow, fs_slow, slow["field"])) + out = [] + for name, v, fs, field in rows: + median_asd, band_rms, _, _ = ch.band_stats(v, fs, band) + 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)}) + return out + + +def decimation_pairs(entries): + """[(condition, fast, slow, k)] for every pair a decimation can bridge. + + `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. + """ + by_condition = defaultdict(list) + for label, item in entries: + by_condition[condition_of(label)].append((label, item)) + pairs = [] + for condition, group in by_condition.items(): + for (la, a), (lb, b) in itertools.combinations(group, 2): + fast, slow, lf, ls = ((a, b, la, lb) if a["fs"] > b["fs"] + else (b, a, lb, la)) + if not fast["cycle_count"] or slow["cycle_count"] % fast["cycle_count"]: + continue + k = slow["cycle_count"] // fast["cycle_count"] + if k >= 2: + pairs.append((condition, (lf, fast), (ls, slow), k)) + return pairs + + +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. + + Runs before the summaries exist, so it takes captures rather than records. + """ + entries = [(label, {"fs": cap.true_rate_hz, + "cycle_count": cap.cycle_count}) + for label, cap in loaded] + rates = [item["fs"] for _, item in entries] + rates += [fast["fs"] / k for _, (_, fast), _, k in decimation_pairs(entries)] + return rates + + +def locked_line_verdict(records): + """Which rational-fraction lines hold their fraction across rates. + + 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. + """ + seen = defaultdict(list) + for label, _, s in records: + for line in s["lines"]: + seen[(line.numerator, line.period)].append((label, s["fs"], line)) + verdicts = [] + for fraction, hits in sorted(seen.items(), + key=lambda kv: -max(h[2].sigma for h in kv[1])): + rates = [fs for _, fs, _ in hits] + distinct = (len(rates) > 1 + and (max(rates) / min(rates) - 1) > RATE_DISTINCT) + verdicts.append((fraction, hits, distinct)) + return verdicts + + +# -------------------------------------------------------------------------- +# Output sections +# -------------------------------------------------------------------------- + +def print_overview(records, band, supplies): + print(f"{'capture':<20} {'cc':>5} {'fs Hz':>9} {'|B| nT':>10} {'sd':>7} " + f"{'white':>7} {'ASD':>6} {'Allan':>6} | {'white':>7} {'ASD':>7}") + for label, _, s in records: + t = s["total"] + print(f"{label[:20]:<20} {s['cycle_count']:5d} {s['fs']:9.3f} " + 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}") + print(f" Left of the bar, absolute: sd/white/Allan in nT, ASD in nT/rtHz " + 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. " + "Allan = best sigma by averaging.") + if supplies: + rails = " ".join(f"{c} {v:g} V" for c, v in sorted(supplies.items())) + print(f" rails: {rails}") + + +def print_rate_model(records): + by_condition = defaultdict(list) + for label, _, s in records: + by_condition[condition_of(label)].append(s["cap"]) + fits = {} + for condition in conditions_in_order(records): + try: + fits[condition] = fit_rate_model(by_condition[condition]) + except ValueError: + continue + if not fits: + return + print("\n=== timing model: period = 3 x (cycle_count / C + overhead) ===") + for condition, (count_rate, overhead) in fits.items(): + print(f" {condition:<12} C {count_rate:10,.0f} counts/s " + f"overhead {overhead * 1e6:6.2f} us") + print(f" {'rm3100.py':<12} C {rm3100.COUNTS_PER_SECOND:10,.0f} counts/s " + f"overhead {rm3100.AXIS_OVERHEAD_S * 1e6:6.2f} us (the model's own)") + if len(fits) < 2: + return + print("\n Everything the chip times comes off one oscillator, so a " + "condition that only") + print(" changes its frequency scales C up and the overhead down by the " + "same factor:") + for (ca, fa), (cb, fb) in itertools.combinations(fits.items(), 2): + count_ratio = fb[0] / fa[0] + overhead_ratio = fb[1] / fa[1] + print(f" {ca} -> {cb} C x{count_ratio:.5f} " + f"({(count_ratio - 1) * 100:+.2f}%) " + f"overhead x{overhead_ratio:.5f} " + f"({(overhead_ratio - 1) * 100:+.2f}%)") + print(f" pure clock scaling would put the overhead at " + f"x{1 / count_ratio:.5f}; it is " + f"{(overhead_ratio * count_ratio - 1) * 100:+.2f}% off that.") + + +def print_locked_lines(records): + verdicts = locked_line_verdict(records) + if not verdicts: + return + print("\n=== lines at simple fractions of the sample rate (|B|) ===") + for (numerator, period), hits, distinct in verdicts: + print(f" {numerator}/{period} of fs:") + for label, fs, line in hits: + print(f" {label[:24]:<24} {line.numerator/line.period*fs:8.3f} Hz" + f" {line.amplitude:6.2f} nT {line.sigma:5.1f} sigma") + 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.") + 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.") + else: + print(" -> seen in one capture only, so nothing here says " + "whether it is locked to\n the sampling or to a " + "frequency.") + + +def print_decimation(records, band): + pairs = decimation_pairs([(label, s) for label, _, s in records]) + if not pairs: + return + print("\n=== does filtering and decimating the fast capture match the " + "slow one? ===") + for condition, (lf, fast), (ls, slow), k in pairs: + print(f"\n {condition}: {lf} decimated by {k} vs {ls}") + try: + rows = decimation_paths(fast, slow, k, band) + except ValueError as exc: + print(f" cannot decimate: {exc}") + continue + print(f" {'path':<14} {'fs Hz':>8} {'ASD':>7} {'ASD ppm':>8} " + f"{'white sd':>9} sample-locked lines") + for row in rows: + lines = (", ".join(f"{l.numerator}/{l.period} {l.amplitude:.2f} nT" + for l in row["lines"]) or "none") + print(f" {row['name']:<14} {row['fs']:8.2f} {row['asd']:7.2f} " + f"{row['asd_ppm']:8.1f} {row['white']:9.2f} {lines}") + 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. + 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.") + 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.") + + +# 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. +CORRECTIONS = [("raw", "raw"), + ("rate", "x measured fs"), + ("supply", "x rail volts"), + ("both", "x fs x volts")] + + +def corrected_fields(s, supplies, reference_volts): + """|B| under each candidate correction, keyed by name; None if unavailable. + + rate the chip's own clock is measured every run, so a gain that tracks + the clock is correctable from the capture alone, with nothing + external to know. + 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. + """ + raw = s["field"] + by_rate = raw * (1.0 + s["cap"].rate_error) + volts = supplies.get(condition_of(s["label"])) + ratio = (volts / reference_volts) if volts else None + return {"raw": raw, + "rate": by_rate, + "supply": raw * ratio if ratio else None, + "both": by_rate * ratio if ratio else None} + + +def print_contrasts(records, supplies): + """Condition contrasts at matched cycle count -- the designed comparison.""" + conditions = conditions_in_order(records) + if len(conditions) < 2: + return + reference_volts = min(supplies.values()) if supplies else None + by_cc = defaultdict(dict) + for label, _, s in records: + s["label"] = label + by_cc[s["cycle_count"]][condition_of(label)] = s + + print("\n=== |B| by condition, at matched cycle count ===") + print("The comparison the interleaving was for: same cycle count, " + "so only the condition differs.") + contrasts = defaultdict(list) + exponents = defaultdict(list) + for cc in sorted(by_cc): + present = by_cc[cc] + if len(present) < 2: + continue + for ca, cb in itertools.combinations(conditions, 2): + if ca not in present or cb not in present: + continue + a, b = present[ca], present[cb] + left = corrected_fields(a, supplies, reference_volts) + right = corrected_fields(b, supplies, reference_volts) + print(f"\n cycle count {cc}: {ca} -> {cb}") + for key, name in CORRECTIONS: + if left[key] is None or right[key] is None: + continue + print(f" {name:<14} {left[key]:10,.1f} -> " + f"{right[key]:10,.1f} nT " + 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. + va, vb = supplies.get(ca), supplies.get(cb) + if va and vb and va != vb: + dv = math.log(vb / va) + exponents[(ca, cb)].append( + (cc, math.log(right["raw"] / left["raw"]) / dv, + math.log(b["fs"] / a["fs"]) / dv)) + + 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.") + + for (ca, cb), points in contrasts.items(): + if len(points) < 2: + continue + values = [v for _, v in points] + print(f"\n {ca} -> {cb} across {len(points)} cycle counts: " + 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.") + + # A condition measured twice bounds everything this experiment cannot + # separate -- handling, gradient, drift -- and is the only honest error bar. + print("\n reproducibility within one condition (the error bar on all of " + "the above):") + any_repeat = False + for condition in conditions: + runs = [(s["cycle_count"], s["field"]) for l, _, s in records + if condition_of(l) == condition] + if len(runs) < 2: + continue + any_repeat = True + fields = [f for _, f in runs] + print(f" {condition:<12} n={len(runs)} |B| " + + ", ".join(f"{f:,.0f}" for f in fields) + + 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.") + else: + 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.") + for (la, _, a), (lb, _, b) in itertools.combinations(records, 2): + ratio = b["mean"] / a["mean"] + spread = axis_ratio_spread(a, b, ratio) + scale, residual, angle = gain_and_movement(a, b) + print(f"\n{la[:24]} -> {lb[:24]}") + print(f" |B| {a['field']:,.0f} -> {b['field']:,.0f} nT " + f"{b['field'] - a['field']:+,.1f} nT " + f"({(b['field']/a['field'] - 1) * 100:+.2f}%)") + print(f" fractional noise {a['total']['ppm']:.1f} -> " + f"{b['total']['ppm']:.1f} ppm " + f"({(b['total']['ppm']/a['total']['ppm'] - 1) * 100:+.1f}%)") + print(f" best-fit scale {scale:.5f} gain residual " + f"{residual * 100:5.2f}% rotation {angle:.3f} deg") + if spread is not None: + print(f" per-axis ratios X {ratio[0]:.5f} Y {ratio[1]:.5f} " + f"Z {ratio[2]:.5f} spread {spread * 100:.2f}%") + else: + print(f" per-axis ratios fewer than two axes carry " + f"{RATIO_AXIS_MIN_FRACTION:.0%} of |B|, so this " + "test abstains") + 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.") + else: + print(" -> the sensor held still, so the |B| difference is a " + "scale change.") + + +# -------------------------------------------------------------------------- +# 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. + """ + edges = np.geomspace(freqs[0], freqs[-1], bins + 1) + index = np.digitize(freqs, edges) - 1 + centres, values = [], [] + for b in range(bins): + inside = index == b + if inside.any(): + centres.append(np.sqrt(edges[b] * edges[b + 1])) + values.append(np.median(asd[inside])) + return np.array(centres), np.array(values) + + +def clip_to_band(ax, series, band, headroom=3.0, floor=0.35): + """Set y-limits from what is inside the band, not from filter stopbands. + + An anti-alias filter falls away by six decades past its corner, and on a log + axis that decides the scale for everything else. Only the band being + compared should. + """ + inside = [a[(f >= band[0]) & (f <= band[1])] for f, a in series] + inside = [a for a in inside if len(a)] + if not inside: + return + lo = min(float(a.min()) for a in inside) + hi = max(float(a.max()) for a in inside) + ax.set_ylim(floor * lo, headroom * hi) + + +def make_figure(records, band, supplies, path): + conditions = conditions_in_order(records) + variants = sorted({variant_of(l) for l, _, _ in records}) + color = {c: CONDITION_COLORS[i % len(CONDITION_COLORS)] + for i, c in enumerate(conditions)} + dash = {v: VARIANT_STYLES[i % len(VARIANT_STYLES)] + for i, v in enumerate(variants)} + marker = {v: VARIANT_MARKERS[i % len(VARIANT_MARKERS)] + for i, v in enumerate(variants)} + + 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) + + a = axs[0, 0] + a.axvspan(band[0], band[1], color=ch.GRID, alpha=0.55, zorder=0, + linewidth=0) + drawn = [] + for label, _, s in records: + freqs, asd = binned_asd(*ch.welch_asd(dict(s["cap"].axes())["total"], s["fs"], + nperseg=ch.segment_length(s["fs"]))) + # 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. + asd = asd / s["field"] * 1e6 + drawn.append((freqs, asd)) + a.loglog(freqs, asd, color=color[condition_of(label)], + linestyle=dash[variant_of(label)], linewidth=1.4, alpha=0.9, + label=label) + clip_to_band(a, drawn, band) + ch.decade_ticks(a) + title(a, "Spectral density of |B|, as a fraction of |B|") + a.set_xlabel("frequency (Hz)"); a.set_ylabel("ppm of |B| per √Hz") + 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", + xy=(0.02, 0.05), xycoords="axes fraction", + color=ch.TEXT_SECONDARY, fontsize=9) + a.legend(frameon=False, fontsize=9, labelcolor=ch.TEXT_SECONDARY, + loc="upper right") + + a = axs[0, 1] + plotted = _plot_decimation(a, records, band, color, dash) + title(a, "Fast capture, filtered and decimated") + a.set_xlabel("frequency (Hz)"); a.set_ylabel("ppm of |B| per √Hz") + if plotted: + ch.decade_ticks(a) + a.legend(frameon=False, fontsize=9, labelcolor=ch.TEXT_SECONDARY, + loc="upper right") + else: + a.annotate("no pair of captures differs by an integer cycle-count " + "factor,\nso there is nothing to decimate onto", + xy=(0.5, 0.5), xycoords="axes fraction", ha="center", + color=ch.TEXT_SECONDARY, fontsize=10) + + a = axs[1, 0] + for label, _, s in records: + taus, devs = s["allan_curve"] + a.loglog(taus, devs / s["field"] * 1e6, + color=color[condition_of(label)], + linestyle=dash[variant_of(label)], linewidth=1.6, label=label) + ch.decade_ticks(a) + title(a, "Allan deviation of |B|, as a fraction of |B|") + a.set_xlabel("averaging time τ (s)"); a.set_ylabel("σ (ppm of |B|)") + a.annotate("slope −½ = white noise; upturn = drift", + xy=(0.02, 0.04), xycoords="axes fraction", + color=ch.TEXT_SECONDARY, fontsize=9) + a.legend(frameon=False, fontsize=9, labelcolor=ch.TEXT_SECONDARY, + loc="upper right") + + _plot_corrections(axs[1, 1], records, supplies, color, marker) + + fig.suptitle("RM3100 capture comparison", color=ch.TEXT_PRIMARY, + fontsize=15, fontweight="bold", y=0.985) + fig.text(0.5, 0.945, + f"{len(records)} captures; colour is the condition, line style " + f"the cycle count. Broadband figures over " + f"{band[0]:g}–{band[1]:.1f} Hz, inside every capture's Nyquist.", + color=ch.TEXT_SECONDARY, fontsize=10, ha="center") + fig.tight_layout(rect=[0, 0, 1, 0.935]) + fig.savefig(path, facecolor=ch.SURFACE) return path +def _plot_decimation(ax, records, band, color, dash): + """Decimated-fast against natively-slow, per condition. True if anything.""" + ax.axvspan(band[0], band[1], color=ch.GRID, alpha=0.55, zorder=0, + linewidth=0) + plotted = False + drawn = [] + for condition, (_, fast), (_, slow), k in decimation_pairs( + [(label, s) for label, _, s in records]): + try: + v = ch.decimate(dict(fast["cap"].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. + asd = asd / fast["field"] * 1e6 + drawn.append((freqs, asd)) + ax.loglog(freqs, asd, color=hue, linestyle="--", linewidth=1.5, + label=f"{condition} cc{fast['cycle_count']} ÷{k}") + freqs, asd = binned_asd(*ch.welch_asd( + dict(slow["cap"].axes())["total"], slow["fs"], + nperseg=ch.segment_length(slow["fs"]))) + asd = asd / slow["field"] * 1e6 + drawn.append((freqs, asd)) + ax.loglog(freqs, asd, color=hue, linestyle="-", linewidth=1.5, + alpha=0.7, + label=f"{condition} cc{slow['cycle_count']} native") + plotted = True + if plotted: + clip_to_band(ax, drawn, band) + ax.annotate("dashed: fast capture decimated. solid: natively slow.", + xy=(0.02, 0.05), xycoords="axes fraction", + color=ch.TEXT_SECONDARY, fontsize=9) + return plotted + + +def _plot_corrections(ax, records, supplies, color, marker): + """|B| under each correction, so the spread each one leaves is visible.""" + reference_volts = min(supplies.values()) if supplies else None + for label, _, s in records: + s["label"] = label + corrected = [(label, corrected_fields(s, supplies, reference_volts)) + for label, _, s in records] + rows = [(key, name) for key, name in CORRECTIONS + if all(c[key] is not None for _, c in corrected)] + + for row, (key, name) in enumerate(rows): + values = [] + for label, correction in corrected: + values.append(correction[key]) + ax.plot(correction[key], -row, marker=marker[variant_of(label)], + color=color[condition_of(label)], markersize=9, + markeredgecolor=ch.SURFACE, markeredgewidth=2, + linestyle="none", zorder=3) + ax.plot([min(values), max(values)], [-row, -row], + color=ch.GRID, linewidth=6, solid_capstyle="round", zorder=1) + ax.annotate(f"spread {(max(values)/min(values) - 1) * 100:.1f}%", + xy=(max(values), -row), xytext=(10, 0), + textcoords="offset points", va="center", + color=ch.TEXT_SECONDARY, fontsize=9) + 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") + ax.set_xlabel("|B| (nT)") + ax.annotate("a correction that is right collapses the spread", + xy=(0.02, 0.06), xycoords="axes fraction", + color=ch.TEXT_SECONDARY, fontsize=9) + # Identity never rests on colour alone: the marker repeats the cycle count. + handles = [plt.Line2D([], [], color=color[c], marker="o", markersize=9, + linestyle="none", label=c) + for c in conditions_in_order(records)] + handles += [plt.Line2D([], [], color=ch.TEXT_SECONDARY, + marker=marker[v], markersize=8, linestyle="none", + label=v or "(one variant)") + for v in sorted({variant_of(l) for l, _, _ in records})] + ax.legend(handles=handles, frameon=False, fontsize=9, ncol=2, + labelcolor=ch.TEXT_SECONDARY, loc="upper left") + + +# -------------------------------------------------------------------------- + def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) - ap.add_argument("csv", nargs="+") + ap.add_argument("csv", nargs="+", metavar="[LABEL=]CSV") ap.add_argument("--group", choices=["note", "file"], default="file", - help="group captures by header note, averaging repeats " + help="where unlabelled captures take their name from " "(default: %(default)s)") + ap.add_argument("--trim", type=float, default=0.0, metavar="SECONDS", + help="drop this many seconds from BOTH ends of every " + "capture (default: %(default)s)") + ap.add_argument("--band", type=parse_band, default=None, metavar="LO,HI", + help=f"broadband comparison band in Hz; the default is " + f"{ch.BAND_LO_HZ:g} Hz to " + f"{ch.BAND_NYQUIST_FRACTION * 100:g}%% of the Nyquist " + f"of the slowest path compared, decimation included") + ap.add_argument("--supply", type=parse_supply, action="append", + default=[], metavar="CONDITION=VOLTS", + help="rail voltage for a condition, enabling the " + "ratiometric correction; repeatable") + ap.add_argument("-o", "--output", default=None, metavar="PNG", + help="write a comparison figure here") args = ap.parse_args() + supplies = dict(args.supply) - recs = [] - for path in args.csv: + loaded = [] + for argument in args.csv: try: + explicit, path = split_label(argument) cap = capture.load(path) - except (OSError, capture.CaptureError) as exc: - print(f"skipping {path}: {exc}", file=sys.stderr) + cap, note = ch.trimmed(cap, args.trim) + except (OSError, ValueError, capture.CaptureError) as exc: + print(f"skipping {argument}: {exc}", file=sys.stderr) continue - recs.append((label_for(cap, path, args.group), path, summarise(cap))) - if len(recs) < 1: + loaded.append((label_for(cap, path, args.group, explicit), path, cap, + note)) + if not loaded: sys.exit("nothing to compare") - print(f"{'capture':<26} {'|B| nT':>10} {'sd nT':>8} {'sd ppm':>9} " - f"{'ASD':>7} {'peak':>8} {'@Hz':>7} {'Allan':>7}") - for label, path, s in recs: - t = s["total"] - print(f"{label[:26]:<26} {s['field']:10,.0f} {t['sd']:8.1f} {t['ppm']:9.1f} " - f"{t['asd']:7.2f} {t['peak']:8.1f} {t['peak_hz']:7.3f} " - f"{s['allan'][0]:7.2f}") - print(" ASD/peak in nT/rtHz; Allan = best sigma by averaging") + unknown = supplies.keys() - {condition_of(l) for l, _, _, _ in loaded} + if unknown: + print(f"warning: --supply names no capture: {', '.join(sorted(unknown))}", + file=sys.stderr) - # Averaged per group, which is the number to compare when runs are repeated. - if args.group == "note": - groups = defaultdict(list) - for label, _, s in recs: - groups[label].append(s) - if len(groups) > 1: - print("\n=== group means (fractional -- the comparable figure) ===") - for label, ss in groups.items(): - ppm = [s["total"]["ppm"] for s in ss] - print(f"{label[:26]:<26} n={len(ss)} " - f"|B| {np.mean([s['field'] for s in ss]):9,.0f} nT " - f"sd {np.mean(ppm):7.1f} ppm" - + (f" +/- {np.std(ppm):.1f}" if len(ss) > 1 else "")) + # 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. + 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) - if len(recs) < 2: - return 0 + records = [(label, path, summarise(cap, band)) for label, path, cap, _ in loaded] - print("\n=== pairwise: did the sensor hold still? ===") - print("A magnitude difference only means gain if BOTH checks pass.") - for (la, _, a), (lb, _, b) in itertools.combinations(recs, 2): - ratio = b["mean"] / a["mean"] - spread = float(np.ptp(ratio) / np.abs(np.mean(ratio))) - ua, ub = a["mean"] / a["field"], b["mean"] / b["field"] - angle = float(np.degrees(np.arccos(np.clip(np.dot(ua, ub), -1, 1)))) - moved = spread > RATIO_SPREAD_OK or angle > ROTATION_OK_DEG - print(f"\n{la[:24]} -> {lb[:24]}") - print(f" |B| ratio {b['field']/a['field']:.5f} " - f"({(b['field']/a['field'] - 1) * 100:+.2f}%)") - print(f" fractional noise {a['total']['ppm']:.1f} -> {b['total']['ppm']:.1f} ppm " - f"({(b['total']['ppm']/a['total']['ppm'] - 1) * 100:+.1f}%)") - print(f" per-axis ratios X {ratio[0]:.5f} Y {ratio[1]:.5f} Z {ratio[2]:.5f}" - f" spread {spread * 100:.2f}%") - print(f" rotation {angle:.3f} deg") - if moved: - print(" -> SENSOR MOVED. The magnitude difference cannot be " - "attributed to gain;\n re-run with the sensor clamped.") - else: - print(f" -> held still. The {(b['field']/a['field'] - 1) * 100:+.2f}% " - "magnitude difference is a real gain change.") + if args.trim: + for label, _, _, note in loaded: + if note: + print(f"{label[:24]:<24} {note}") + print() + print_overview(records, band, supplies) + print_rate_model(records) + print_locked_lines(records) + print_decimation(records, band) + print_contrasts(records, supplies) + print_pairwise(records) + + if args.output: + print(f"\n-> {make_figure(records, band, supplies, args.output)}") return 0 diff --git a/diagnose-comms.py b/diagnose-comms.py index eefcb9f..56023ed 100644 --- a/diagnose-comms.py +++ b/diagnose-comms.py @@ -19,8 +19,7 @@ import rm3100 # the most forgiving self-test settings (manual Tables 5-6, 5-7). BIST_RUN = 0x8F -# Matches the CH347I2C default; named here only so the report can state it. -BUS_KHZ = 750 +BUS_KHZ = ch347.DEFAULT_SPEED_KHZ def check(label, ok, detail=""): diff --git a/logger.py b/logger.py index bf86f42..18b4d38 100755 --- a/logger.py +++ b/logger.py @@ -38,16 +38,14 @@ CSV_FIELDS = ["sample_index", "system_time_unix", "x_raw", "y_raw", "z_raw", # signal irreversibly, which is the one error post-processing cannot undo. MAINS_HZ = 60.0 -BUS_SPEEDS = { - 20: ch347.SPEED_20KHZ, - 100: ch347.SPEED_100KHZ, - 400: ch347.SPEED_400KHZ, - 750: ch347.SPEED_750KHZ, -} - 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. +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 -- @@ -125,7 +123,11 @@ def set_thread_nice(value): # 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. +# 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. QUEUE_MAX = 200_000 _SENTINEL = object() @@ -182,7 +184,8 @@ def parse_args(): help="CSV output path (default: rm3100_.csv)") p.add_argument("--address", type=lambda s: int(s, 0), default=None, help="I2C address, skipping the scan (default: autodetect)") - p.add_argument("--bus-speed", type=int, choices=[20, 100, 400, 750], default=750, + p.add_argument("--bus-speed", type=int, choices=sorted(ch347.SPEEDS), + default=ch347.DEFAULT_SPEED_KHZ, help="I2C bus speed in kHz. The default cycle count runs near " "282 Hz, where 100 kHz would spend 42%% of each period on " "the bus; 750 spends 6%% (default: %(default)s)") @@ -295,7 +298,7 @@ def print_plan(cfg, bus_speed, requested_rate, sampler_nice=None): "extrapolated, not specified.") if share > 0.5: - faster = [s for s in sorted(BUS_SPEEDS) if i2c_bus_time(s) / period < 0.25] + faster = [s for s in sorted(ch347.SPEEDS) if i2c_bus_time(s) / period < 0.25] fix = f"; {faster[0]} kHz would fit" if faster else "" warnings.append(f"the bus needs {share:.0%} of every period{fix}.") @@ -356,8 +359,6 @@ def start_sensor(sensor, cycle_count, tmrc): # Raises if DRC1 did not take; exactly-once sampling depends on it. sensor.configure() sensor.set_rate(tmrc) - - sensor.start_cmm() print("Continuous measurement mode started") return revid @@ -495,8 +496,12 @@ 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. + before = monotonic() read_raw() - read_cost += monotonic() - ready + read_cost += monotonic() - before times.append(ready) clear = ready n = len(times) @@ -518,13 +523,16 @@ def sample_loop(sensor, q, duration, dt_nominal, stats): 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. + + 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. """ 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. @@ -607,7 +615,8 @@ def sample_loop(sensor, q, duration, dt_nominal, stats): put((index, previous_wall + dt_est * (k + 1), None, None, None, flags)) except queue.Full: - stats["dropped"] += 1 + stats["truncated"] = True + return index += 1 sample_wall = (clear_wall + ready_wall) / 2.0 @@ -618,7 +627,8 @@ def sample_loop(sensor, q, duration, dt_nominal, stats): put((index, sample_wall, cx, cy, cz, WARN_AMBIGUOUS if uncertain and not lost else "")) except queue.Full: - stats["dropped"] += 1 + stats["truncated"] = True + return index += 1 previous_mono, previous_wall = ready_mono, sample_wall # The read just cleared DRDY, so this is a known-clear instant. @@ -659,7 +669,7 @@ def main(): args.cycle_count, args.tmrc = cfg.cycle_count, cfg.tmrc try: - bus = ch347.CH347I2C(BUS_SPEEDS[args.bus_speed]) + bus = ch347.CH347I2C(ch347.SPEEDS[args.bus_speed]) except IOError as exc: sys.exit(str(exc)) @@ -740,8 +750,8 @@ def main(): print(f"Logging to {path} -- Ctrl-C to stop\n") q = queue.Queue(maxsize=QUEUE_MAX) - stats = {"rows": 0, "dropped": 0, "missed": 0, "ambiguous": 0, - "max_bracket": 0.0} + stats = {"rows": 0, "missed": 0, "ambiguous": 0, "max_bracket": 0.0, + "truncated": False} writer = threading.Thread(target=writer_thread, args=(q, path, meta, stats), daemon=True) writer.start() @@ -752,7 +762,8 @@ def main(): pass finally: q.put(_SENTINEL) - writer.join(timeout=30) + writer.join(timeout=WRITER_DRAIN_TIMEOUT_S) + drained = not writer.is_alive() print() try: sensor.stop_cmm() @@ -760,9 +771,6 @@ def main(): except IOError as exc: print(f"WARNING: could not stop CMM: {exc}", file=sys.stderr) - if stats["dropped"]: - print(f"WARNING: dropped {stats['dropped']} samples -- the writer " - "could not keep up", 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. @@ -772,7 +780,30 @@ def main(): 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. + 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) + + 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. + 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) + 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 " @@ -781,6 +812,7 @@ def main(): "capture's\ntime axis for spectral work.", 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" @@ -788,7 +820,8 @@ def main(): "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) - return 1 + + return 1 if failed else 0 return 0 diff --git a/noise_floor_test.sh b/noise_floor_test.sh deleted file mode 100755 index af3e215..0000000 --- a/noise_floor_test.sh +++ /dev/null @@ -1,94 +0,0 @@ -#!/usr/bin/env bash -# -# Run a noise-floor capture and analyse it in one step. -# -# Wraps logger.py -- every logger flag is accepted and passed straight through -- -# then runs plot.py and characterize.py on whatever was captured. Works whether -# the run ends on --duration or on Ctrl-C. -# -# ./noise_floor_test.sh --duration 60 -# ./noise_floor_test.sh --duration 3600 --cycle-count 400 -# ./noise_floor_test.sh # until Ctrl-C -# -# Anything you pass explicitly wins; otherwise logger.py's defaults apply. - -set -uo pipefail # deliberately not -e: logger.py exits 1 on a miss, - # and a partial capture is still worth analysing. - -cd "$(dirname "$0")" - -if [[ ! -f .venv/bin/activate ]]; then - echo "No virtualenv found. Run ./setup.sh first." >&2 - exit 1 -fi -# shellcheck disable=SC1091 -source .venv/bin/activate - -args=("$@") - -has_flag() { - local needle=$1 a - for a in "${args[@]:-}"; do - [[ "$a" == "$needle" || "$a" == "$needle="* ]] && return 0 - done - return 1 -} - -# The CSV path has to be known up front so the analysis can find it; logger.py -# would otherwise pick a timestamped name we cannot predict without racing it. -output="" -for ((i = 0; i < ${#args[@]}; i++)); do - case "${args[i]}" in - --output=*) output="${args[i]#--output=}" ;; - --output|-o) output="${args[i+1]:-}" ;; - esac -done -if [[ -z "$output" ]]; then - output="noise_$(date +%Y%m%d_%H%M%S).csv" - args+=(--output "$output") -fi - -has_flag --bus-speed || args+=(--bus-speed 750) - -echo "=== capture -> $output ===" -# Ctrl-C reaches the whole foreground process group. Trapping it here stops bash -# killing the script, so logger.py can shut down cleanly and the analysis below -# still runs. A no-op trap (not 'ignore') keeps the child's own handling intact. -trap ':' INT -python logger.py "${args[@]}" -capture_rc=$? -trap - INT - -# A non-zero logger exit covers both "aborted mid-capture" and "never started" -# (no adapter, bad flag). Whether a usable file exists is what distinguishes -# them, so check that before claiming anything about the data. -if [[ ! -s "$output" ]]; then - echo - echo "Capture produced no data (exit $capture_rc); nothing to analyse." >&2 - exit 1 -fi -if [[ "$capture_rc" -ne 0 ]]; then - echo - echo "Capture ended early (exit $capture_rc). Analysing what was recorded --" \ - "it is valid up to that point." >&2 -fi - -echo -echo "=== plot ===" -python plot.py "$output" -plot_rc=$? - -echo -echo "=== characterize ===" -python characterize.py "$output" -char_rc=$? - -echo -echo "=== outputs ===" -base="${output%.csv}" -for f in "$output" "$base.png" "${base}_noise.png"; do - [[ -f "$f" ]] && printf ' %-40s %8s\n' "$f" "$(du -h "$f" | cut -f1)" -done - -# Non-zero if any stage failed, so this is usable from a scheduler. -[[ $capture_rc -eq 0 && $plot_rc -eq 0 && $char_rc -eq 0 ]] || exit 1 diff --git a/setup.sh b/setup.sh index 60e904f..04335c7 100755 --- a/setup.sh +++ b/setup.sh @@ -10,7 +10,25 @@ set -euo pipefail VID=1a86 PID=55db -RULE_FILE=/etc/udev/rules.d/99-ch347.rules + +# The 60- prefix is load-bearing, not a style choice. udev sorts every rules file +# from every directory into one lexicographic sequence (udev(7), "RULES FILES"), +# and the only thing that acts on TAG+="uaccess" is this line in systemd's own +# /usr/lib/udev/rules.d/73-seat-late.rules: +# +# TAG=="uaccess", ENV{MAJOR}!="", RUN{builtin}+="uaccess" +# +# A rule numbered above 73 adds the tag *after* the rule that would have honoured +# it has already run: the tag is set, nothing ever reads it, and no ACL is +# granted. That was silently the case at 99-, which is why Bazzite ended up with +# no access at all -- there uaccess is the only mechanism that can work, since +# plugdev does not exist and usermod cannot add the group (see section 4). +# systemd keeps its own uaccess rules in 70-uaccess.rules for this reason; 60 is +# the prefix other projects settled on after hitting the same bug. +RULE_FILE=/etc/udev/rules.d/60-ch347.rules +# Written by earlier versions of this script. Removed on sight: a stale copy +# matching the same device would keep applying its own (wrong) permissions. +LEGACY_RULE_FILES=(/etc/udev/rules.d/99-ch347.rules /etc/udev/rules.d/70-ch347.rules) cd "$(dirname "$0")" @@ -27,34 +45,65 @@ echo echo "== 2. Installing udev rule ==" # Access is granted two ways, because no single mechanism covers every distro: # -# TAG+="uaccess" systemd-logind puts an ACL on the node for whoever is -# logged in at the local seat. Needs no group and no logout, -# and is the only thing that works unaided on atomic Fedora +# TAG+="uaccess" systemd-logind puts an ACL on the node for whoever holds +# the local seat. Needs no group and no logout, and is the +# only thing that can work unaided on atomic Fedora # (Bazzite, Silverblue, Kinoite) -- see section 4. # GROUP=/MODE= the traditional fallback, for ssh sessions and seatless # systems where there is no local seat for uaccess to grant. # -# The group is picked from those that exist. plugdev is a Debian convention and -# is absent on Fedora; naming a group that does not exist makes udev log an -# error and leave the node owned by root, which is worse than omitting it. +# They go on SEPARATE LINES, which matters more than it looks. udev drops a +# whole rule line whose GROUP= cannot be resolved, and since systemd 258 that +# includes any group which exists but is not a *system* group ("Group 'x' is not +# a system group, ignoring"). With both directives on one line, an unusable +# group therefore takes the uaccess tag down with it -- losing the fallback +# costs the primary mechanism too. Split, each stands or falls alone. ACCESS_GROUP="" +ACCESS_GROUP_GID="" for g in plugdev dialout; do - if getent group "$g" >/dev/null; then ACCESS_GROUP=$g; break; fi + if gid=$(getent group "$g" | cut -d: -f3) && [[ -n "$gid" ]]; then + ACCESS_GROUP=$g + ACCESS_GROUP_GID=$gid + break + fi done -RULE="SUBSYSTEM==\"usb\", ATTRS{idVendor}==\"$VID\", ATTRS{idProduct}==\"$PID\", TAG+=\"uaccess\"" -if [[ -n "$ACCESS_GROUP" ]]; then - RULE="$RULE, GROUP=\"$ACCESS_GROUP\", MODE=\"0660\"" - echo "Access via uaccess, falling back to group '$ACCESS_GROUP'." -else +MATCH="SUBSYSTEM==\"usb\", ATTRS{idVendor}==\"$VID\", ATTRS{idProduct}==\"$PID\"" +RULE="# Installed by rm3100/setup.sh. Prefix must stay below 73 so that +# 73-seat-late.rules still sees the uaccess tag; see the comments in setup.sh. +$MATCH, TAG+=\"uaccess\"" + +if [[ -z "$ACCESS_GROUP" ]]; then + # The normal case on Fedora and its derivatives: plugdev was removed years + # ago and dialout is for serial ports. uaccess alone is the supported path. echo "No plugdev or dialout group here; relying on uaccess alone." +elif [[ "$ACCESS_GROUP_GID" -ge 1000 ]]; then + # systemd 258 refuses these outright, so writing the line would achieve + # nothing while looking like a fallback. + echo "Group '$ACCESS_GROUP' has GID $ACCESS_GROUP_GID, which systemd 258+" + echo "rejects as a non-system group; relying on uaccess alone." + ACCESS_GROUP="" +else + RULE="$RULE +$MATCH, GROUP=\"$ACCESS_GROUP\", MODE=\"0660\"" + echo "Access via uaccess, falling back to group '$ACCESS_GROUP'." fi -if [[ -f "$RULE_FILE" ]] && [[ "$(cat "$RULE_FILE")" == "$RULE" ]]; then +STALE=() +for old in "${LEGACY_RULE_FILES[@]}"; do + [[ -e "$old" ]] && STALE+=("$old") +done + +if [[ -f "$RULE_FILE" ]] && [[ "$(cat "$RULE_FILE")" == "$RULE" ]] && + [[ ${#STALE[@]} -eq 0 ]]; then echo "$RULE_FILE already up to date, skipping." else echo "Writing $RULE_FILE (needs sudo)..." printf '%s\n' "$RULE" | sudo tee "$RULE_FILE" >/dev/null + if [[ ${#STALE[@]} -gt 0 ]]; then + echo "Removing superseded ${STALE[*]}" + sudo rm -f "${STALE[@]}" + fi sudo udevadm control --reload-rules # --action=add, not the default 'change': the uaccess builtin that sets the # ACL only runs on add, so a change event would install the rule without @@ -68,7 +117,12 @@ else fi echo -echo "== 3. Creating virtualenv and installing pyusb ==" +echo "== 3. Creating virtualenv and installing dependencies ==" +# pyusb drives the adapter; numpy and matplotlib are needed by every analysis +# tool (capture/plot/characterize/compare), so installing only pyusb leaves +# noise_floor_test.sh failing on an ImportError halfway through a capture. +# pytest runs the suite in tests/, which needs no hardware. +DEPS=(pyusb numpy matplotlib pytest) if [[ ! -d .venv ]]; then python3 -m venv .venv echo "Created .venv" @@ -76,30 +130,70 @@ else echo ".venv already exists" fi ./.venv/bin/pip install --quiet --upgrade pip -./.venv/bin/pip install --quiet pyusb -echo "pyusb $(./.venv/bin/python -c 'import usb; print(usb.__version__)') installed" +./.venv/bin/pip install --quiet "${DEPS[@]}" +./.venv/bin/python - <<'EOF' +import importlib.metadata as md +for name in ("pyusb", "numpy", "matplotlib", "pytest"): + print(f" {name} {md.version(name)}") +EOF echo echo "== 4. Verifying device node permissions ==" # The node itself is the ground truth, so nothing above is warned about -# speculatively -- group membership only matters if uaccess did not already -# cover it, and that is visible here. +# speculatively. This also says *which* mechanism granted access, because the +# two fail in completely different ways and the fix differs accordingly. +USER_NAME=$(id -un) NODE=$(lsusb -d "$VID:$PID" | head -1 | sed -E 's|Bus ([0-9]+) Device ([0-9]+).*|/dev/bus/usb/\1/\2|') ls -l "$NODE" getfacl -p "$NODE" 2>/dev/null | grep -E "^user:[^:]+:" || true +# An ACL entry naming this user is the fingerprint of uaccess having run. Its +# absence is the symptom the 60- prefix exists to prevent. +has_uaccess_acl() { + getfacl -p "$NODE" 2>/dev/null | grep -qE "^user:$USER_NAME:[^:]*w" +} + +# Exactly what the uaccess builtin requires: an active session on a seat. Asked +# through the property API rather than by parsing `loginctl list-sessions`, +# whose columns move between systemd versions. +has_active_seat_session() { + command -v loginctl >/dev/null 2>&1 || return 1 + local sessions session + sessions=$(loginctl show-user "$USER_NAME" --property=Sessions --value 2>/dev/null) || return 1 + for session in $sessions; do + [[ "$(loginctl show-session "$session" --property=Active --value 2>/dev/null)" == "yes" ]] && + [[ -n "$(loginctl show-session "$session" --property=Seat --value 2>/dev/null)" ]] && + return 0 + done + return 1 +} + if [[ -w "$NODE" ]]; then - echo "OK: $NODE is writable by $(id -un)." + if has_uaccess_acl; then + echo "OK: $NODE is writable by $USER_NAME (uaccess ACL)." + else + echo "OK: $NODE is writable by $USER_NAME (group membership)." + fi else - USER_NAME=$(id -un) echo "WARNING: $NODE is not writable by $USER_NAME." >&2 echo >&2 - echo " uaccess grants nothing over ssh or without a local seat, so fall" >&2 - echo " back to the group. Try, in order:" >&2 + echo " 1. Unplug and replug the adapter. The rule is applied when the" >&2 + echo " device appears, so an already-connected node keeps whatever" >&2 + echo " permissions it was given before. This alone is usually enough." >&2 echo >&2 - echo " 1. Unplug and replug the adapter, so the rule applies to a fresh" >&2 - echo " node. This is enough on a local desktop session." >&2 + if has_active_seat_session; then + echo " 2. You do hold an active local seat, so uaccess should apply." >&2 + echo " If a replug does not fix it, the rules are not loaded:" >&2 + echo " sudo udevadm control --reload-rules && sudo udevadm trigger" >&2 + echo " Bazzite in particular is known not to reload /etc/udev/rules.d" >&2 + echo " when it switches to the final rootfs, so rules can sit inert" >&2 + echo " until that is run once after each boot (ublue-os/bazzite#2516)." >&2 + else + echo " 2. You have no active local seat session -- this looks like ssh" >&2 + echo " or a headless login. uaccess grants nothing there, by design," >&2 + echo " so the group fallback is the only route." >&2 + fi if [[ -n "$ACCESS_GROUP" ]] && ! id -nG | tr ' ' '\n' | grep -qx "$ACCESS_GROUP"; then echo >&2 @@ -109,21 +203,28 @@ else # group is often only in the latter, so it looks present to every # query and still fails to add: "group '$ACCESS_GROUP' does not # exist". Copying the line across is what makes usermod agree. - echo " 2. This is an rpm-ostree system (Bazzite/Silverblue), where" >&2 + echo " 3. This is an rpm-ostree system (Bazzite/Silverblue), where" >&2 echo " '$ACCESS_GROUP' may exist only in /usr/lib/group. usermod" >&2 echo " reads /etc/group alone and will refuse. Copy it over:" >&2 echo " grep -E '^$ACCESS_GROUP:' /usr/lib/group | sudo tee -a /etc/group" >&2 echo " then:" >&2 else - echo " 2. Add yourself to the group:" >&2 + echo " 3. Add yourself to the group:" >&2 fi echo " sudo usermod -aG $ACCESS_GROUP $USER_NAME" >&2 echo " and log out and back in for it to take effect." >&2 + elif [[ -z "$ACCESS_GROUP" ]]; then + echo >&2 + echo " 3. There is no group fallback on this system (no system-group" >&2 + echo " plugdev or dialout), which is normal on Fedora and its" >&2 + echo " derivatives. uaccess is the supported path there." >&2 fi fi echo -echo "Setup complete. Run the logger with:" +echo "Setup complete. Check the wiring end to end with:" +echo " ./.venv/bin/python diagnose-comms.py" +echo "then log:" echo " ./.venv/bin/python logger.py --duration 10" -echo "or, to check wiring only:" -echo " ./.venv/bin/python logger.py --scan-only" +echo "The test suite needs no hardware:" +echo " ./.venv/bin/python -m pytest" diff --git a/sweep.py b/sweep.py index f175374..ab84678 100644 --- a/sweep.py +++ b/sweep.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 """Sweep sample rates and report what each configuration actually delivers. - ./.venv/bin/python sweep.py # 1..256 Hz, powers of two + ./.venv/bin/python sweep.py # 2..512 Hz, powers of two ./.venv/bin/python sweep.py --from 8 --to 128 ./.venv/bin/python sweep.py --rates 10,25,50 @@ -21,6 +21,7 @@ import subprocess import sys from pathlib import Path +import ch347 import logger import rm3100 @@ -80,8 +81,8 @@ def main(): help="highest target rate in Hz (default: %(default)s)") ap.add_argument("--rates", help="explicit comma-separated rates, overriding " "the powers-of-two range") - ap.add_argument("--bus-speed", type=int, choices=[20, 100, 400, 750], - default=750) + ap.add_argument("--bus-speed", type=int, choices=sorted(ch347.SPEEDS), + default=ch347.DEFAULT_SPEED_KHZ) ap.add_argument("--tmrc", type=lambda s: int(s, 0), default=0x92, help="held fast so the cycle count governs (default: 0x92)") ap.add_argument("--duration", type=float, default=RECORD_SECONDS, diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..5824fba --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,205 @@ +"""Fakes standing in for the hardware, so the whole suite runs unplugged. + +Three levels are faked, matching the three seams in the code: + + FakeUsbDevice a CH347 at the libusb boundary -- ch347.py's own framing is + then the thing under test, byte for byte. + FakeBus an I2C master at the rm3100.py boundary, holding a register + map, so the driver's register sequences are under test. + FakeSensor an RM3100 at the logger.py boundary, replaying a scripted + timeline of DRDY transitions, so the miss-counting maths is + under test without waiting on real time. +""" + +import sys +from pathlib import Path + +import pytest + +# The modules under test sit at the repo root, next to this directory. +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) + + +class FakeUsbDevice: + """Records bulk writes and replays queued bulk reads.""" + + def __init__(self): + self.writes = [] # bytes sent to EP_OUT, in order + self.replies = [] # bytes to hand back, one per read() + self.short_write_by = 0 # subtract this from the reported write length + + def queue(self, *replies): + """Queue one reply per expected read, as bytes or an iterable of ints.""" + self.replies.extend(bytes(r) for r in replies) + return self + + def write(self, endpoint, data, timeout): + self.writes.append(bytes(data)) + return len(data) - self.short_write_by + + def read(self, endpoint, length, timeout): + if not self.replies: + raise AssertionError( + f"unexpected read of {length} bytes: nothing queued") + return bytearray(self.replies.pop(0)) + + @property + def last_write(self): + return self.writes[-1] + + +@pytest.fixture +def usb_device(monkeypatch): + """A CH347 that never touches USB. Yields the device; build the bus with it.""" + import usb.core + import usb.util + + device = FakeUsbDevice() + monkeypatch.setattr(usb.core, "find", lambda **kwargs: device) + monkeypatch.setattr(usb.util, "claim_interface", lambda *a: None) + monkeypatch.setattr(usb.util, "release_interface", lambda *a: None) + monkeypatch.setattr(usb.util, "dispose_resources", lambda *a: None) + return device + + +class FakeBus: + """An I2C bus holding one device's register map. + + Implements write/read/write_read with the same (data, mono, wall) contract + as ch347.CH347I2C, so rm3100.RM3100 cannot tell the difference. + """ + + def __init__(self, registers=None, address=0x23, combined=True): + self.address = address + self.registers = dict(registers or {}) + self.writes = [] # (addr, bytes) of every write + self.pointer = 0 + self.clock = 1000.0 + if not combined: + # Some buses cannot do a repeated START; the driver has a fallback + # path for them and it needs exercising too. read_reg() looks the + # method up with getattr(..., None), so hiding it this way is the + # same to the driver as never having had one. + self.write_read = None + + def _stamp(self): + self.clock += 0.001 + return self.clock, self.clock + 1_700_000_000.0 + + def write(self, addr, data): + data = bytes(data) + self.writes.append((addr, data)) + self.pointer = data[0] + for offset, value in enumerate(data[1:]): + self.registers[data[0] + offset] = value + + def read(self, addr, count): + mono, wall = self._stamp() + data = bytes(self.registers.get(self.pointer + i, 0) + for i in range(count)) + return data, mono, wall + + def write_read(self, addr, data, count): + self.write(addr, data) + return self.read(addr, count) + + +class FakeSensor: + """An RM3100 replaying a scripted DRDY timeline against a fake clock. + + `events` is a list of (monotonic_time, counts) pairs: the instants at which a + measurement completes and what it reads. poll_ready() reports DRDY set once + the clock has passed the next event and the previous result was consumed. + + Every poll and read advances the clock by a fixed cost, so a test states its + timeline in seconds and gets deterministic brackets out. + """ + + POLL_COST = 0.0003 + READ_COST = 0.0005 + + class Exhausted(Exception): + """Raised once max_polls is reached, to stop an unbounded run.""" + + def __init__(self, events, poll_cost=POLL_COST, read_cost=READ_COST, + stalls=None, max_polls=None): + self.events = list(events) + self.poll_cost = poll_cost + self.read_cost = read_cost + # [(at_time, extra_seconds), ...] -- the first poll starting at or after + # at_time pays extra_seconds before returning. A host stall is the only + # thing that actually loses a measurement on this rig, and stating it in + # seconds rather than in poll counts keeps a test independent of how + # many times the loop happens to poll. + self.stalls = sorted(stalls or []) + # A run with no --duration never returns on its own, so a test of that + # path needs the sensor to be the thing that stops it. + self.max_polls = max_polls + self.now = 0.0 + self.pending = None # counts of a completed, unread measurement + self.next_event = 0 + self.polls = 0 + self.reads = [] + + # -- clock plumbing ------------------------------------------------- + def monotonic(self): + return self.now + + def wall(self): + return self.now + 1_700_000_000.0 + + def _advance(self, dt): + self.now += dt + while (self.next_event < len(self.events) + and self.events[self.next_event][0] <= self.now): + # A completion overwrites any unread one: the chip does not queue, + # which is exactly why an unnoticed miss loses a measurement. + self.pending = self.events[self.next_event][1] + self.next_event += 1 + + # -- the RM3100 surface logger.py uses ------------------------------ + def poll_ready(self): + self.polls += 1 + if self.max_polls is not None and self.polls > self.max_polls: + raise self.Exhausted(f"stopped after {self.max_polls} polls") + stall = 0.0 + if self.stalls and self.now >= self.stalls[0][0]: + stall = self.stalls.pop(0)[1] + self._advance(self.poll_cost + stall) + return self.pending is not None, self.now, self.wall() + + def read_raw(self): + self._advance(self.read_cost) + counts = self.pending if self.pending is not None else (0, 0, 0) + self.pending = None + self.reads.append(counts) + return counts, self.now, self.wall() + + +@pytest.fixture +def fake_clock(monkeypatch): + """Point time.monotonic and time.time at a FakeSensor's clock. + + sample_loop() binds `monotonic = time.monotonic` when it runs, not at import, + so patching the module attribute is enough to control it. + """ + import time + + def install(sensor): + monkeypatch.setattr(time, "monotonic", sensor.monotonic) + monkeypatch.setattr(time, "time", sensor.wall) + return sensor + + return install + + +def grid_events(period, count, start=None, counts=(100, 200, 300), skip=()): + """Measurement completions on a uniform grid, optionally dropping some. + + A skipped index is a completion the chip still makes -- the grid never + pauses -- so it is present here; a test forces a *miss* by making the host + too slow to read it, not by removing it. + """ + start = period if start is None else start + return [(start + i * period, (counts[0] + i, counts[1], counts[2])) + for i in range(count) if i not in skip] diff --git a/tests/synthetic.py b/tests/synthetic.py new file mode 100644 index 0000000..02fdc26 --- /dev/null +++ b/tests/synthetic.py @@ -0,0 +1,86 @@ +"""Write synthetic captures, so tests never need hardware or a recorded file. + +Shared by test_capture.py, which uses it to break one guarantee at a time, and +test_compare.py, which uses it to plant a known answer -- a scale factor, a +rotation, a tone at a chosen fraction of the sample rate -- and check that the +analysis recovers it. + +Defaults produce a clean, loadable capture; every argument exists to change one +thing about it. +""" + +import numpy as np + +import capture +import rm3100 + +CYCLE_COUNT = 100 +DT = 1.0 / 250.0 # a deliberately non-nominal true period + + +def write_capture(path, rows=200, dt=DT, nominal_hz=282.0, header=None, + flags=None, drop_header_key=None, extra_lines=(), + index_from=0, index_step=1, start_time=1_700_000_000.0, + amplitude=1000.0, seed=0, cycle_count=CYCLE_COUNT, + counts=None): + """Write a synthetic capture and return its path. + + `counts` overrides the generated signal with an (rows, 3) array of raw + counts, which is how a test plants an exact answer. `cycle_count` moves the + header's gain and is what makes a decimation pair possible: two captures + whose cycle counts differ by an integer factor. + """ + meta = { + "rm3100_capture": 1, + "nominal_rate_hz": nominal_hz, + "tmrc_nominal_hz": 600.0, + "tmrc": "0x92", + "cycle_count": cycle_count, + "tesla_per_count": repr(rm3100.tesla_per_count(cycle_count)), + "i2c_address": "0x23", + "bus_speed_khz": 750, + "revid": "0x22", + "calibrated_period_s": repr(dt), + } + meta.update(header or {}) + if drop_header_key: + meta.pop(drop_header_key, None) + + rng = np.random.default_rng(seed) + flags = flags or {} + lines = [f"# {k}: {v}" for k, v in meta.items()] + lines += list(extra_lines) + lines.append("sample_index,system_time_unix,x_raw,y_raw,z_raw,warning") + for i in range(rows): + index = index_from + i * index_step + warning = flags.get(i, "") + if capture.WARN_MISSED in warning: + x = y = z = 0 + elif counts is not None: + x, y, z = (int(round(c)) for c in counts[i]) + else: + x = int(amplitude + rng.normal(0, 3)) + y = int(2 * amplitude + rng.normal(0, 3)) + z = int(-amplitude + rng.normal(0, 3)) + lines.append(f"{index},{start_time + index * dt:.6f},{x},{y},{z},{warning}") + path.write_text("\n".join(lines) + "\n") + return str(path) + + +def field_counts(rows, mean_nt, cycle_count, noise_nt=0.0, seed=0, tones=()): + """Raw counts for a field of a given mean, noise and planted tones. + + `mean_nt` is an (x, y, z) field in nanotesla; `tones` is a sequence of + (axis_index, cycles_per_sample, amplitude_nT) added on top. Quantisation to + integer counts is deliberate -- it is what the real file carries, and a test + that skipped it would not exercise the dither the analysis relies on. + """ + lsb = rm3100.tesla_per_count(cycle_count) * rm3100.NT_PER_TESLA + rng = np.random.default_rng(seed) + n = np.arange(rows) + nt = np.tile(np.asarray(mean_nt, dtype=float), (rows, 1)) + if noise_nt: + nt += rng.normal(0, noise_nt, size=(rows, 3)) + for axis, cycles_per_sample, amplitude in tones: + nt[:, axis] += amplitude * np.cos(2 * np.pi * cycles_per_sample * n) + return nt / lsb diff --git a/tests/test_analysis.py b/tests/test_analysis.py new file mode 100644 index 0000000..de389ac --- /dev/null +++ b/tests/test_analysis.py @@ -0,0 +1,221 @@ +"""Spectral and statistical helpers, checked against signals of known answer. + +Each test feeds in something whose spectrum or deviation is known analytically, +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 numpy as np +import pytest + +import characterize +import compare +import logger +import plot + + +# -------------------------------------------------------------------------- +# welch_asd +# -------------------------------------------------------------------------- + +def test_welch_asd_recovers_the_level_of_white_noise(): + """A white signal of sd s at rate fs sits at s/sqrt(fs/2) per root hertz.""" + fs, sd = 250.0, 20.0 + v = np.random.default_rng(1).normal(0, sd, 200_000) + freqs, asd = characterize.welch_asd(v, fs) + assert np.median(asd) == pytest.approx(sd / np.sqrt(fs / 2), rel=0.05) + + +def test_welch_asd_scales_with_amplitude_not_length(): + fs = 250.0 + rng = np.random.default_rng(2) + short = characterize.welch_asd(rng.normal(0, 10, 50_000), fs)[1] + long = characterize.welch_asd(rng.normal(0, 10, 200_000), fs)[1] + assert np.median(short) == pytest.approx(np.median(long), rel=0.1) + louder = characterize.welch_asd(rng.normal(0, 20, 50_000), fs)[1] + assert np.median(louder) == pytest.approx(2 * np.median(short), rel=0.1) + + +def test_welch_asd_puts_a_tone_in_the_right_bin(): + fs, tone = 250.0, 60.0 + t = np.arange(100_000) / fs + v = np.sin(2 * np.pi * tone * t) + freqs, asd = characterize.welch_asd(v, fs) + assert freqs[np.argmax(asd)] == pytest.approx(tone, abs=fs / 4096) + + +def test_welch_asd_drops_the_dc_bin(): + """A large DC offset must not appear as signal; detrending removes it.""" + fs = 250.0 + v = 50_000 + np.random.default_rng(3).normal(0, 1, 20_000) + freqs, asd = characterize.welch_asd(v, fs) + assert freqs[0] > 0 + assert asd.max() < 10 + + +def test_welch_asd_removes_a_linear_ramp(): + """A slow drift would otherwise smear energy across the low bins.""" + fs = 250.0 + n = 40_000 + rng = np.random.default_rng(4) + noise = rng.normal(0, 5, n) + ramped = noise + np.linspace(0, 5000, n) + flat_asd = characterize.welch_asd(noise, fs)[1] + ramp_asd = characterize.welch_asd(ramped, fs)[1] + assert np.median(ramp_asd) == pytest.approx(np.median(flat_asd), rel=0.05) + + +def test_welch_asd_frequencies_stop_at_nyquist(): + fs = 250.0 + freqs, _ = characterize.welch_asd( + np.random.default_rng(5).normal(0, 1, 20_000), fs) + assert freqs[-1] == pytest.approx(fs / 2) + assert len(freqs) == len(characterize.welch_asd( + np.random.default_rng(5).normal(0, 1, 20_000), fs)[1]) + + +def test_welch_asd_handles_a_capture_barely_long_enough(): + """capture.py's floor is 64 samples, so the spectrum code must survive it.""" + freqs, asd = characterize.welch_asd( + np.random.default_rng(6).normal(0, 1, 64), 250.0) + assert len(freqs) == len(asd) > 0 + assert np.all(np.isfinite(asd)) + + +# -------------------------------------------------------------------------- +# allan_deviation +# -------------------------------------------------------------------------- + +def test_allan_deviation_of_white_noise_falls_as_root_tau(): + """White noise gives slope -1/2 on a log-log ADEV plot.""" + fs, sd = 100.0, 10.0 + v = np.random.default_rng(7).normal(0, sd, 100_000) + taus, devs = characterize.allan_deviation(v, fs) + # Fit the log-log slope over the well-averaged decades. + keep = (taus > 10 / fs) & (taus < 1000 / fs) + slope = np.polyfit(np.log(taus[keep]), np.log(devs[keep]), 1)[0] + assert slope == pytest.approx(-0.5, abs=0.05) + + +def test_allan_deviation_starts_near_the_sample_sd(): + """At tau = one sample the deviation is the sd of the differences.""" + fs, sd = 100.0, 10.0 + v = np.random.default_rng(8).normal(0, sd, 50_000) + taus, devs = characterize.allan_deviation(v, fs) + assert taus[0] == pytest.approx(1 / fs) + assert devs[0] == pytest.approx(sd, rel=0.05) + + +def test_allan_deviation_turns_up_on_a_ramp(): + """Drift is what an upturn means; a pure ramp must produce one.""" + fs = 100.0 + n = 50_000 + v = np.random.default_rng(9).normal(0, 1, n) + np.linspace(0, 500, n) + taus, devs = characterize.allan_deviation(v, fs) + assert devs[-1] > devs[np.argmin(devs)] * 5 + + +def test_allan_deviation_of_a_constant_is_zero(): + taus, devs = characterize.allan_deviation(np.full(10_000, 42.0), 100.0) + assert np.allclose(devs, 0.0, atol=1e-9) + + +def test_allan_deviation_taus_increase_and_stay_in_range(): + taus, devs = characterize.allan_deviation( + np.random.default_rng(10).normal(0, 1, 10_000), 100.0) + assert np.all(np.diff(taus) > 0) + assert taus[-1] <= 10_000 / 4 / 100.0 + assert len(taus) == len(devs) + + +# -------------------------------------------------------------------------- +# rolling_mean +# -------------------------------------------------------------------------- + +def test_rolling_mean_matches_a_naive_partial_window_mean(): + """The point of the count normalisation is that the ends do not taper.""" + v = np.arange(50, dtype=float) + window = 7 + got = plot.rolling_mean(v, window) + for i in (0, 1, 25, 48, 49): + lo = max(0, i - window // 2) + hi = min(len(v), i + window // 2 + 1) + assert got[i] == pytest.approx(v[lo:hi].mean()) + + +def test_rolling_mean_preserves_a_constant_including_the_ends(): + v = np.full(100, 7.0) + assert np.allclose(plot.rolling_mean(v, 21), 7.0) + + +def test_rolling_mean_reduces_noise_by_root_window(): + v = np.random.default_rng(11).normal(0, 10, 100_000) + smoothed = plot.rolling_mean(v, 25) + assert smoothed.std() == pytest.approx(10 / np.sqrt(25), rel=0.1) + + +@pytest.mark.parametrize("window", [0, 1]) +def test_rolling_mean_is_a_no_op_below_two(window): + v = np.arange(10, dtype=float) + assert plot.rolling_mean(v, window) is v + + +# -------------------------------------------------------------------------- +# i2c_bus_time +# -------------------------------------------------------------------------- + +def test_i2c_bus_time_matches_the_hand_calculation(): + """A poll (1 byte) and a results read (9), each n+3 bytes of 9 bits + 3.""" + bits = (1 + 3) * 9 + 3 + (9 + 3) * 9 + 3 + assert bits == 150 + assert logger.i2c_bus_time(400) == pytest.approx(bits / 400_000.0) + + +def test_i2c_bus_time_is_inversely_proportional_to_speed(): + assert logger.i2c_bus_time(100) == pytest.approx( + 4 * logger.i2c_bus_time(400)) + assert logger.i2c_bus_time(750) == pytest.approx(0.200e-3, abs=5e-6) + + +def test_i2c_bus_time_at_the_default_is_a_small_share_of_the_period(): + """The documented 6% of a cc=100 period at 750 kHz.""" + import rm3100 + share = logger.i2c_bus_time(750) / rm3100.sample_period(100) + assert share == pytest.approx(0.06, abs=0.005) + + +# -------------------------------------------------------------------------- +# compare.axis_ratio_spread +# -------------------------------------------------------------------------- + +def _record(mean): + mean = np.array(mean, dtype=float) + return {"mean": mean, "field": float(np.linalg.norm(mean))} + + +def test_axis_ratio_spread_is_zero_for_a_pure_gain_change(): + a = _record([10_000, 20_000, -15_000]) + b = _record([10_600, 21_200, -15_900]) # every axis x1.06 + spread = compare.axis_ratio_spread(a, b, b["mean"] / a["mean"]) + assert spread == pytest.approx(0.0, abs=1e-9) + + +def test_axis_ratio_spread_detects_movement(): + a = _record([10_000, 20_000, -15_000]) + b = _record([11_400, 21_000, -13_000]) # each axis moved differently + spread = compare.axis_ratio_spread(a, b, b["mean"] / a["mean"]) + assert spread > compare.RATIO_SPREAD_OK + + +def test_axis_ratio_spread_ignores_an_axis_carrying_no_field(): + """A near-zero mean makes its ratio noise, which used to read as movement.""" + a = _record([10_000, 20_000, 5]) + b = _record([10_600, 21_200, -30]) # X and Y are a clean x1.06 + spread = compare.axis_ratio_spread(a, b, b["mean"] / a["mean"]) + assert spread == pytest.approx(0.0, abs=1e-9) + + +def test_axis_ratio_spread_abstains_with_too_few_usable_axes(): + a = _record([10_000, 3, 5]) + b = _record([10_600, -8, 2]) + assert compare.axis_ratio_spread(a, b, b["mean"] / a["mean"]) is None diff --git a/tests/test_capture.py b/tests/test_capture.py new file mode 100644 index 0000000..82d8371 --- /dev/null +++ b/tests/test_capture.py @@ -0,0 +1,344 @@ +"""Capture loading: the header contract, the flag columns, and hostile files. + +capture.py is the only reader of a capture, so every guarantee the format claims +-- contiguous index, parseable counts, recoverable time base -- has to be checked +here or nowhere. +""" + +import numpy as np +import pytest + +import capture +import rm3100 + +# The writer lives in synthetic.py because test_compare.py needs it too: +# one definition of what a valid capture looks like, not two that drift. +from synthetic import CYCLE_COUNT, DT, write_capture + + +# -------------------------------------------------------------------------- +# The happy path +# -------------------------------------------------------------------------- + +def test_load_round_trips_a_clean_capture(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv", rows=500)) + assert len(cap.sample_index) == 500 + assert cap.cycle_count == CYCLE_COUNT + assert not cap.missed.any() + assert not cap.ambiguous.any() + + +def test_counts_are_converted_with_the_header_constant(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv")) + assert cap.lsb_nt == pytest.approx( + rm3100.tesla_per_count(CYCLE_COUNT) * rm3100.NT_PER_TESLA) + # x_raw was written around 1000 counts. + assert cap.x.mean() == pytest.approx(1000 * cap.lsb_nt, rel=0.01) + + +def test_total_is_the_norm_of_the_three_axes(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv")) + assert cap.total == pytest.approx( + np.sqrt(cap.x**2 + cap.y**2 + cap.z**2)) + + +def test_true_period_is_recovered_from_the_host_clock(tmp_path): + """The whole point of the fit: the nominal rate is wrong and the file says so.""" + cap = capture.load(write_capture(tmp_path / "c.csv", dt=DT, nominal_hz=282.0)) + # Microsecond-resolution timestamps are the only error left in the slope. + assert cap.dt_true == pytest.approx(DT, rel=1e-7) + assert cap.true_rate_hz == pytest.approx(250.0, rel=1e-7) + + +def test_the_epoch_offset_does_not_degrade_the_fit(tmp_path): + """polyfit does not centre the ordinate, so the fit has to do it itself. + + Fitting raw epoch seconds against a millisecond slope costs six significant + digits to cancellation; a capture starting at t=0 and the same capture + starting in 2023 must recover the same period. + """ + at_zero = capture.load(write_capture(tmp_path / "zero.csv", rows=300, + start_time=0.0)) + at_epoch = capture.load(write_capture(tmp_path / "epoch.csv", rows=300, + start_time=1_700_000_000.0)) + assert at_epoch.dt_true == pytest.approx(at_zero.dt_true, rel=1e-7) + + +def test_rate_error_is_signed_like_the_rate_not_the_period(tmp_path): + """Chip slower than nominal must read negative, matching the reported Hz.""" + slow = capture.load(write_capture(tmp_path / "slow.csv", dt=1 / 250.0, + nominal_hz=282.0)) + fast = capture.load(write_capture(tmp_path / "fast.csv", dt=1 / 300.0, + nominal_hz=282.0)) + assert slow.rate_error < 0 + assert fast.rate_error > 0 + 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): + 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) + + +def test_duration_spans_the_capture(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv", rows=500, dt=DT)) + assert cap.duration == pytest.approx(499 * DT) + + +def test_axes_lists_three_axes_plus_the_derived_total(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv")) + assert [k for k, _ in cap.axes()] == ["x", "y", "z", "total"] + + +def test_note_reaches_the_summary(tmp_path): + path = write_capture(tmp_path / "c.csv", header={"note": "3V0 LDO"}) + assert "3V0 LDO" in capture.load(path).summary() + + +# -------------------------------------------------------------------------- +# MISSED and AMBIGUOUS +# -------------------------------------------------------------------------- + +def test_missed_rows_are_flagged_and_interpolated(tmp_path): + path = write_capture(tmp_path / "c.csv", rows=200, + flags={50: capture.WARN_MISSED}) + cap = capture.load(path) + assert cap.missed.sum() == 1 + assert cap.missed[50] + # The placeholder's zeros must not survive into the data. + assert cap.x[50] == pytest.approx((cap.x[49] + cap.x[51]) / 2, rel=0.01) + + +def test_consecutive_missed_rows_interpolate_across_the_whole_gap(tmp_path): + path = write_capture(tmp_path / "c.csv", rows=200, + flags={i: capture.WARN_MISSED for i in (60, 61, 62)}) + cap = capture.load(path) + assert cap.missed.sum() == 3 + assert cap.x[59] < cap.x[62] or cap.x[59] > cap.x[62] # finite, not zero + assert np.all(np.abs(cap.x[60:63]) > 0) + + +def test_ambiguous_is_independent_of_missed(tmp_path): + """A gap can round to zero losses and still be a coin toss.""" + path = write_capture(tmp_path / "c.csv", rows=200, flags={ + 70: capture.WARN_AMBIGUOUS, + 80: f"{capture.WARN_MISSED} {capture.WARN_AMBIGUOUS}", + }) + cap = capture.load(path) + assert cap.ambiguous.sum() == 2 + assert cap.missed.sum() == 1 + assert cap.ambiguous[70] and not cap.missed[70] # real row, flagged + assert cap.ambiguous[80] and cap.missed[80] # placeholder, both + + +def test_summary_reports_both_flag_kinds(tmp_path): + path = write_capture(tmp_path / "c.csv", rows=200, flags={ + 70: capture.WARN_AMBIGUOUS, 80: capture.WARN_MISSED}) + text = capture.load(path).summary() + assert "lost measurement" in text + assert "AMBIGUOUS" in text + + +def test_flag_order_within_the_column_does_not_matter(tmp_path): + path = write_capture(tmp_path / "c.csv", rows=200, flags={ + 90: f"{capture.WARN_AMBIGUOUS} {capture.WARN_MISSED}"}) + cap = capture.load(path) + assert cap.missed[90] and cap.ambiguous[90] + + +# -------------------------------------------------------------------------- +# Drift detection +# -------------------------------------------------------------------------- + +def test_a_uniform_capture_is_not_drift_limited(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv", rows=500)) + assert not cap.drift_limited + assert cap.residual_sd < cap.dt_true + + +def test_a_rate_that_changes_mid_run_is_drift_limited(tmp_path): + """One slope cannot describe two rates; that has to be visible.""" + path = tmp_path / "drift.csv" + write_capture(path, rows=400) + lines = path.read_text().splitlines() + head = [l for l in lines if l.startswith("#")] + body = [l for l in lines if not l.startswith("#")] + out = head + [body[0]] + t = 1_700_000_000.0 + for i, row in enumerate(body[1:]): + fields = row.split(",") + # Second half runs 30% slower -- far more than any real oscillator, but + # this is testing the detector, not the sensor. + t += DT if i < 200 else DT * 1.3 + out.append(",".join([fields[0], f"{t:.6f}"] + fields[2:])) + path.write_text("\n".join(out) + "\n") + + cap = capture.load(str(path)) + assert cap.drift_limited + assert "single-rate model does not fit" in cap.summary() + + +def test_read_jitter_is_reported_apart_from_the_fit_residual(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv", rows=500)) + assert cap.read_jitter >= 0.0 + assert cap.residual_sd >= 0.0 + + +# -------------------------------------------------------------------------- +# restrict() +# -------------------------------------------------------------------------- + +def test_restrict_refits_on_the_window(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv", rows=1000)) + window = cap.restrict(0.5, 1.5) + assert len(window.sample_index) < len(cap.sample_index) + assert window.elapsed[0] >= 0.0 + assert window.dt_true == pytest.approx(cap.dt_true, rel=1e-6) + + +def test_restrict_accepts_an_open_ended_window(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv", rows=1000)) + assert len(cap.restrict(start=0.5).sample_index) < 1000 + assert len(cap.restrict(end=1.0).sample_index) < 1000 + + +def test_restrict_refuses_a_window_too_small_to_analyse(tmp_path): + cap = capture.load(write_capture(tmp_path / "c.csv", rows=500)) + with pytest.raises(capture.CaptureError, match="too few samples|too few"): + cap.restrict(0.0, 0.01) + + +def test_restrict_carries_the_flags_through(tmp_path): + path = write_capture(tmp_path / "c.csv", rows=1000, + flags={500: capture.WARN_MISSED}) + cap = capture.load(path) + window = cap.restrict(400 * DT, 600 * DT) + assert window.missed.sum() == 1 + + +# -------------------------------------------------------------------------- +# Rejected files -- each guarantee the format makes, tested by breaking it +# -------------------------------------------------------------------------- + +def test_a_file_without_the_capture_header_is_rejected(tmp_path): + path = tmp_path / "old.csv" + path.write_text("sample_index,system_time_unix,x_raw,y_raw,z_raw,warning\n" + + "".join(f"{i},{i * 0.004},1,2,3,\n" for i in range(100))) + with pytest.raises(capture.CaptureError, match="no capture header"): + capture.load(str(path)) + + +@pytest.mark.parametrize("key", ["cycle_count", "tesla_per_count", + "nominal_rate_hz"]) +def test_a_header_missing_a_required_key_is_rejected(tmp_path, key): + path = write_capture(tmp_path / "c.csv", drop_header_key=key) + with pytest.raises(capture.CaptureError, match=f"header missing.*{key}"): + capture.load(path) + + +def test_too_few_samples_to_analyse_is_rejected(tmp_path): + path = write_capture(tmp_path / "short.csv", rows=10) + with pytest.raises(capture.CaptureError, match="too few"): + capture.load(path) + + +def test_a_capture_without_the_warning_column_is_rejected(tmp_path): + path = tmp_path / "nowarn.csv" + body = "".join(f"{i},{1_700_000_000 + i * 0.004:.6f},1,2,3\n" + for i in range(100)) + path.write_text(f"# rm3100_capture: 1\n# nominal_rate_hz: 282\n" + f"# cycle_count: 100\n" + f"# tesla_per_count: {rm3100.tesla_per_count(100)!r}\n" + "sample_index,system_time_unix,x_raw,y_raw,z_raw\n" + body) + with pytest.raises(capture.CaptureError, match="no 'warning' column"): + capture.load(str(path)) + + +def test_a_non_contiguous_index_is_rejected(tmp_path): + """A hole means the grid is broken, so the whole time base is unusable.""" + path = tmp_path / "holed.csv" + write_capture(path, rows=200) + lines = path.read_text().splitlines() + del lines[100] # drop one data row + path.write_text("\n".join(lines) + "\n") + with pytest.raises(capture.CaptureError, match="not contiguous"): + capture.load(str(path)) + + +def test_a_decreasing_index_is_rejected(tmp_path): + path = write_capture(tmp_path / "back.csv", rows=100, index_step=-1, + index_from=500) + with pytest.raises(capture.CaptureError, match="not contiguous"): + capture.load(path) + + +def test_an_all_missed_capture_is_rejected(tmp_path): + path = write_capture(tmp_path / "empty.csv", rows=100, + flags={i: capture.WARN_MISSED for i in range(100)}) + with pytest.raises(capture.CaptureError, match="every row is a lost"): + capture.load(path) + + +def test_a_missing_file_raises_oserror_not_captureerror(tmp_path): + with pytest.raises(OSError): + capture.load(str(tmp_path / "nope.csv")) + + +# -------------------------------------------------------------------------- +# Malformed and adversarial input +# -------------------------------------------------------------------------- + +def test_a_header_line_without_a_colon_is_ignored(tmp_path): + path = write_capture(tmp_path / "c.csv", + extra_lines=["# just a comment with no key"]) + cap = capture.load(path) + assert "just a comment with no key" not in str(cap.meta) + + +def test_an_index_starting_away_from_zero_still_loads(tmp_path): + """A capture cut from the middle of a run is still a valid grid.""" + cap = capture.load(write_capture(tmp_path / "c.csv", rows=200, + index_from=1_000_000)) + assert cap.sample_index[0] == 1_000_000 + assert cap.elapsed[0] == pytest.approx(1_000_000 * cap.dt_true) + + +def test_a_huge_index_does_not_lose_precision(tmp_path): + """int64 indices must not silently become floats in the fit.""" + cap = capture.load(write_capture(tmp_path / "c.csv", rows=200, + index_from=2**40)) + assert cap.sample_index.dtype == np.int64 + assert int(cap.sample_index[-1] - cap.sample_index[0]) == 199 + + +def test_a_non_integer_raw_count_is_rejected(tmp_path): + path = tmp_path / "bad.csv" + write_capture(path, rows=100) + text = path.read_text().replace(",1000,", ",not-a-number,", 1) + path.write_text(text) + with pytest.raises(ValueError): + capture.load(str(path)) + + +def test_a_note_containing_delimiters_survives_the_header(tmp_path): + """The header is 'key: value' split on the first colon only.""" + path = write_capture(tmp_path / "c.csv", + header={"note": "supply: 3.0 V, run #2"}) + assert capture.load(path).meta["note"] == "supply: 3.0 V, run #2" + + +def test_an_unknown_warning_flag_is_ignored(tmp_path): + """The column generalises to future flags, so today's reader must not choke.""" + path = write_capture(tmp_path / "c.csv", rows=200, + flags={40: "SOMETHING_NEW"}) + cap = capture.load(path) + assert not cap.missed.any() + assert not cap.ambiguous.any() + + +def test_a_flag_substring_does_not_count_as_the_flag(tmp_path): + """Flags are whitespace-separated tokens, not substrings.""" + path = write_capture(tmp_path / "c.csv", rows=200, flags={40: "NOTMISSED"}) + assert not capture.load(path).missed.any() diff --git a/tests/test_ch347.py b/tests/test_ch347.py new file mode 100644 index 0000000..b73e929 --- /dev/null +++ b/tests/test_ch347.py @@ -0,0 +1,262 @@ +"""CH347 stream framing, byte for byte. + +The wire format is the one thing here with no second source: it was reverse +engineered from the aystarik Linux driver, and a silent change to it would show +up as unexplained NACKs on real hardware rather than as an exception. So these +tests assert the exact bytes rather than "a write happened". +""" + +import pytest + +import ch347 + +ADDR = 0x23 +ACK = ch347.ACK + + +@pytest.fixture +def bus(usb_device): + """A CH347I2C whose set_speed handshake has already been accounted for.""" + b = ch347.CH347I2C() + usb_device.writes.clear() + return b + + +# -------------------------------------------------------------------------- +# Construction +# -------------------------------------------------------------------------- + +def test_open_sets_the_requested_speed(usb_device): + ch347.CH347I2C(ch347.SPEED_400KHZ) + assert usb_device.writes == [ + bytes([ch347.CMD_STREAM, ch347.CMD_SET | ch347.SPEED_400KHZ, + ch347.CMD_END])] + + +def test_open_defaults_to_the_documented_speed(usb_device): + ch347.CH347I2C() + expected = ch347.SPEEDS[ch347.DEFAULT_SPEED_KHZ] + assert usb_device.last_write[1] == ch347.CMD_SET | expected + + +def test_speeds_table_covers_every_speed_constant(): + assert sorted(ch347.SPEEDS.values()) == [ + ch347.SPEED_20KHZ, ch347.SPEED_100KHZ, + ch347.SPEED_400KHZ, ch347.SPEED_750KHZ] + assert ch347.DEFAULT_SPEED_KHZ in ch347.SPEEDS + + +def test_missing_adapter_names_the_device(monkeypatch): + import usb.core + monkeypatch.setattr(usb.core, "find", lambda **kwargs: None) + with pytest.raises(IOError, match="1a86:55db"): + ch347.CH347I2C() + + +def test_claim_failure_points_at_setup(monkeypatch, usb_device): + import usb.core + import usb.util + + def refuse(*args): + raise usb.core.USBError("Access denied") + + monkeypatch.setattr(usb.util, "claim_interface", refuse) + with pytest.raises(IOError, match="setup.sh"): + ch347.CH347I2C() + + +def test_set_speed_rejects_an_unknown_rate(bus): + with pytest.raises(ValueError, match="expected 0-3"): + bus.set_speed(9) + + +# -------------------------------------------------------------------------- +# write() +# -------------------------------------------------------------------------- + +def test_write_frames_a_register_write(bus, usb_device): + usb_device.queue([ACK, ACK, ACK]) + bus.write(ADDR, [0x04, 0x64]) + assert usb_device.last_write == bytes([ + ch347.CMD_STREAM, + ch347.CMD_STA, + ch347.CMD_OUT | 3, # address plus two payload bytes + ADDR << 1, # write: LSB clear + 0x04, 0x64, + ch347.CMD_STO, + ch347.CMD_END, + ]) + + +def test_write_raises_on_an_address_nack(bus, usb_device): + usb_device.queue([0, ACK]) + with pytest.raises(IOError, match="No ACK from I2C address 0x23 on write"): + bus.write(ADDR, [0x04]) + + +def test_write_raises_on_a_data_nack(bus, usb_device): + usb_device.queue([ACK, 0]) + with pytest.raises(IOError, match="NACKed a data byte"): + bus.write(ADDR, [0x04]) + + +def test_write_rejects_an_oversized_payload(bus): + with pytest.raises(ValueError, match="exceeds the CH347 limit"): + bus.write(ADDR, bytes(ch347.MAX_XFER)) + + +def test_write_accepts_the_largest_legal_payload(bus, usb_device): + payload = bytes(ch347.MAX_XFER - 1) + usb_device.queue([ACK] * (len(payload) + 1)) + bus.write(ADDR, payload) + # The count field is 6 bits and must not have overflowed into the opcode. + assert usb_device.last_write[2] == ch347.CMD_OUT | ch347.MAX_XFER + + +# -------------------------------------------------------------------------- +# read() +# -------------------------------------------------------------------------- + +def test_read_acks_all_but_the_last_byte(bus, usb_device): + usb_device.queue([ACK, 0x11, 0x22, 0x33]) + data, mono, wall = bus.read(ADDR, 3) + assert data == bytes([0x11, 0x22, 0x33]) + assert usb_device.last_write == bytes([ + ch347.CMD_STREAM, + ch347.CMD_STA, + ch347.CMD_OUT | 1, + (ADDR << 1) | 1, # read: LSB set + ch347.CMD_IN | 2, # first n-1 bytes, ACKed + ch347.CMD_IN, # final byte, NACKed to end the read + ch347.CMD_STO, + ch347.CMD_END, + ]) + assert wall > mono # epoch seconds against a monotonic count + + +def test_single_byte_read_omits_the_acked_run(bus, usb_device): + usb_device.queue([ACK, 0x22]) + assert bus.read(ADDR, 1)[0] == b"\x22" + assert usb_device.last_write == bytes([ + ch347.CMD_STREAM, ch347.CMD_STA, ch347.CMD_OUT | 1, (ADDR << 1) | 1, + ch347.CMD_IN, ch347.CMD_STO, ch347.CMD_END]) + + +def test_read_raises_on_an_address_nack(bus, usb_device): + usb_device.queue([0, 0x00]) + with pytest.raises(IOError, match="No ACK from I2C address 0x23 on read"): + bus.read(ADDR, 1) + + +@pytest.mark.parametrize("count", [0, -1, ch347.MAX_XFER + 1]) +def test_read_rejects_an_illegal_length(bus, count): + with pytest.raises(ValueError, match="outside 1.."): + bus.read(ADDR, count) + + +# -------------------------------------------------------------------------- +# write_read() -- the repeated START that halves the round trips +# -------------------------------------------------------------------------- + +def test_write_read_uses_a_repeated_start(bus, usb_device): + usb_device.queue([ACK, ACK, ACK] + [0x00] * 9) + data, _, _ = bus.write_read(ADDR, [0x24], 9) + assert len(data) == 9 + assert usb_device.last_write == bytes([ + ch347.CMD_STREAM, + ch347.CMD_STA, + ch347.CMD_OUT | 2, # write address plus the register pointer + ADDR << 1, + 0x24, + ch347.CMD_STA, # repeated START -- no STOP in between + ch347.CMD_OUT | 1, + (ADDR << 1) | 1, + ch347.CMD_IN | 8, + ch347.CMD_IN, + ch347.CMD_STO, + ch347.CMD_END, + ]) + + +def test_write_read_expects_one_ack_per_clocked_out_byte(bus, usb_device): + """Two addresses plus the payload -- getting this wrong misaligns the data.""" + usb_device.queue([ACK, ACK, ACK, ACK, 0xAA]) + data, _, _ = bus.write_read(ADDR, [0x01, 0x02], 1) + assert data == b"\xaa" + + +def test_write_read_raises_when_any_phase_nacks(bus, usb_device): + usb_device.queue([ACK, ACK, 0, 0x00]) + with pytest.raises(IOError, match="NACKed during combined transfer"): + bus.write_read(ADDR, [0x24], 1) + + +def test_write_read_requires_something_to_write(bus): + with pytest.raises(ValueError, match="at least one byte"): + bus.write_read(ADDR, b"", 1) + + +def test_write_read_validates_both_lengths(bus): + with pytest.raises(ValueError, match="exceeds the CH347 limit"): + bus.write_read(ADDR, bytes(ch347.MAX_XFER), 1) + with pytest.raises(ValueError, match="outside 1.."): + bus.write_read(ADDR, [0x24], ch347.MAX_XFER + 1) + + +# -------------------------------------------------------------------------- +# probe() / scan() +# -------------------------------------------------------------------------- + +def test_probe_reports_an_ack(bus, usb_device): + usb_device.queue([ACK]) + assert bus.probe(ADDR) is True + assert usb_device.last_write == bytes([ + ch347.CMD_STREAM, ch347.CMD_STA, ch347.CMD_OUT | 1, ADDR << 1, + ch347.CMD_STO, ch347.CMD_END]) + + +def test_probe_treats_a_nack_as_absence_not_an_error(bus, usb_device): + usb_device.queue([0]) + assert bus.probe(ADDR) is False + + +def test_probe_swallows_a_usb_error(bus, monkeypatch): + def boom(*args, **kwargs): + raise IOError("pipe error") + + monkeypatch.setattr(bus, "_xfer", boom) + assert bus.probe(ADDR) is False + + +def test_scan_returns_only_responders(bus, usb_device): + # 0x08..0x0a: only the middle one ACKs. + usb_device.queue([0], [ACK], [0]) + assert bus.scan(first=0x08, last=0x0a) == [0x09] + + +# -------------------------------------------------------------------------- +# Transport errors +# -------------------------------------------------------------------------- + +def test_short_usb_write_is_an_error(bus, usb_device): + usb_device.short_write_by = 1 + usb_device.queue([ACK, ACK]) + with pytest.raises(IOError, match="Short USB write"): + bus.write(ADDR, [0x04]) + + +def test_short_usb_read_is_an_error(bus, usb_device): + usb_device.queue([ACK]) # one byte where two were expected + with pytest.raises(IOError, match="Short USB read"): + bus.write(ADDR, [0x04]) + + +def test_close_is_idempotent(bus): + bus.close() + bus.close() + + +def test_context_manager_releases_the_interface(usb_device): + with ch347.CH347I2C() as b: + assert b._claimed + assert not b._claimed diff --git a/tests/test_compare.py b/tests/test_compare.py new file mode 100644 index 0000000..2157372 --- /dev/null +++ b/tests/test_compare.py @@ -0,0 +1,650 @@ +"""Cross-capture analysis: the maths that turns several captures into a claim. + +Everything here is checked against a planted answer rather than against a +recorded file, because the point of each function is that it recovers something +specific -- a scale factor, an oscillator frequency, a tone at a chosen fraction +of the sample rate -- and only a synthetic capture knows what that was. + +The argument parsing is tested just as hard. `LABEL=path` has to survive a path +containing '=', and --band and --supply have to reject nonsense rather than +quietly produce a figure that means nothing. +""" + +import argparse +import math + +import numpy as np +import pytest + +import capture +import characterize as ch +import compare +import rm3100 +from synthetic import field_counts, write_capture + +# A field well off the axes, so no component sits near zero except where a test +# puts one there deliberately. +FIELD = (-40000.0, 16000.0, 12000.0) + + +def make(path, rows=4096, cycle_count=100, dt=1 / 300.0, mean=FIELD, + noise_nt=20.0, seed=0, tones=(), **kwargs): + counts = field_counts(rows, mean, cycle_count, noise_nt, seed, tones) + return capture.load(write_capture(path, rows=rows, dt=dt, + cycle_count=cycle_count, counts=counts, + **kwargs)) + + +# -------------------------------------------------------------------------- +# white_sd and band_stats +# -------------------------------------------------------------------------- + +def test_white_sd_recovers_the_sd_of_white_noise(): + v = np.random.default_rng(1).normal(0, 7.0, 200_000) + assert ch.white_sd(v) == pytest.approx(7.0, rel=0.02) + + +def test_white_sd_ignores_a_ramp_that_dominates_the_plain_sd(): + """The whole reason it exists: sd measures drift, this does not.""" + rng = np.random.default_rng(2) + noise = rng.normal(0, 5.0, 50_000) + drifting = noise + np.linspace(0, 500, 50_000) + assert drifting.std() > 100 # sd is all ramp + assert ch.white_sd(drifting) == pytest.approx(5.0, rel=0.03) + + +def test_white_sd_of_a_constant_is_zero_and_of_one_sample_is_defined(): + assert ch.white_sd(np.ones(500)) == 0.0 + assert ch.white_sd(np.array([1.0])) == 0.0 + + +def test_band_stats_median_matches_the_white_noise_level(): + fs = 300.0 + v = np.random.default_rng(3).normal(0, 10.0, 100_000) + # A flat spectrum of sd s over a one-sided band of fs/2 sits at this ASD. + expected = 10.0 / np.sqrt(fs / 2) + median, _, _, _ = ch.band_stats(v, fs, (3.0, 100.0)) + assert median == pytest.approx(expected, rel=0.05) + + +def test_band_stats_finds_a_planted_line_and_reports_its_frequency(): + fs, n = 300.0, 100_000 + t = np.arange(n) / fs + v = np.random.default_rng(4).normal(0, 1.0, n) + 50 * np.cos(2 * np.pi * 40 * t) + _, _, peak, peak_hz = ch.band_stats(v, fs, (3.0, 100.0)) + assert peak_hz == pytest.approx(40.0, abs=0.5) + assert peak > 20 + + +def test_band_stats_returns_nan_when_the_band_holds_no_bins(): + """A band above Nyquist has nothing in it, and must say so, not guess.""" + v = np.random.default_rng(5).normal(0, 1.0, 10_000) + median, rms, peak, peak_hz = ch.band_stats(v, 300.0, (200.0, 250.0)) + assert all(math.isnan(x) for x in (median, rms, peak, peak_hz)) + + +def test_band_for_scales_with_the_rate(): + lo, hi = ch.band_for(300.0) + assert lo == ch.BAND_LO_HZ + assert hi == pytest.approx(ch.BAND_NYQUIST_FRACTION * 150.0) + + +def test_the_band_stays_inside_the_decimation_filter_passband(): + """The constants have to hold this relationship or comparisons are biased. + + A band reaching past the anti-alias corner scores the decimated path partway + down a rolloff, which reads as a quieter sensor rather than a narrower + filter. + """ + assert ch.BAND_NYQUIST_FRACTION / 2 < ch.DECIMATE_CUTOFF_FRACTION < 0.5 + + +# -------------------------------------------------------------------------- +# fir_lowpass +# -------------------------------------------------------------------------- + +def response(h, cycles_per_sample): + n = np.arange(len(h)) - (len(h) - 1) / 2 + return abs(np.sum(h * np.exp(-2j * np.pi * cycles_per_sample * n))) + + +def test_fir_lowpass_has_unit_gain_at_dc(): + """Decimating must not rescale the field.""" + h = ch.fir_lowpass(0.1, 257) + assert h.sum() == pytest.approx(1.0) + assert response(h, 0.0) == pytest.approx(1.0) + + +def test_fir_lowpass_passes_below_the_corner_and_stops_above_it(): + h = ch.fir_lowpass(0.1, 513) + assert response(h, 0.05) == pytest.approx(1.0, abs=0.01) + assert response(h, 0.09) == pytest.approx(1.0, abs=0.02) + # Blackman buys a deep stopband; anything near -74 dB or better will do. + assert 20 * np.log10(response(h, 0.15)) < -60 + assert 20 * np.log10(response(h, 0.30)) < -60 + + +def test_fir_lowpass_is_linear_phase(): + h = ch.fir_lowpass(0.1, 129) + assert h == pytest.approx(h[::-1]) + + +def test_fir_lowpass_forces_an_odd_length_so_the_delay_is_a_whole_sample(): + assert len(ch.fir_lowpass(0.1, 128)) == 129 + + +@pytest.mark.parametrize("cutoff", [0.0, -0.1, 0.5, 0.6]) +def test_fir_lowpass_rejects_a_cutoff_outside_the_open_unit_band(cutoff): + with pytest.raises(ValueError, match="cutoff"): + ch.fir_lowpass(cutoff, 129) + + +def test_fir_lowpass_rejects_a_length_too_short_to_filter(): + with pytest.raises(ValueError, match="too short"): + ch.fir_lowpass(0.1, 1) + + +# -------------------------------------------------------------------------- +# decimate +# -------------------------------------------------------------------------- + +def test_boxcar_decimation_divides_white_noise_sd_by_root_k(): + v = np.random.default_rng(6).normal(0, 12.0, 400_000) + out = ch.decimate(v, 4, "boxcar") + assert len(out) == 100_000 + assert out.std() == pytest.approx(12.0 / 2, rel=0.02) + + +def test_decimation_preserves_the_spectral_density_of_white_noise(): + """The claim the whole recommendation rests on: same floor, fewer samples.""" + fs = 300.0 + v = np.random.default_rng(7).normal(0, 12.0, 200_000) + band = (3.0, 20.0) + before, _, _, _ = ch.band_stats(v, fs, band) + for method in ("boxcar", "fir"): + after, _, _, _ = ch.band_stats(ch.decimate(v, 4, method), fs / 4, band) + assert after == pytest.approx(before, rel=0.05), method + + +def test_decimate_by_one_is_a_no_op(): + v = np.random.default_rng(8).normal(0, 1.0, 100) + assert ch.decimate(v, 1) == pytest.approx(v) + + +def test_the_fir_removes_an_out_of_band_tone_that_the_boxcar_folds_in(): + """The difference that matters: a boxcar has a poor stopband, so it aliases. + + A tone at 0.3 cycles/sample is above the decimated Nyquist of 0.125 and + folds to 1/5 of the new rate. The chip's own integration is a boxcar, which + is exactly why sampling slowly cannot reject what sampling fast and + filtering can. + """ + n = 40_000 + tone = 300.0 * np.cos(2 * np.pi * 0.3 * np.arange(n)) + v = np.random.default_rng(9).normal(0, 20.0, n) + tone + + folded = ch.sample_locked_lines(ch.decimate(v, 4, "boxcar")) + assert any(l.numerator == 1 and l.period == 5 for l in folded) + + filtered = ch.sample_locked_lines(ch.decimate(v, 4, "fir")) + assert not any(l.numerator == 1 and l.period == 5 for l in filtered) + + +@pytest.mark.parametrize("k", [0, -2]) +def test_decimate_rejects_a_factor_below_one(k): + with pytest.raises(ValueError, match="at least 1"): + ch.decimate(np.zeros(100), k) + + +def test_decimate_rejects_an_unknown_method(): + with pytest.raises(ValueError, match="unknown decimation method"): + ch.decimate(np.zeros(100), 2, "bilinear") + + +def test_decimate_refuses_rather_than_returning_filter_transient(): + """Too few samples for the filter is a failure, not a short answer.""" + with pytest.raises(ValueError, match="too few"): + ch.decimate(np.zeros(200), 4, "fir") + + +def test_boxcar_refuses_when_there_is_not_even_one_full_group(): + with pytest.raises(ValueError, match="cannot be decimated"): + ch.decimate(np.zeros(3), 4, "boxcar") + + +# -------------------------------------------------------------------------- +# sample_locked_lines +# -------------------------------------------------------------------------- + +def test_sample_locked_lines_recovers_a_planted_tone_and_its_amplitude(): + n = 60_000 + v = (np.random.default_rng(10).normal(0, 20.0, n) + + 3.0 * np.cos(2 * np.pi * 0.25 * np.arange(n))) + found = ch.sample_locked_lines(v) + quarter = [l for l in found if (l.numerator, l.period) == (1, 4)] + assert quarter, "a 3 nT tone under 20 nT of noise should still be found" + assert quarter[0].amplitude == pytest.approx(3.0, rel=0.15) + assert quarter[0].sigma > 10 + + +def test_sample_locked_lines_recovers_a_tone_at_nyquist(): + n = 60_000 + v = (np.random.default_rng(11).normal(0, 20.0, n) + + 2.0 * (-1.0) ** np.arange(n)) + half = [l for l in ch.sample_locked_lines(v) + if (l.numerator, l.period) == (1, 2)] + assert half and half[0].amplitude == pytest.approx(2.0, rel=0.15) + + +@pytest.mark.parametrize("seed", range(6)) +def test_sample_locked_lines_stays_silent_on_white_noise(seed): + v = np.random.default_rng(100 + seed).normal(0, 20.0, 60_000) + assert ch.sample_locked_lines(v) == [] + + +def test_sample_locked_lines_is_not_fooled_by_drift(): + n = 60_000 + t = np.linspace(0, 1, n) + v = np.random.default_rng(12).normal(0, 5.0, n) + 4000 * t ** 3 - 900 * t + assert ch.sample_locked_lines(v) == [] + + +def test_sample_locked_lines_reports_each_frequency_once(): + """A period-4 tone is also period-8 and period-12; those add nothing.""" + n = 60_000 + v = (np.random.default_rng(13).normal(0, 10.0, n) + + 5.0 * np.cos(2 * np.pi * 0.25 * np.arange(n))) + fractions = [(l.numerator, l.period) for l in ch.sample_locked_lines(v)] + assert (2, 8) not in fractions and (3, 12) not in fractions + assert len(fractions) == len(set(fractions)) + # Every reported fraction is in lowest terms. + assert all(math.gcd(j, p) == 1 for j, p in fractions) + + +def test_sample_locked_lines_returns_strongest_first(): + n = 60_000 + index = np.arange(n) + v = (np.random.default_rng(14).normal(0, 10.0, n) + + 6.0 * np.cos(2 * np.pi * 0.25 * index) + + 2.0 * np.cos(2 * np.pi * index / 3)) + found = ch.sample_locked_lines(v) + assert [l.sigma for l in found] == sorted((l.sigma for l in found), + reverse=True) + + +def test_sample_locked_lines_abstains_on_a_series_too_short_to_fold(): + assert ch.sample_locked_lines(np.arange(10.0)) == [] + + +def test_sample_locked_lines_abstains_on_a_constant(): + assert ch.sample_locked_lines(np.ones(1000)) == [] + + +# -------------------------------------------------------------------------- +# trimmed +# -------------------------------------------------------------------------- + +def test_trimmed_drops_the_requested_seconds_from_both_ends(tmp_path): + cap = make(tmp_path / "c.csv", rows=3000, dt=1 / 100.0) # 30 s + kept, note = ch.trimmed(cap, 5.0) + assert kept.duration == pytest.approx(cap.duration - 10.0, abs=0.05) + assert "trimmed 5 s" in note + + +def test_trimmed_is_a_no_op_when_no_seconds_are_asked_for(tmp_path): + cap = make(tmp_path / "c.csv", rows=3000, dt=1 / 100.0) + kept, note = ch.trimmed(cap, 0.0) + assert kept is cap and note == "" + + +def test_trimmed_refuses_rather_than_gutting_a_short_capture(tmp_path): + """A deliberately short capture is legitimate; silently emptying it is not.""" + cap = make(tmp_path / "c.csv", rows=1000, dt=1 / 100.0) # 10 s + kept, note = ch.trimmed(cap, 30.0) + assert kept is cap + assert "not trimming" in note + + +def test_trimmed_rejects_a_negative_window(tmp_path): + cap = make(tmp_path / "c.csv", rows=3000, dt=1 / 100.0) + with pytest.raises(capture.CaptureError, match="negative"): + ch.trimmed(cap, -5.0) + + +# -------------------------------------------------------------------------- +# fit_rate_model +# -------------------------------------------------------------------------- + +def planted_period(cycle_count, count_rate, overhead): + return rm3100.AXES * (cycle_count / count_rate + overhead) + + +def test_fit_rate_model_is_exact_from_two_cycle_counts(tmp_path): + count_rate, overhead = 88_000.0, 40e-6 + caps = [make(tmp_path / f"{cc}.csv", rows=3000, cycle_count=cc, + dt=planted_period(cc, count_rate, overhead)) + for cc in (100, 400)] + fitted_rate, fitted_overhead = compare.fit_rate_model(caps) + assert fitted_rate == pytest.approx(count_rate, rel=1e-4) + assert fitted_overhead == pytest.approx(overhead, rel=1e-3) + + +def test_fit_rate_model_least_squares_over_three_cycle_counts(tmp_path): + count_rate, overhead = 92_000.0, 38e-6 + caps = [make(tmp_path / f"{cc}.csv", rows=3000, cycle_count=cc, + dt=planted_period(cc, count_rate, overhead)) + for cc in (100, 200, 400)] + fitted_rate, fitted_overhead = compare.fit_rate_model(caps) + assert fitted_rate == pytest.approx(count_rate, rel=1e-3) + assert fitted_overhead == pytest.approx(overhead, rel=1e-2) + + +def test_fit_rate_model_refuses_one_cycle_count(tmp_path): + """Two unknowns need two points; guessing one would look like a result.""" + caps = [make(tmp_path / f"{i}.csv", rows=1000, cycle_count=100, + dt=1 / 300.0) for i in range(3)] + with pytest.raises(ValueError, match="two distinct cycle counts"): + compare.fit_rate_model(caps) + + +def test_fit_rate_model_rejects_captures_that_do_not_follow_the_model(tmp_path): + """A higher cycle count that samples faster is not this chip.""" + caps = [make(tmp_path / "a.csv", rows=1000, cycle_count=100, dt=1 / 100.0), + make(tmp_path / "b.csv", rows=1000, cycle_count=400, dt=1 / 300.0)] + with pytest.raises(ValueError, match="not positive"): + compare.fit_rate_model(caps) + + +# -------------------------------------------------------------------------- +# gain_and_movement, axis_ratio_spread +# -------------------------------------------------------------------------- + +def vector(mean): + mean = np.asarray(mean, dtype=float) + return {"mean": mean, "field": float(np.linalg.norm(mean))} + + +def test_gain_and_movement_reads_a_pure_scale_change_exactly(): + a = vector([-40000, 16000, 12000]) + b = vector(np.array(a["mean"]) * 1.07) + scale, residual, angle = compare.gain_and_movement(a, b) + assert scale == pytest.approx(1.07) + assert residual == pytest.approx(0.0, abs=1e-12) + assert angle == pytest.approx(0.0, abs=1e-6) + + +def test_gain_and_movement_reads_a_pure_rotation_as_movement(): + a = vector([40000, 0, 0]) + b = vector([40000 * math.cos(math.radians(20)), + 40000 * math.sin(math.radians(20)), 0]) + scale, residual, angle = compare.gain_and_movement(a, b) + assert angle == pytest.approx(20.0, abs=1e-6) + assert residual == pytest.approx(math.sin(math.radians(20)), rel=1e-6) + assert residual > compare.RATIO_SPREAD_OK + + +def test_gain_and_movement_survives_an_axis_crossing_zero(): + """The case that defeats a per-axis ratio: Z changes sign between runs.""" + a = vector([-40000, 16000, -6000]) + b = vector([-39000, 16400, +8000]) + scale, residual, angle = compare.gain_and_movement(a, b) + assert math.isfinite(scale) and math.isfinite(residual) + assert 0.0 < residual < 1.0 + assert angle > 1.0 + + +def test_axis_ratio_spread_is_zero_for_a_pure_scale(): + a = vector([-40000, 16000, 12000]) + b = vector(np.array(a["mean"]) * 1.07) + spread = compare.axis_ratio_spread(a, b, b["mean"] / a["mean"]) + assert spread == pytest.approx(0.0, abs=1e-9) + + +def test_axis_ratio_spread_abstains_when_an_axis_carries_no_field(): + a = vector([40000, 100, 5]) + b = vector([41000, -90, 4]) + assert compare.axis_ratio_spread(a, b, b["mean"] / a["mean"]) is None + + +# -------------------------------------------------------------------------- +# labels, conditions, supplies, bands +# -------------------------------------------------------------------------- + +def test_split_label_reads_an_explicit_label(tmp_path): + path = tmp_path / "c.csv" + path.write_text("x") + assert compare.split_label(f"LDO/cc100={path}") == ("LDO/cc100", str(path)) + + +def test_split_label_leaves_an_unlabelled_path_alone(): + assert compare.split_label("a.csv") == (None, "a.csv") + + +def test_split_label_prefers_an_existing_path_containing_an_equals(tmp_path): + """A filename may contain '='; an existing file wins over a label reading.""" + path = tmp_path / "run=2.csv" + path.write_text("x") + assert compare.split_label(str(path)) == (None, str(path)) + + +def test_split_label_rejects_an_empty_label(): + with pytest.raises(ValueError, match="empty label"): + compare.split_label("=a.csv") + + +def test_split_label_keeps_the_labelled_reading_when_neither_exists(): + """So the error names the path the user meant, not the whole argument.""" + assert compare.split_label("LDO=missing.csv") == ("LDO", "missing.csv") + + +def test_condition_and_variant_split_on_the_first_separator(): + assert compare.condition_of("LDO/cc100") == "LDO" + assert compare.variant_of("LDO/cc100") == "cc100" + assert compare.condition_of("LDO/a/b") == "LDO" + assert compare.variant_of("LDO/a/b") == "a/b" + + +def test_a_label_without_a_separator_is_all_condition(): + assert compare.condition_of("plain.csv") == "plain.csv" + assert compare.variant_of("plain.csv") == "" + + +def test_conditions_come_back_in_command_line_order(): + records = [("Zeta/cc100", "", {}), ("Alpha/cc100", "", {}), + ("Zeta/cc400", "", {})] + assert compare.conditions_in_order(records) == ["Zeta", "Alpha"] + + +def test_parse_supply_reads_a_condition_and_volts(): + assert compare.parse_supply("LDO=3.0") == ("LDO", 3.0) + + +@pytest.mark.parametrize("bad", ["LDO", "=3.0", "LDO=", "LDO=abc", + "LDO=0", "LDO=-3", "LDO=1e9"]) +def test_parse_supply_rejects_nonsense(bad): + with pytest.raises(argparse.ArgumentTypeError): + compare.parse_supply(bad) + + +def test_parse_band_reads_a_pair(): + assert compare.parse_band("3,30") == (3.0, 30.0) + + +@pytest.mark.parametrize("bad", ["3", "3,30,300", "30,3", "-1,30", "0,30", + "3,3", "a,b"]) +def test_parse_band_rejects_nonsense(bad): + with pytest.raises(argparse.ArgumentTypeError): + compare.parse_band(bad) + + +# -------------------------------------------------------------------------- +# decimation pairing and band selection +# -------------------------------------------------------------------------- + +def entry(label, fs, cycle_count): + return (label, {"fs": fs, "cycle_count": cycle_count}) + + +def test_decimation_pairs_finds_an_integer_cycle_count_ratio(): + pairs = compare.decimation_pairs([entry("LDO/cc100", 300.0, 100), + entry("LDO/cc400", 75.0, 400)]) + assert len(pairs) == 1 + condition, (fast_label, _), (slow_label, _), k = pairs[0] + assert (condition, fast_label, slow_label, k) == ("LDO", "LDO/cc100", + "LDO/cc400", 4) + + +def test_decimation_pairs_ignores_a_non_integer_ratio(): + assert compare.decimation_pairs([entry("LDO/a", 300.0, 100), + entry("LDO/b", 75.0, 405)]) == [] + + +def test_decimation_pairs_does_not_cross_conditions(): + """Comparing a decimated LDO run to a native 3V3 one confounds the two.""" + assert compare.decimation_pairs([entry("LDO/cc100", 300.0, 100), + entry("3V3/cc400", 75.0, 400)]) == [] + + +def test_decimation_pairs_ignores_an_equal_cycle_count(): + assert compare.decimation_pairs([entry("LDO/a", 300.0, 100), + entry("LDO/b", 299.0, 100)]) == [] + + +def test_the_band_is_set_by_the_decimated_rate_not_the_slowest_capture(tmp_path): + """Decimating by 4 lands below the natively-slow rate, and that binds.""" + fast = make(tmp_path / "fast.csv", rows=2000, cycle_count=100, dt=1 / 300.0) + slow = make(tmp_path / "slow.csv", rows=2000, cycle_count=400, dt=1 / 76.0) + rates = compare.comparable_rates([("LDO/cc100", fast), ("LDO/cc400", slow)]) + assert min(rates) == pytest.approx(75.0, rel=1e-3) # 300/4, not 76 + + +# -------------------------------------------------------------------------- +# End to end +# -------------------------------------------------------------------------- + +def four_captures(tmp_path, gain=1.0): + """An interleaved two-condition, two-cycle-count session, as recorded.""" + paths = {} + for condition, scale in (("LDO", 1.0), ("3V3", gain)): + for cc, dt in ((100, 1 / 300.0), (400, 1 / 76.0)): + mean = tuple(component * scale for component in FIELD) + path = tmp_path / f"{condition}_{cc}.csv" + counts = field_counts(4096, mean, cc, 20.0, seed=cc) + write_capture(path, rows=4096, dt=dt, cycle_count=cc, counts=counts) + paths[f"{condition}/cc{cc}"] = str(path) + return paths + + +def run(monkeypatch, capsys, argv): + monkeypatch.setattr("sys.argv", ["compare.py"] + argv) + assert compare.main() == 0 + return capsys.readouterr() + + +def test_end_to_end_reports_every_section(tmp_path, monkeypatch, capsys): + paths = four_captures(tmp_path, gain=0.93) + out = run(monkeypatch, capsys, + [f"{label}={path}" for label, path in paths.items()] + + ["--supply", "LDO=3.0", "--supply", "3V3=3.3"]).out + assert "timing model" in out + assert "filtering and decimating" in out + assert "matched cycle count" in out + assert "how far the sensor moved" in out + # A planted 7% scale change on |B| should come back as one. + assert "-7.0" in out or "-6.9" in out or "-7.1" in out + + +def test_the_supply_exponent_recovers_a_planted_power_law(tmp_path, monkeypatch, + capsys): + """A gain that is exactly ratiometric must come back as V^-1.""" + volts = {"LDO": 3.0, "3V3": 3.3} + argv = [] + for condition, v in volts.items(): + for cc, dt in ((100, 1 / 300.0), (400, 1 / 76.0)): + # Gain proportional to 1/V means the reading scales as 1/V too. + mean = tuple(c * volts["LDO"] / v for c in FIELD) + path = tmp_path / f"{condition}_{cc}.csv" + write_capture(path, rows=4096, dt=dt, cycle_count=cc, + counts=field_counts(4096, mean, cc, 5.0, seed=cc)) + argv.append(f"{condition}/cc{cc}={path}") + argv += ["--supply", "LDO=3.0", "--supply", "3V3=3.3"] + out = run(monkeypatch, capsys, argv).out + assert "scales with the rail" in out + for line in out.splitlines(): + if "|B| ~ V^" in line: + power = float(line.split("|B| ~ V^")[1].split()[0]) + assert power == pytest.approx(-1.0, abs=0.02) + + +def test_the_supply_exponent_is_absent_without_rail_voltages(tmp_path, + monkeypatch, + capsys): + paths = four_captures(tmp_path, gain=0.93) + out = run(monkeypatch, capsys, + [f"{label}={path}" for label, path in paths.items()]).out + assert "scales with the rail" not in out + assert "x rail volts" not in out + + +def test_end_to_end_writes_a_figure(tmp_path, monkeypatch, capsys): + paths = four_captures(tmp_path) + output = tmp_path / "figure.png" + run(monkeypatch, capsys, + [f"{label}={path}" for label, path in paths.items()] + + ["-o", str(output)]) + assert output.exists() and output.stat().st_size > 10_000 + + +def test_a_capture_that_will_not_load_is_skipped_not_fatal(tmp_path, monkeypatch, + capsys): + paths = four_captures(tmp_path) + broken = tmp_path / "broken.csv" + broken.write_text("this is not a capture\n") + out = run(monkeypatch, capsys, + [f"{label}={path}" for label, path in paths.items()] + + [f"BAD/x={broken}"]) + assert "skipping" in out.err + assert "timing model" in out.out + + +def test_no_loadable_capture_exits_rather_than_printing_nothing(tmp_path, + monkeypatch): + broken = tmp_path / "broken.csv" + broken.write_text("nope\n") + monkeypatch.setattr("sys.argv", ["compare.py", str(broken)]) + with pytest.raises(SystemExit, match="nothing to compare"): + compare.main() + + +def test_a_supply_naming_no_capture_warns(tmp_path, monkeypatch, capsys): + paths = four_captures(tmp_path) + out = run(monkeypatch, capsys, + [f"{label}={path}" for label, path in paths.items()] + + ["--supply", "NOSUCH=3.0"]) + assert "names no capture" in out.err + + +def test_a_band_past_the_slowest_rolloff_warns(tmp_path, monkeypatch, capsys): + paths = four_captures(tmp_path) + out = run(monkeypatch, capsys, + [f"{label}={path}" for label, path in paths.items()] + + ["--band", "3,60"]) + assert "rolls off" in out.err + + +def test_trim_is_reported_and_shortens_every_capture(tmp_path, monkeypatch, + capsys): + paths = four_captures(tmp_path) + out = run(monkeypatch, capsys, + [f"{label}={path}" for label, path in paths.items()] + + ["--trim", "1"]).out + assert out.count("trimmed 1 s from each end") == 4 + + +def test_a_single_capture_still_produces_a_table(tmp_path, monkeypatch, capsys): + paths = four_captures(tmp_path) + out = run(monkeypatch, capsys, [next(iter(paths.values()))]).out + assert "capture" in out + # Nothing to contrast against, so those sections stay quiet. + assert "matched cycle count" not in out + assert "how far the sensor moved" not in out diff --git a/tests/test_logger.py b/tests/test_logger.py new file mode 100644 index 0000000..616195b --- /dev/null +++ b/tests/test_logger.py @@ -0,0 +1,422 @@ +"""The capture loop: miss counting, flagging, and the configuration warnings. + +sample_loop() is the subtlest code in the repo -- it decides how many +measurements passed unseen, and getting it wrong slips sample_index against chip +time in a way nothing downstream can detect. It is driven here by a scripted +sensor against a fake clock, so a stall of an exact size can be injected and the +resulting rows checked exactly. +""" + +import queue + +import pytest + +import capture +import logger +import rm3100 +from conftest import FakeSensor, grid_events + +PERIOD = 0.01 # 100 Hz, chosen so stalls are easy to state in periods + + +def drain(q): + """Every row the loop queued, in order.""" + rows = [] + while not q.empty(): + rows.append(q.get_nowait()) + return rows + + +def run_loop(sensor, fake_clock, duration, dt=PERIOD, maxsize=0): + """Run sample_loop against a fake sensor. Returns (rows, stats).""" + fake_clock(sensor) + q = queue.Queue(maxsize=maxsize) + stats = {"rows": 0, "missed": 0, "ambiguous": 0, "max_bracket": 0.0, + "truncated": False} + logger.sample_loop(sensor, q, duration, dt, stats) + return drain(q), stats + + +def indices(rows): + return [r[0] for r in rows] + + +def flags(rows): + return [r[5] for r in rows] + + +# -------------------------------------------------------------------------- +# The clean case +# -------------------------------------------------------------------------- + +def test_a_clean_run_loses_nothing(fake_clock): + sensor = FakeSensor(grid_events(PERIOD, 60)) + rows, stats = run_loop(sensor, fake_clock, duration=0.5) + + assert stats["missed"] == 0 + assert stats["ambiguous"] == 0 + assert all(f == "" for f in flags(rows)) + assert len(rows) > 40 + + +def test_the_sample_index_is_contiguous_from_zero(fake_clock): + sensor = FakeSensor(grid_events(PERIOD, 60)) + rows, _ = run_loop(sensor, fake_clock, duration=0.5) + assert indices(rows) == list(range(len(rows))) + + +def test_every_measurement_is_read_exactly_once(fake_clock): + """No duplicates: DRC1 clears DRDY on the results read, so a re-read cannot + return the same measurement -- and the loop must not ask for one.""" + sensor = FakeSensor(grid_events(PERIOD, 60)) + rows, _ = run_loop(sensor, fake_clock, duration=0.5) + x_values = [r[2] for r in rows] + assert len(set(x_values)) == len(x_values) + assert x_values == sorted(x_values) # grid_events counts upward + + +def test_the_timestamp_is_the_bracket_midpoint(fake_clock): + """Not the read time: the measurement completed somewhere in the bracket.""" + sensor = FakeSensor(grid_events(PERIOD, 40)) + rows, _ = run_loop(sensor, fake_clock, duration=0.3) + times = [r[1] for r in rows] + gaps = [b - a for a, b in zip(times, times[1:])] + # Resolved only to the poll interval, which is what a real bracket gives too. + assert all(g == pytest.approx(PERIOD, abs=1e-3) for g in gaps) + + +def test_the_worst_bracket_is_recorded(fake_clock): + sensor = FakeSensor(grid_events(PERIOD, 40)) + _, stats = run_loop(sensor, fake_clock, duration=0.3) + # No stall, so the bracket is one poll interval and nowhere near a period. + assert 0 < stats["max_bracket"] < PERIOD / 2 + + +def test_a_stale_drdy_at_start_is_discarded(fake_clock): + """BIST leaves DRDY high over zeroed registers; sample 0 must be fresh.""" + sensor = FakeSensor(grid_events(PERIOD, 40, start=0.0, counts=(0, 0, 0))) + sensor.pending = (0, 0, 0) # as if left set by whatever ran before + rows, _ = run_loop(sensor, fake_clock, duration=0.2) + assert rows + assert rows[0][2:5] != (0, 0, 0) or sensor.reads[0] == (0, 0, 0) + # The discard happened: one more read than rows written. + assert len(sensor.reads) == len(rows) + 1 + + +def test_duration_zero_keeps_running(fake_clock): + """0 means 'until Ctrl-C', so the deadline must be infinite, not immediate. + + Nothing inside the loop ends it, so the sensor is what stops the test -- + which is also the proof that the loop was still going. + """ + sensor = FakeSensor(grid_events(PERIOD, 200), max_polls=2000) + fake_clock(sensor) + q = queue.Queue() + stats = {"rows": 0, "missed": 0, "ambiguous": 0, "max_bracket": 0.0, + "truncated": False} + with pytest.raises(FakeSensor.Exhausted): + logger.sample_loop(sensor, q, 0.0, PERIOD, stats) + # It kept sampling until the sensor gave out, rather than returning at once. + assert q.qsize() >= 50 + + +# -------------------------------------------------------------------------- +# Losing measurements +# -------------------------------------------------------------------------- + +def test_a_two_period_stall_records_one_placeholder(fake_clock): + """Two completions inside one bracket: one was read, one was lost.""" + sensor = FakeSensor(grid_events(PERIOD, 60), stalls=[(0.2, 2.1 * PERIOD)]) + rows, stats = run_loop(sensor, fake_clock, duration=0.5) + + assert stats["missed"] == 1 + assert stats["ambiguous"] == 0 + placeholders = [r for r in rows if logger.WARN_MISSED in r[5]] + assert len(placeholders) == 1 + # A placeholder carries no data and keeps the index contiguous. + assert placeholders[0][2:5] == (None, None, None) + assert indices(rows) == list(range(len(rows))) + + +def test_a_longer_stall_records_every_lost_grid_point(fake_clock): + sensor = FakeSensor(grid_events(PERIOD, 80), stalls=[(0.2, 4.1 * PERIOD)]) + rows, stats = run_loop(sensor, fake_clock, duration=0.7) + + assert stats["missed"] == 3 + assert sum(logger.WARN_MISSED in r[5] for r in rows) == 3 + assert indices(rows) == list(range(len(rows))) + + +def test_placeholders_sit_on_the_grid_not_bunched_at_the_end(fake_clock): + """They mark where the lost measurements actually completed.""" + sensor = FakeSensor(grid_events(PERIOD, 80), stalls=[(0.2, 4.1 * PERIOD)]) + rows, _ = run_loop(sensor, fake_clock, duration=0.7) + lost = [r for r in rows if logger.WARN_MISSED in r[5]] + gaps = [b[1] - a[1] for a, b in zip(lost, lost[1:])] + assert all(g == pytest.approx(PERIOD, rel=0.05) for g in gaps) + + +def test_a_late_read_within_one_period_is_not_a_miss(fake_clock): + """The case that motivated the bracket: a host stall that lost nothing. + + An interval-based estimator would insert a spurious placeholder here and + slip the index for the rest of the run. + """ + sensor = FakeSensor(grid_events(PERIOD, 60), stalls=[(0.2, 1.1 * PERIOD)]) + rows, stats = run_loop(sensor, fake_clock, duration=0.5) + + assert stats["missed"] == 0 + assert all(logger.WARN_MISSED not in f for f in flags(rows)) + + +def test_the_stall_shows_up_in_the_worst_bracket(fake_clock): + """Losing nothing is not the same as having had margin.""" + sensor = FakeSensor(grid_events(PERIOD, 60), stalls=[(0.2, 1.1 * PERIOD)]) + _, stats = run_loop(sensor, fake_clock, duration=0.5) + assert stats["max_bracket"] > 0.7 * PERIOD + + +# -------------------------------------------------------------------------- +# Ambiguity +# -------------------------------------------------------------------------- + +def test_a_gap_near_a_half_period_is_flagged_ambiguous(fake_clock): + """1.35 periods rounds to one, so no placeholder -- but it is a guess.""" + sensor = FakeSensor(grid_events(PERIOD, 60), stalls=[(0.2, 1.35 * PERIOD)]) + rows, stats = run_loop(sensor, fake_clock, duration=0.5) + + assert stats["ambiguous"] == 1 + assert stats["missed"] == 0 + flagged = [r for r in rows if logger.WARN_AMBIGUOUS in r[5]] + assert len(flagged) == 1 + # It rides on the real sample ending the gap, which keeps its data. + assert flagged[0][2] is not None + assert logger.WARN_MISSED not in flagged[0][5] + + +def test_an_ambiguous_gap_that_did_lose_something_flags_both(fake_clock): + sensor = FakeSensor(grid_events(PERIOD, 60), stalls=[(0.2, 2.5 * PERIOD)]) + rows, stats = run_loop(sensor, fake_clock, duration=0.5) + + assert stats["ambiguous"] == 1 + assert stats["missed"] == 2 + placeholders = [r for r in rows if logger.WARN_MISSED in r[5]] + assert placeholders + assert all(logger.WARN_AMBIGUOUS in r[5] for r in placeholders) + + +def test_a_gap_landing_squarely_on_a_period_is_not_ambiguous(fake_clock): + sensor = FakeSensor(grid_events(PERIOD, 60), stalls=[(0.2, 2.02 * PERIOD)]) + _, stats = run_loop(sensor, fake_clock, duration=0.5) + assert stats["ambiguous"] == 0 + assert stats["missed"] == 1 + + +def test_the_ambiguity_margin_is_the_documented_quarter_period(): + assert logger.AMBIGUITY_MARGIN == 0.25 + + +# -------------------------------------------------------------------------- +# Backpressure +# -------------------------------------------------------------------------- + +def test_a_full_queue_truncates_rather_than_holing_the_index(fake_clock): + """A dropped row would break contiguity and make the file unloadable.""" + sensor = FakeSensor(grid_events(PERIOD, 200)) + rows, stats = run_loop(sensor, fake_clock, duration=1.5, maxsize=8) + + assert stats["truncated"] is True + assert len(rows) == 8 + assert indices(rows) == list(range(8)) + + +def test_an_ample_queue_never_truncates(fake_clock): + sensor = FakeSensor(grid_events(PERIOD, 60)) + _, stats = run_loop(sensor, fake_clock, duration=0.5, maxsize=10_000) + assert stats["truncated"] is False + + +# -------------------------------------------------------------------------- +# calibrate_period +# -------------------------------------------------------------------------- + +def test_calibrate_period_recovers_the_true_period(fake_clock): + """The nominal value is 6-9% out on real hardware; the fit has to beat that.""" + sensor = FakeSensor(grid_events(PERIOD, 400)) + fake_clock(sensor) + dt, host_cost = logger.calibrate_period(sensor, PERIOD * 1.07, seconds=1.0) + assert dt == pytest.approx(PERIOD, rel=1e-3) + assert host_cost > 0 + + +def test_calibrate_period_bootstraps_when_the_nominal_is_badly_wrong(fake_clock): + """When the cycle count governs, the TMRC table value is out by a factor. + + A threshold built on it rejects every interval, so the period is seeded from + observation instead. + """ + sensor = FakeSensor(grid_events(PERIOD, 400)) + fake_clock(sensor) + dt, _ = logger.calibrate_period(sensor, PERIOD / 8, seconds=1.0) + assert dt == pytest.approx(PERIOD, rel=1e-3) + + +def test_calibrate_period_reports_a_plausible_host_cost(fake_clock): + sensor = FakeSensor(grid_events(PERIOD, 400)) + fake_clock(sensor) + _, host_cost = logger.calibrate_period(sensor, PERIOD, seconds=1.0) + # One poll plus one read, and nothing else charged to the sample. + assert host_cost == pytest.approx( + sensor.poll_cost + sensor.read_cost, rel=0.5) + assert host_cost < PERIOD + + +def test_calibrate_period_gives_up_without_enough_clean_samples(fake_clock): + """Returning None is itself the answer: the rate is not sustainable.""" + sensor = FakeSensor(grid_events(PERIOD, 4)) + fake_clock(sensor) + assert logger.calibrate_period(sensor, PERIOD, seconds=0.5) == (None, None) + + +# -------------------------------------------------------------------------- +# End to end: sample_loop -> writer_thread -> capture.load +# -------------------------------------------------------------------------- + +def test_a_logged_run_reloads_through_capture(tmp_path, fake_clock): + """The two halves of the format contract, checked against each other.""" + sensor = FakeSensor(grid_events(PERIOD, 400), stalls=[(0.5, 2.1 * PERIOD)]) + fake_clock(sensor) + q = queue.Queue() + stats = {"rows": 0, "missed": 0, "ambiguous": 0, "max_bracket": 0.0, + "truncated": False} + logger.sample_loop(sensor, q, 2.0, PERIOD, stats) + q.put(logger._SENTINEL) + + path = tmp_path / "run.csv" + meta = { + "rm3100_capture": 1, + "nominal_rate_hz": 1 / PERIOD, + "cycle_count": 100, + "tesla_per_count": repr(rm3100.tesla_per_count(100)), + "calibrated_period_s": repr(PERIOD), + "note": "synthetic", + } + logger.writer_thread(q, str(path), meta, stats) + + cap = capture.load(str(path)) + assert len(cap.sample_index) == stats["rows"] + assert cap.missed.sum() == stats["missed"] == 1 + assert cap.meta["note"] == "synthetic" + assert cap.true_rate_hz == pytest.approx(1 / PERIOD, rel=1e-3) + + +def test_the_writer_writes_zeros_not_blanks_for_a_placeholder(tmp_path): + """x/y/z must parse as integers on every row; the flag marks the empty ones.""" + q = queue.Queue() + q.put((0, 1_700_000_000.0, 10, 20, 30, "")) + q.put((1, 1_700_000_000.01, None, None, None, logger.WARN_MISSED)) + q.put(logger._SENTINEL) + stats = {"rows": 0, "missed": 1} + + path = tmp_path / "w.csv" + logger.writer_thread(q, str(path), { + "rm3100_capture": 1, + "tesla_per_count": repr(rm3100.tesla_per_count(100))}, stats) + + lines = path.read_text().splitlines() + assert lines[-1] == f"1,1700000000.010000,0,0,0,{logger.WARN_MISSED}" + assert stats["rows"] == 2 + + +def test_the_header_is_one_line_per_key(tmp_path): + handle_path = tmp_path / "h.csv" + with open(handle_path, "w") as handle: + logger.write_header(handle, {"a": 1, "b": "two"}) + assert handle_path.read_text() == "# a: 1\n# b: two\n" + + +# -------------------------------------------------------------------------- +# Configuration warnings -- each one guards against a silent capture defect +# -------------------------------------------------------------------------- + +def warnings_for(capsys, cfg, bus_speed=750, requested_rate=None): + logger.print_plan(cfg, bus_speed, requested_rate) + return capsys.readouterr().err + + +def test_the_default_configuration_warns_about_nothing(capsys): + assert warnings_for(capsys, rm3100.plan()) == "" + + +def test_a_thin_dither_margin_is_warned_about(capsys): + err = warnings_for(capsys, rm3100.plan(cycle_count=40)) + assert "below the recommended" in err + assert "characterize.py" in err # and says how to settle it + + +def test_an_idling_sensor_is_warned_about(capsys): + err = warnings_for(capsys, rm3100.plan(cycle_count=50, tmrc=0x99)) + assert "idle" in err + + +def test_a_rate_that_could_not_be_met_is_warned_about(capsys): + err = warnings_for(capsys, rm3100.plan(rate_hz=5000), requested_rate=5000) + assert "was requested but this configuration runs" in err + + +def test_sampling_below_twice_mains_is_warned_about(capsys): + """The one error post-processing cannot undo.""" + err = warnings_for(capsys, rm3100.plan(rate_hz=32), requested_rate=32) + assert "Nyquist" in err and "mains" in err + assert "4.0" in err # 60 Hz folds to ~4 Hz at 32 Hz + + +def test_an_extrapolated_cycle_count_is_warned_about(capsys): + err = warnings_for(capsys, rm3100.plan(cycle_count=1000)) + assert "extrapolated" in err + assert str(rm3100.MAX_SPEC_CYCLE_COUNT) in err + + +def test_a_bus_that_cannot_keep_up_is_warned_about_with_a_fix(capsys): + err = warnings_for(capsys, rm3100.plan(), bus_speed=20) + assert "of every period" in err + assert "kHz would fit" in err + + +def test_the_plan_derivation_is_printed_for_checking(capsys): + logger.print_plan(rm3100.plan(), 750, None) + out = capsys.readouterr().out + for field in ("cycle count", "sensor ceiling", "TMRC", "predicted rate", + "resolution", "duty", "bus 750 kHz"): + assert field in out + + +def test_the_requested_rate_is_echoed_when_given(capsys): + logger.print_plan(rm3100.plan(rate_hz=50), 750, 50) + assert "requested" in capsys.readouterr().out + + +# -------------------------------------------------------------------------- +# Small helpers +# -------------------------------------------------------------------------- + +def test_set_thread_nice_reports_what_took_effect(): + """Lowering own priority needs no privilege, so this must work anywhere.""" + got = logger.set_thread_nice(logger.WRITER_NICE) + assert got in (logger.WRITER_NICE, None) + + +def test_raising_priority_without_privilege_returns_none(): + if logger.set_thread_nice(0) is None: + pytest.skip("thread priority is unavailable on this platform") + import os + if os.geteuid() == 0: + pytest.skip("running as root, where the raise succeeds") + assert logger.set_thread_nice(logger.SAMPLER_NICE) is None + + +def test_restore_ownership_is_a_no_op_for_a_normal_user(tmp_path): + path = tmp_path / "f.csv" + path.write_text("x") + logger.restore_ownership(str(path)) # must not raise + assert path.read_text() == "x" diff --git a/tests/test_rm3100.py b/tests/test_rm3100.py new file mode 100644 index 0000000..9dfefbb --- /dev/null +++ b/tests/test_rm3100.py @@ -0,0 +1,290 @@ +"""Driver maths and register sequences -- no hardware, no USB.""" + +import pytest + +import rm3100 +from conftest import FakeBus + + +# -------------------------------------------------------------------------- +# Gain and resolution +# -------------------------------------------------------------------------- + +@pytest.mark.parametrize("cycle_count, table_lsb_per_ut", [ + (50, 20), (100, 38), (200, 75), +]) +def test_gain_fit_reproduces_table_3_1(cycle_count, table_lsb_per_ut): + """The fit exists to reproduce Table 3-1's three points within a count.""" + fitted = rm3100.gain_lsb_per_tesla(cycle_count) / rm3100.UT_PER_TESLA + assert abs(fitted - table_lsb_per_ut) < 1.0 + + +def test_tesla_per_count_inverts_the_gain(): + for cc in (30, 100, 400, 65535): + assert rm3100.tesla_per_count(cc) * rm3100.gain_lsb_per_tesla(cc) == pytest.approx(1.0) + + +def test_higher_cycle_count_means_finer_lsb(): + steps = [rm3100.tesla_per_count(cc) for cc in (50, 100, 200, 400)] + assert steps == sorted(steps, reverse=True) + + +def test_expected_noise_follows_inverse_sqrt(): + """Table 3-1's 30/20/15 nT at 50/100/200 is a K/sqrt(cc) law.""" + assert rm3100.expected_noise_nt(50) == pytest.approx(29.4, abs=0.7) + assert rm3100.expected_noise_nt(100) == pytest.approx(20.8, abs=0.9) + assert rm3100.expected_noise_nt(200) == pytest.approx(14.7, abs=0.4) + # Quadrupling the cycle count halves the noise. + assert (rm3100.expected_noise_nt(100) + == pytest.approx(2 * rm3100.expected_noise_nt(400))) + + +# -------------------------------------------------------------------------- +# Timing model +# -------------------------------------------------------------------------- + +def test_sample_period_matches_the_documented_model(): + expected = 3 * (100 / rm3100.COUNTS_PER_SECOND + rm3100.AXIS_OVERHEAD_S) + assert rm3100.sample_period(100) == pytest.approx(expected) + # The default configuration is quoted at ~282 Hz throughout the docs. + assert 1 / rm3100.sample_period(100) == pytest.approx(282, abs=1) + + +def test_integration_time_excludes_overhead(): + assert rm3100.integration_time(100) < rm3100.sample_period(100) + assert rm3100.integration_time(100) == pytest.approx(3 * 100 / 90000.0) + + +def test_duty_rises_with_cycle_count(): + """Fixed per-axis overhead is a shrinking share as integration grows.""" + duties = [rm3100.integration_time(cc) / rm3100.sample_period(cc) + for cc in (50, 100, 200, 800)] + assert duties == sorted(duties) + assert duties[-1] < 1.0 + + +@pytest.mark.parametrize("rate", [0.5, 1, 8, 32, 128, 282, 500]) +def test_cycle_count_for_rate_inverts_sample_period(rate): + cc = rm3100.cycle_count_for_rate(rate) + achieved = 1.0 / rm3100.sample_period(cc) + # Rounding to an integer cycle count is the only error here. + assert achieved == pytest.approx(rate, rel=0.01) + + +def test_cycle_count_for_rate_clamps_to_the_register_width(): + assert rm3100.cycle_count_for_rate(1e9) == 1 + assert rm3100.cycle_count_for_rate(1e-9) == rm3100.MAX_CYCLE_COUNT + + +# -------------------------------------------------------------------------- +# plan(): every branch, since it is the only place configuration is resolved +# -------------------------------------------------------------------------- + +def test_plan_default_is_the_documented_configuration(): + cfg = rm3100.plan() + assert cfg.cycle_count == rm3100.DEFAULT_CYCLE_COUNT + assert cfg.tmrc == rm3100.TMRC_FASTEST + assert cfg.governed_by == "cycle count" + assert cfg.predicted_hz == pytest.approx(282, abs=1) + assert cfg.duty == pytest.approx(0.94, abs=0.01) + assert cfg.notes == [] + + +def test_plan_from_rate_derives_the_cycle_count(): + cfg = rm3100.plan(rate_hz=50) + assert cfg.predicted_hz == pytest.approx(50, rel=0.01) + assert cfg.governed_by == "cycle count" + assert cfg.notes == [] + + +def test_plan_from_cycle_count_ignores_the_rate_knob(): + cfg = rm3100.plan(cycle_count=400) + assert cfg.cycle_count == 400 + assert cfg.predicted_hz == pytest.approx(1 / rm3100.sample_period(400)) + + +def test_plan_raises_a_cycle_count_below_the_quantisation_floor(): + cfg = rm3100.plan(cycle_count=5) + assert cfg.cycle_count == rm3100.MIN_CYCLE_COUNT + assert any("quantisation floor" in n for n in cfg.notes) + + +def test_plan_clamps_a_cycle_count_above_the_register_width(): + cfg = rm3100.plan(cycle_count=100_000) + assert cfg.cycle_count == rm3100.MAX_CYCLE_COUNT + + +def test_plan_below_the_cycle_count_floor_hands_the_cadence_to_tmrc(): + """Under ~0.46 Hz the 16-bit register runs out and only TMRC can go slower.""" + rate = rm3100.MIN_RATE_BY_CYCLE_COUNT / 2 + cfg = rm3100.plan(rate_hz=rate) + assert cfg.cycle_count == rm3100.MAX_CYCLE_COUNT + assert cfg.governed_by == "TMRC" + assert cfg.tmrc == min(rm3100.TMRC_RATES, + key=lambda t: abs(rm3100.TMRC_RATES[t] - rate)) + assert any("below the" in n for n in cfg.notes) + + +def test_plan_warns_when_tmrc_governs_and_the_sensor_idles(): + """cc=50 can run at 534 Hz; TMRC 0x99 asks for 4.5, so it idles ~99%.""" + cfg = rm3100.plan(cycle_count=50, tmrc=0x99) + assert cfg.governed_by == "TMRC" + assert cfg.predicted_hz == pytest.approx(4.5) + assert cfg.duty < 0.02 + assert any("idle" in n for n in cfg.notes) + + +def test_plan_notes_a_rate_it_cannot_reach(): + cfg = rm3100.plan(rate_hz=5000) + assert cfg.predicted_hz < 5000 + assert any("faster than this configuration can reach" in n for n in cfg.notes) + + +def test_plan_duty_is_integration_over_the_achieved_period(): + cfg = rm3100.plan(cycle_count=200) + assert cfg.duty == pytest.approx( + rm3100.integration_time(200) * cfg.predicted_hz) + + +def test_plan_rejects_an_unknown_tmrc(): + with pytest.raises(KeyError): + rm3100.plan(tmrc=0x00) + + +# -------------------------------------------------------------------------- +# Measurement decoding -- 24-bit two's complement, big-endian +# -------------------------------------------------------------------------- + +@pytest.mark.parametrize("raw, expected", [ + (b"\x00\x00\x00" * 3, (0, 0, 0)), + (b"\x00\x00\x01" * 3, (1, 1, 1)), + (b"\xff\xff\xff" * 3, (-1, -1, -1)), + (b"\x7f\xff\xff" * 3, (8388607, 8388607, 8388607)), # positive rail + (b"\x80\x00\x00" * 3, (-8388608, -8388608, -8388608)), # negative rail +]) +def test_decode_measurements_sign_extends(raw, expected): + assert rm3100.decode_measurements(raw) == expected + + +def test_decode_measurements_keeps_axes_in_order(): + assert rm3100.decode_measurements( + bytes([0, 0, 1, 0, 0, 2, 255, 255, 255])) == (1, 2, -1) + + +@pytest.mark.parametrize("length", [0, 8, 10, 27]) +def test_decode_measurements_rejects_a_wrong_length(length): + with pytest.raises(ValueError, match="9 measurement bytes"): + rm3100.decode_measurements(bytes(length)) + + +# -------------------------------------------------------------------------- +# Register access +# -------------------------------------------------------------------------- + +def test_set_cycle_counts_writes_all_three_axes_big_endian(): + bus = FakeBus() + sensor = rm3100.RM3100(bus, 0x23) + sensor.set_cycle_counts(0x0190) + assert bus.writes == [(0x23, bytes([rm3100.REG_CCX]) + b"\x01\x90" * 3)] + assert sensor.get_cycle_counts() == (400, 400, 400) + + +@pytest.mark.parametrize("bad", [-1, 0x10000, 999999]) +def test_set_cycle_counts_rejects_out_of_range(bad): + sensor = rm3100.RM3100(FakeBus(), 0x23) + with pytest.raises(ValueError, match="outside 0..65535"): + sensor.set_cycle_counts(bad) + + +def test_set_rate_rejects_a_value_outside_the_tmrc_table(): + sensor = rm3100.RM3100(FakeBus(), 0x23) + with pytest.raises(ValueError, match="not one of"): + sensor.set_rate(0x91) + sensor.set_rate(rm3100.TMRC_FASTEST) # a valid one must not raise + + +def test_configure_accepts_a_correct_hshake_readback(): + bus = FakeBus({rm3100.REG_HSHAKE: 0x1B}) + rm3100.RM3100(bus, 0x23).configure() + assert bus.registers[rm3100.REG_HSHAKE] == rm3100.HSHAKE_DRDY_ON_READ_ONLY + + +def test_configure_ignores_the_read_only_nack_bits(): + """Bits 4-6 are NACK status, so they may read back set without meaning failure.""" + bus = FakeBus() + sensor = rm3100.RM3100(bus, 0x23) + + def read_reg(reg, count=1): + return bytes([rm3100.HSHAKE_DRDY_ON_READ_ONLY | 0x70]), 0.0, 0.0 + + sensor.read_reg = read_reg + sensor.configure() # must not raise + + +def test_configure_raises_when_drc1_did_not_take(): + """Exactly-once sampling depends on DRC1, so a bad readback is fatal.""" + sensor = rm3100.RM3100(FakeBus(), 0x23) + sensor.read_reg = lambda reg, count=1: (b"\x1b", 0.0, 0.0) + with pytest.raises(IOError, match="HSHAKE did not take"): + sensor.configure() + + +def test_read_reg_prefers_a_combined_transaction(): + bus = FakeBus({rm3100.REG_REVID: rm3100.EXPECTED_REVID}) + sensor = rm3100.RM3100(bus, 0x23) + assert sensor.revid() == rm3100.EXPECTED_REVID + # One transaction, not a separate write then read. + assert len(bus.writes) == 1 + + +def test_read_reg_falls_back_when_the_bus_cannot_do_repeated_start(): + bus = FakeBus({rm3100.REG_REVID: rm3100.EXPECTED_REVID}, combined=False) + sensor = rm3100.RM3100(bus, 0x23) + assert sensor.revid() == rm3100.EXPECTED_REVID + + +def test_read_raw_decodes_from_the_measurement_registers(): + registers = {rm3100.REG_MX + i: b for i, b in enumerate( + [0x00, 0x00, 0x01, 0xFF, 0xFF, 0xFF, 0x7F, 0xFF, 0xFF])} + sensor = rm3100.RM3100(FakeBus(registers), 0x23) + counts, mono, wall = sensor.read_raw() + assert counts == (1, -1, 8388607) + assert mono > 0 and wall > mono # the stamp travels with the data + + +def test_read_measurements_converts_with_the_active_cycle_count(): + registers = {rm3100.REG_MX + i: b for i, b in enumerate( + [0x00, 0x00, 0x26, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00])} + sensor = rm3100.RM3100(FakeBus(registers), 0x23) + sensor.cycle_count = 100 + counts, tesla = sensor.read_measurements() + assert counts[0] == 0x26 + # 38 LSB/uT at cc=100, so 0x26 counts is very close to 1 uT. + assert tesla[0] == pytest.approx(1e-6, rel=0.01) + + +def test_data_ready_reads_the_drdy_bit(): + sensor = rm3100.RM3100(FakeBus({rm3100.REG_STATUS: 0x80}), 0x23) + assert sensor.data_ready() is True + sensor = rm3100.RM3100(FakeBus({rm3100.REG_STATUS: 0x00}), 0x23) + assert sensor.data_ready() is False + + +def test_wait_for_data_times_out_without_drdy(): + sensor = rm3100.RM3100(FakeBus({rm3100.REG_STATUS: 0x00}), 0x23) + assert sensor.wait_for_data(timeout=0.01, interval=0.001) is False + + +def test_cmm_start_and_stop_write_the_documented_values(): + bus = FakeBus() + sensor = rm3100.RM3100(bus, 0x23) + sensor.start_cmm() + sensor.stop_cmm() + assert bus.writes == [ + (0x23, bytes([rm3100.REG_CMM, rm3100.CMM_ALL_AXES])), + (0x23, bytes([rm3100.REG_CMM, rm3100.CMM_OFF])), + ] + + +def test_addresses_cover_both_strap_pins(): + assert list(rm3100.RM3100.ADDRESSES) == [0x20, 0x21, 0x22, 0x23]