diff --git a/README.md b/README.md index e810fac..2ea6fab 100644 --- a/README.md +++ b/README.md @@ -147,21 +147,186 @@ system-wide `pip install`. | Flag | Default | Meaning | |---|---|---| -| `--cycle-count` | 200 | gain/resolution vs. speed; 200 → 75 LSB/µT | -| `--tmrc` | 0x96 | continuous-mode rate; 0x92 (~600 Hz) to 0x9D (~0.3 Hz) | +| `--cycle-count` | 50 | 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* | | `--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 80% of the period on the bus at cc=50 | | `--scan-only` | — | scan and exit | -CSV columns: `timestamp_iso, elapsed_s, x_raw, y_raw, z_raw, x_uT, y_uT, z_uT, -magnitude_uT`. +Exit status is non-zero if an interval was missed. The partial capture is kept +and is valid up to that point. + +### 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. + +``` +# rm3100_capture: 1 +# nominal_rate_hz: 37.0 +# tmrc: 0x96 +# cycle_count: 200 +# tesla_per_count: 1.3347570742124934e-08 +# calibrated_period_s: 0.028891621... +# i2c_address: 0x23 +# bus_speed_khz: 400 +# revid: 0x22 +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 +13,1755930857.041244,0,0,0,MISSED AMBIGUOUS +21,1755930857.272108,-2761,769,377,AMBIGUOUS +``` + +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. + +| flag | meaning | +|---|---| +| `MISSED` | placeholder row; the measurement was never read, so it has no data. The row exists to keep `sample_index` contiguous | +| `AMBIGUOUS` | the gap ending at this row was of uncertain length, so the index may have slipped from here on | + +**The two are independent.** A gap measuring 1.35 periods rounds to one, so no +placeholder is written — yet it sits far enough from an integer to distrust, and +that case carries `AMBIGUOUS` on the *real* sample ending the gap, which keeps +its data. A placeholder inside an uncertain gap carries both. + +The flag is a confidence measure, not a claim that 1.35 and 1.65 are equally +likely — 1.35 probably is one period and 1.65 probably two. Deciding which needs +neighbouring timestamps and assumptions and can still be wrong, which is exactly +why these rows are flagged rather than silently resolved. + +`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. + +All times are unix epoch seconds; no ISO-8601 anywhere in the data. ```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 ``` +### One-step capture and analysis + +```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. + +**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. + +**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. +- *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. + +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. + +**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: + +``` +elapsed_nominal = N * dt_nominal # uniform, but ~6% 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%. + ### Sample rate **Cycle count, not TMRC, sets the ceiling.** Measured with `--bus-speed 400`: @@ -202,11 +367,16 @@ silently losing data. | 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 | +| `logger.py` | CLI: scan → identify → configure → log, aborting on a missed interval | +| `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 | -| `diagnose.py` | walks USB → bus → identity → registers → BIST → live read, stopping at the first failure | +| `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 | +| `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 | The adapter/sensor split is deliberate: `ch347.py` is a general I2C master usable with any device, and `rm3100.py` needs only a bus object exposing @@ -270,8 +440,12 @@ Symptoms map cleanly onto causes, so work down this list. | Symptom | Cause | |---|---| -| `Cannot claim CH347 interface 2: Access denied` | udev rule missing — run `./setup.sh` | +| `Cannot claim CH347 interface 2: Access denied` | udev rule missing — run `./setup.sh`. Do **not** use `sudo`: it works as your normal user, and running as root leaves root-owned capture files | | `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 | @@ -318,7 +492,424 @@ 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. -### 60 s capture at maximum rate +### 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: + +```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 +``` + +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 50 at 750 kHz**, which is what the defaults now do. That runs the +sensor at its own ceiling of ~534 Hz, with TMRC held fast so the cycle count +governs and the duty cycle stays near 100%. + +```bash +./.venv/bin/python logger.py --duration 60 # cc=50, 750 kHz, ~534 Hz +``` + +Against cycle count 200 this costs **4.4%** in post-filter noise — 89% duty +against 97%, and nothing else — while giving **3.7x the spectrum**, 267 Hz of +Nyquist against 73 Hz. That trade is worth taking because **aliased interference +cannot be filtered out afterwards at any cycle count**, so seeing it beats a few +percent of noise. At 534 Hz mains and its first four harmonics all sit in band +and can be notched; at 145 Hz only the 60 Hz fundamental does, and the rest fold +irreversibly onto signal. + +750 kHz is not optional at this cycle count. The period is 1.873 ms, and the +irreducible I2C traffic per sample is: + +| bus | traffic | share of period | +|---|---|---| +| 100 kHz | 1.500 ms | **80%** | +| 400 kHz | 0.375 ms | 20% | +| 750 kHz | 0.200 ms | **11%** | + +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 100 or 200 if dither is thin.** The one assumption +worth checking is that cc=50 stays dithered: its intrinsic noise is 0.58 LSB, +and simulation puts the failure point below ~0.2 LSB, so there is roughly 3x of +margin. `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 100 (1.17 LSB) or 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% | +| 100 | 281 Hz | 141 Hz | 26.17 | 0.80 | 93.7% | +| 200 | 145 Hz | 73 Hz | 13.35 | 1.10 | 97.0% | + +### 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 @@ -359,6 +950,18 @@ Three things worth noting: - [x] Register write/read verified — cycle counts written and read back - [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] Verified on hardware: exact grid, −6.11% drift recovered, 0.017% run-to-run +- [x] 3.0 V LDO built and compared — noise indistinguishable once normalised +- [ ] 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 -Initial communication is complete. A proper Python API is planned; these modules -are the initial-communication milestone, deliberately kept simple. +> **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. + +A proper Python API is planned; these modules are the initial-communication +milestone, deliberately kept simple. diff --git a/capture.py b/capture.py new file mode 100644 index 0000000..c7591a1 --- /dev/null +++ b/capture.py @@ -0,0 +1,243 @@ +"""Read an RM3100 capture and reconstruct everything the file does not store. + +A capture holds only irreducible facts -- sample count, host clock, raw counts -- +plus a header of configuration. This module turns that back into field values and +a time base, and is the single place either is derived. + +Two time bases come out of `load()`: + + elapsed_nominal N * dt_nominal, using the TMRC table rate + elapsed N * dt_true, using the rate actually observed + +The chip's RC oscillator is accurate to about +/-7% (manual section 5.2.1) and +our unit runs ~6% slow, so `elapsed_nominal` is uniform but wrongly scaled. +`dt_true` comes from regressing the host clock on the sample index: the chip +supplies short-term regularity, the host supplies long-term rate calibration. +""" + +import csv + +import numpy as np + +import rm3100 + +# Flags logger.py writes in the `warning` column, space separated. Duplicated +# rather than imported so that reading a capture never pulls in the USB stack. +WARN_MISSED = "MISSED" # placeholder row: no data, keeps the index contiguous +WARN_AMBIGUOUS = "AMBIGUOUS" # gap ending here was of uncertain length + + +class CaptureError(Exception): + pass + + +class Capture: + """A loaded capture: raw counts, nanotesla axes, and two time bases.""" + + def __init__(self, path, meta, sample_index, system_time, counts, + missed=None, ambiguous=None): + self.path = path + self.missed = (np.zeros(len(sample_index), dtype=bool) + if missed is None else missed) + # A subset of missed: gaps whose length could not be counted + # confidently, so sample_index may have slipped across them. + self.ambiguous = (np.zeros(len(sample_index), dtype=bool) + if ambiguous is None else ambiguous) + self.meta = meta + self.sample_index = sample_index + self.system_time = system_time + self.counts = counts + + self.cycle_count = int(meta["cycle_count"]) + self.tesla_per_count = float(meta["tesla_per_count"]) + self.lsb_nt = self.tesla_per_count * rm3100.NT_PER_TESLA + self.nominal_rate_hz = float(meta["nominal_rate_hz"]) + self.dt_nominal = 1.0 / self.nominal_rate_hz + + 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) + 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.elapsed_nominal = sample_index * self.dt_nominal + self.elapsed = sample_index * self.dt_true + + # Two different things live in the timing error, and conflating them is + # misleading. Read latency is local: the scatter of one read interval + # about the next. Fit residual is global: how far the whole capture + # departs from a single straight line, which on a long run is dominated + # by the oscillator's rate drifting with temperature, not by the host. + self.read_jitter = float(np.diff(system_time).std()) if len(system_time) > 1 else 0.0 + self.residual_sd = float(self.residuals.std()) + self.residual_span = float(np.ptp(self.residuals)) + # One sample period of accumulated error means the single slope is no + # longer describing the capture. + self.drift_limited = self.residual_sd > self.dt_true + + @property + def true_rate_hz(self): + return 1.0 / self.dt_true + + @property + def duration(self): + return float(self.elapsed[-1] - self.elapsed[0]) + + def axes(self): + """(key, array) for the three axes plus the derived total.""" + return [("x", self.x), ("y", self.y), ("z", self.z), + ("total", self.total)] + + def restrict(self, start=None, end=None): + """Return a new Capture covering a window of drift-corrected seconds. + + Re-fits on the subset, so a window's rate is its own rather than + inherited -- which is what makes it usable for spotting drift across a + long run by comparing windows. + """ + keep = np.ones(len(self.sample_index), dtype=bool) + if start is not None: + keep &= self.elapsed >= start + if end is not None: + keep &= self.elapsed <= end + if keep.sum() < 64: + raise CaptureError( + f"{self.path}: window [{start}, {end}] leaves " + f"{int(keep.sum())} samples, too few to analyse") + return Capture(self.path, self.meta, self.sample_index[keep], + self.system_time[keep], + {a: self.counts[a][keep] for a in "xyz"}, + self.missed[keep], self.ambiguous[keep]) + + def summary(self): + lines = [ + f"{self.path}: {len(self.sample_index):,} samples over " + f"{self.duration:.2f} s", + f" cycle count {self.cycle_count}, 1 LSB = {self.lsb_nt:.2f} nT", + f" rate {self.true_rate_hz:.3f} Hz measured vs " + f"{self.nominal_rate_hz:g} Hz nominal ({self.rate_error * 100:+.2f}%)", + f" read jitter {self.read_jitter * 1e3:.3f} ms sd | " + f"fit residual {self.residual_sd * 1e3:.1f} ms sd, " + f"{self.residual_span:.2f} s span", + f" Nyquist {self.true_rate_hz / 2:.2f} Hz", + ] + if self.drift_limited: + lines += ["", + f" WARNING: the single-rate model does not fit this capture " + f"(residual sd {self.residual_sd * 1e3:.0f} ms", + f" against a {self.dt_true * 1e3:.1f} ms period)."] + if self.missed.any(): + # Each lost measurement is an independent chance to insert one + # placeholder too many or too few, and the error accumulates in + # sample_index. That is the likelier cause here than thermal drift. + lines += [ + f" With {int(self.missed.sum()):,} lost measurement(s) the likely cause is " + "miscounted placeholders:", + " how many grid points passed unseen can only be estimated, so the index", + " slips by about one per miss. The grid is exact only in a loss-free " + "capture.", + " Re-record at a lower rate rather than trusting this one for " + "spectral work.", + ] + else: + lines += [ + " No measurements were lost, so this is the chip's oscillator drifting --", + " expect thermal variation over a long run. Frequencies are scaled by an", + " average rate and will be smeared; analyse shorter windows", + " (capture.restrict) for spectral work.", + ] + if self.missed.any(): + n = int(self.missed.sum()) + lines.insert(1, f" {n:,} lost measurement(s) " + f"({n / len(self.missed) * 100:.3f}%), interpolated") + if self.ambiguous.any(): + lines.insert(1, f" {int(self.ambiguous.sum()):,} AMBIGUOUS gap(s) " + "-- length uncertain, index may have slipped") + if self.meta.get("note"): + lines.insert(1, f" note: {self.meta['note']}") + return "\n".join(lines) + + +def _parse_header(handle): + """Consume leading '# key: value' lines, leaving the reader at the CSV.""" + meta = {} + while True: + position = handle.tell() + line = handle.readline() + if not line: + break + if not line.startswith("#"): + handle.seek(position) + break + key, _, value = line[1:].partition(":") + if value: + meta[key.strip()] = value.strip() + return meta + + +def load(path): + """Load a capture written by logger.py.""" + with open(path, newline="") as handle: + meta = _parse_header(handle) + rows = list(csv.DictReader(handle)) + + if "rm3100_capture" not in meta: + raise CaptureError( + f"{path}: no capture header found. Files written before the header " + "format was introduced cannot be read -- re-record them.") + missing = {"cycle_count", "tesla_per_count", "nominal_rate_hz"} - meta.keys() + if missing: + raise CaptureError(f"{path}: header missing {', '.join(sorted(missing))}") + if len(rows) < 64: + raise CaptureError(f"{path}: only {len(rows)} samples, too few to analyse") + + # logger.py writes "MISSED" in the count fields for a measurement it could + # not read, keeping sample_index contiguous so the chip-time grid stays + # valid across the gap. + if "warning" not in rows[0]: + raise CaptureError( + f"{path}: no 'warning' column. Captures predating it cannot be " + "read -- re-record them.") + flags = [set(r["warning"].split()) for r in rows] + missed = np.array([WARN_MISSED in f for f in flags]) + # Independent of missed: a gap can be uncertain yet round to zero losses, + # in which case the flag rides on the real sample that ends it. + ambiguous = np.array([WARN_AMBIGUOUS in f for f in flags]) + if missed.all(): + 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. + gaps = np.diff(sample_index) + if np.any(gaps != 1): + bad = int(sample_index[np.argmax(gaps != 1)]) + raise CaptureError( + f"{path}: sample_index is not contiguous (breaks after {bad}). " + "The chip grid is only valid for an unbroken index.") + + # Lost measurements carry no data. They are linearly interpolated so the + # uniform grid the spectra depend on is preserved, and counted so the + # substitution is never silent. + counts = {} + good = ~missed + for a in "xyz": + # Every row parses as an integer -- placeholders carry zeros, and the + # warning column is what marks them as having no data. + v = np.array([int(r[f"{a}_raw"]) for r in rows], dtype=np.float64) + if missed.any(): + v[missed] = np.interp(np.flatnonzero(missed), np.flatnonzero(good), + v[good]) + counts[a] = v + + return Capture( + path, meta, sample_index, + np.array([float(r["system_time_unix"]) for r in rows]), + counts, missed, ambiguous, + ) diff --git a/ch347.py b/ch347.py index b7d247b..3604435 100644 --- a/ch347.py +++ b/ch347.py @@ -13,6 +13,8 @@ aystarik/ch347-i2c-spi-gpio Linux driver, whose USB id table entry USB_DEVICE_INTERFACE_NUMBER(0x1a86, 0x55db, 0x02) matches this device. """ +import time + import usb.core import usb.util @@ -34,7 +36,11 @@ CMD_OUT = 0x80 # write; low 6 bits are the byte count CMD_IN = 0xC0 # read; low 6 bits are the byte count CMD_SET = 0x60 # set bus speed; low bits select the rate -# Bus speeds, as the low nibble of CMD_SET. +# Bus speeds, as the low nibble of CMD_SET. The fastest is the default: bus +# speed does not set the sample rate -- the sensor does -- but it does set how +# long a read takes, and so how tightly an event on the far side of the bus can +# be timestamped. Measured host cost per sample is bus time plus ~0.6 ms of +# fixed USB round-trip latency: 2.13 ms at 100 kHz, 0.92 at 400, 0.79 at 750. SPEED_20KHZ = 0 SPEED_100KHZ = 1 SPEED_400KHZ = 2 @@ -54,7 +60,7 @@ class CH347I2C: followed by the payload. """ - def __init__(self, speed=SPEED_100KHZ): + def __init__(self, speed=SPEED_750KHZ): self._dev = usb.core.find(idVendor=VENDOR_ID, idProduct=PRODUCT_ID) if self._dev is None: raise IOError( @@ -90,21 +96,32 @@ class CH347I2C: return False def _xfer(self, out_bytes, in_len): - """Send one command stream, then read in_len status/data bytes.""" + """Send one command stream, then read in_len status/data bytes. + + Returns (data, mono, wall), the clocks read the instant the reply + lands -- before any parsing or unwinding of the call stack. A caller + timing an event on the far side of the bus (DRDY going high, say) wants + that moment, not one several Python frames later: the frames add both a + systematic lag and jitter from whatever the interpreter does between. + + The stamp travels with its own reply rather than being left on the + instance, so it cannot go stale or be picked up by the wrong call. + """ written = self._dev.write(EP_OUT, bytes(out_bytes), TIMEOUT_MS) if written != len(out_bytes): raise IOError( f"Short USB write to CH347: sent {written} of {len(out_bytes)} bytes" ) if in_len == 0: - return b"" + return b"", time.monotonic(), time.time() reply = bytes(self._dev.read(EP_IN, in_len, TIMEOUT_MS)) + mono, wall = time.monotonic(), time.time() if len(reply) != in_len: raise IOError( f"Short USB read from CH347: got {len(reply)} of {in_len} bytes" ) - return reply + return reply, mono, wall def set_speed(self, speed): """Select the I2C clock rate (one of the SPEED_* constants).""" @@ -126,7 +143,7 @@ class CH347I2C: packet += [CMD_STO, CMD_END] # One status byte per clocked-out byte: the address plus the payload. - reply = self._xfer(packet, len(data) + 1) + reply, _, _ = self._xfer(packet, len(data) + 1) if reply[0] != ACK: raise IOError(f"No ACK from I2C address 0x{addr:02x} on write") if any(byte != ACK for byte in reply[1:]): @@ -136,7 +153,7 @@ class CH347I2C: ) def read(self, addr, count): - """Read count bytes from a 7-bit address. Raises IOError on NACK.""" + """Read count bytes from an address. Returns (data, mono, wall).""" if not 1 <= count <= MAX_XFER: raise ValueError(f"Read of {count} bytes outside 1..{MAX_XFER}") @@ -147,16 +164,58 @@ class CH347I2C: packet.append(CMD_IN | (count - 1)) packet += [CMD_IN, CMD_STO, CMD_END] - reply = self._xfer(packet, count + 1) + reply, mono, wall = self._xfer(packet, count + 1) if reply[0] != ACK: raise IOError(f"No ACK from I2C address 0x{addr:02x} on read") - return reply[1:] + return reply[1:], mono, wall + + def write_read(self, addr, data, count): + """Write data then read count bytes in one transaction. + + Returns (data, mono, wall) -- see _xfer() on why the stamp is returned + rather than stored. + + Uses a repeated START rather than STOP-then-START, so the whole + exchange is a single USB round trip instead of two. Measured at + 400 kHz this halves a one-byte register read (0.65 -> 0.35 ms), which + matters because the host round trip, not the bus, sets the sample-rate + ceiling. + + The RM3100's manual draws its register reads with a STOP between the + pointer write and the read (sections 4.5.2, 5.8.4), but the part + accepts a repeated START -- verified against REVID and the measurement + registers. + """ + data = bytes(data) + if not data: + raise ValueError("write_read needs at least one byte to write") + if len(data) > MAX_XFER - 1: + raise ValueError(f"Write of {len(data)} bytes exceeds the CH347 limit") + if not 1 <= count <= MAX_XFER: + raise ValueError(f"Read of {count} bytes outside 1..{MAX_XFER}") + + packet = [CMD_STREAM, CMD_STA, CMD_OUT | (len(data) + 1), addr << 1] + packet += data + packet += [CMD_STA, CMD_OUT | 1, (addr << 1) | 1] + if count > 1: + packet.append(CMD_IN | (count - 1)) + packet += [CMD_IN, CMD_STO, CMD_END] + + # One status byte per clocked-out byte: the write address, the payload, + # and the read address -- then the payload itself. + acks = len(data) + 2 + reply, mono, wall = self._xfer(packet, acks + count) + if any(byte != ACK for byte in reply[:acks]): + raise IOError( + f"I2C address 0x{addr:02x} NACKed during combined transfer " + f"(status {reply[:acks].hex(' ')})") + return reply[acks:], mono, wall def probe(self, addr): """Return True if a device ACKs its address. A NACK is not an error.""" packet = [CMD_STREAM, CMD_STA, CMD_OUT | 1, addr << 1, CMD_STO, CMD_END] try: - return self._xfer(packet, 1)[0] == ACK + return self._xfer(packet, 1)[0][0] == ACK except IOError: return False diff --git a/characterize.py b/characterize.py index e6b5f32..1388a3e 100644 --- a/characterize.py +++ b/characterize.py @@ -12,18 +12,17 @@ Produces a four-panel figure and a text summary: over -- the honest measure of a noise floor. Residual distribution after removing a slow trend, so a non-Gaussian tail or quantisation shows up. - Sample interval whether the timing supports spectral analysis - at all. + Read latency host-side diagnostics only -- the measurement + grid itself is uniform regardless. -Timing caveat: in continuous measurement mode the sensor samples on its own -internal schedule, so the true sample instants are near-uniform even when our -reads are jittery. The spectral estimates assume uniform spacing at the mean -observed rate. That assumption holds only if no samples were missed or read -twice -- which is exactly what the sample-interval panel is there to check. +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. """ import argparse -import csv import sys import matplotlib @@ -31,6 +30,8 @@ matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np +import capture + SURFACE = "#fcfcfb" TEXT_PRIMARY = "#0b0b0b" TEXT_SECONDARY = "#52514e" @@ -40,27 +41,16 @@ REFERENCE = "#8a8880" # Categorical slots 1-3; three peer axes, validated all-pairs in light mode. AXES = [("x", "X", "#2a78d6"), ("y", "Y", "#eb6834"), ("z", "Z", "#1baf7a")] +# Total field is derived from the three axes rather than a fourth peer, so it +# takes neutral ink instead of a competing hue -- which also keeps the +# categorical set at the three slots that validate. +TOTAL = ("total", "Total", TEXT_PRIMARY) +SERIES = AXES + [TOTAL] + # Table 3-1: "Noise Density @ Max. Single-Axis Sample Rate". SPEC_ASD_NT = 1.2 -def load(path): - t, cols, gains = [], {a: [] for a, _, _ in AXES}, [] - with open(path, newline="") as fh: - for row in csv.DictReader(fh): - t.append(float(row["elapsed_s"])) - for a, _, _ in AXES: - ut = float(row[f"{a}_uT"]) - cols[a].append(ut) - if abs(ut) > 1.0: - gains.append(float(row[f"{a}_raw"]) / ut) - if len(t) < 64: - sys.exit(f"{path} has too few samples to characterise") - return (np.array(t), - {a: np.array(v) * 1000.0 for a, v in cols.items()}, # work in nT - float(np.median(gains))) - - def welch_asd(v, fs, nperseg=4096): """Amplitude spectral density in units/sqrt(Hz) via Welch's method.""" nperseg = min(nperseg, len(v) // 4 * 2 or len(v)) @@ -115,32 +105,39 @@ def main(): help="ignore samples after this elapsed time (s)") args = ap.parse_args() - t, data, gain = load(args.csv) - keep = t >= args.start - if args.end is not None: - keep &= t <= args.end - t, data = t[keep], {a: v[keep] for a, v in data.items()} - if len(t) < 64: - sys.exit("selected window has too few samples") + try: + cap = capture.load(args.csv) + if args.start or args.end is not None: + cap = cap.restrict(args.start or None, args.end) + except (OSError, capture.CaptureError) as exc: + sys.exit(str(exc)) - duration = t[-1] - t[0] - fs = (len(t) - 1) / duration - intervals = np.diff(t) - lsb_nt = 1000.0 / gain + # Total field is derived from the three axes; its noise is dominated by + # whichever axis carries most of the field, weighted by direction cosines, + # so it tracks that axis rather than being an independent measurement. + data = dict(cap.axes()) + t = cap.elapsed + # Rate measured against the host clock, not the nominal table value: a 6% + # error would put every spectral feature 6% off. + fs = cap.true_rate_hz + lsb_nt = cap.lsb_nt + # Host read latency, kept apart from the sample grid it does not affect. + latency = np.diff(cap.system_time) # local read spacing, not the fit residual - print(f"{args.csv}: {len(t):,} samples over {duration:.2f} s") - print(f" mean rate {fs:.1f} Hz interval {intervals.mean()*1000:.3f} ms " - f"+/- {intervals.std()*1000:.3f} ms max {intervals.max()*1000:.1f} ms") - print(f" gain {gain:.2f} LSB/uT -> 1 LSB = {lsb_nt:.1f} nT") - print(f" Nyquist {fs/2:.1f} Hz\n") - print("axis sd p2p sd/LSB white-noise ASD median ASD") + # Normalising by |B| makes noise comparable across captures whose scale + # differs -- a gain change carries the noise with it, so absolute sd alone + # will read as a noise difference when only the scale moved. + field = data["total"].mean() + + print(cap.summary() + "\n") + print("series sd sd/|B| p2p sd/LSB white-noise ASD median ASD") fig, axs = plt.subplots(2, 2, figsize=(13.5, 9), dpi=150) fig.patch.set_facecolor(SURFACE) for ax in axs.flat: ax.set_facecolor(SURFACE) - for key, label, color in AXES: + for key, label, color in SERIES: v = data[key] sd = v.std() # A flat (white) spectrum of this sd would sit at this level. @@ -159,7 +156,7 @@ def main(): axs[1, 0].hist(resid, bins=120, histtype="step", linewidth=1.4, color=color, label=label, density=True) - print(f"{label:4s} {sd:8.1f} {v.max()-v.min():9.1f} nT " + 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") @@ -186,15 +183,18 @@ def main(): a.set_xlabel("nT"); a.set_ylabel("density") a = axs[1, 1] - a.hist(intervals * 1000, bins=120, color=REFERENCE) + a.hist(latency * 1000, bins=120, color=REFERENCE) a.set_yscale("log") - a.axvline(1000 / fs, color=TEXT_PRIMARY, linestyle="--", linewidth=1.2) - a.annotate(f"mean {1000/fs:.2f} ms", xy=(1000 / fs, 1), + a.axvline(cap.dt_true * 1000, color=TEXT_PRIMARY, linestyle="--", linewidth=1.2) + a.annotate(f"grid {cap.dt_true * 1000:.2f} ms", xy=(cap.dt_true * 1000, 1), xytext=(6, 0), textcoords="offset points", color=TEXT_PRIMARY, fontsize=9) - a.set_title("Sample interval", loc="left", + a.set_title("Host read latency", loc="left", color=TEXT_PRIMARY, fontsize=12, fontweight="bold", pad=8) - a.set_xlabel("ms"); a.set_ylabel("count") + a.annotate("spread is host-side only;\nthe sample grid is exact", + xy=(0.98, 0.94), xycoords="axes fraction", ha="right", va="top", + color=TEXT_SECONDARY, fontsize=9) + a.set_xlabel("interval between reads (ms)"); a.set_ylabel("count") for ax in axs.flat: ax.grid(True, which="both", color=GRID, linewidth=0.7) @@ -206,15 +206,19 @@ def main(): ax.tick_params(colors=TEXT_SECONDARY, labelsize=9, length=0) ax.xaxis.label.set_color(TEXT_SECONDARY) ax.yaxis.label.set_color(TEXT_SECONDARY) - for ax in (axs[0, 0], axs[0, 1], axs[1, 0]): - ax.legend(frameon=False, fontsize=9, labelcolor=TEXT_SECONDARY) + # Explicit placement: "best" puts the ASD legend on top of the spec-line + # annotation in the lower left. + for ax, loc in ((axs[0, 0], "upper right"), (axs[0, 1], "upper right"), + (axs[1, 0], "upper right")): + ax.legend(frameon=False, fontsize=9, labelcolor=TEXT_SECONDARY, loc=loc) fig.suptitle("RM3100 noise characterisation", color=TEXT_PRIMARY, fontsize=15, fontweight="bold", y=0.985) fig.text(0.5, 0.945, - f"{len(t):,} samples, {duration:.1f} s at {fs:.0f} Hz, " - f"1 LSB = {lsb_nt:.1f} nT. Spectra assume uniform sampling at the " - f"mean rate (see interval panel).", + 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.", color=TEXT_SECONDARY, fontsize=10, ha="center") fig.tight_layout(rect=[0, 0, 1, 0.935]) diff --git a/compare.py b/compare.py new file mode 100644 index 0000000..8e11892 --- /dev/null +++ b/compare.py @@ -0,0 +1,144 @@ +#!/usr/bin/env python3 +"""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 + +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. + +The second trap is attributing a scale change to the variable under test when +the sensor simply moved. Two diagnostics separate them: + + 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. + +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. +""" + +import argparse +import itertools +import sys +from collections import defaultdict + +import numpy as np + +import capture +import characterize as ch + +# 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 + + +def summarise(cap): + 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} + 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) + out[key] = { + "sd": v.std(), + "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"), + } + taus, devs = ch.allan_deviation(d["total"], fs) + i = int(np.argmin(devs)) + out["allan"] = (float(devs[i]), float(taus[i])) + return out + + +def label_for(cap, path, group_by): + if group_by == "note": + return cap.meta.get("note", "(no note)") + return path + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("csv", nargs="+") + ap.add_argument("--group", choices=["note", "file"], default="file", + help="group captures by header note, averaging repeats " + "(default: %(default)s)") + args = ap.parse_args() + + recs = [] + for path in args.csv: + try: + cap = capture.load(path) + except (OSError, capture.CaptureError) as exc: + print(f"skipping {path}: {exc}", file=sys.stderr) + continue + recs.append((label_for(cap, path, args.group), path, summarise(cap))) + if len(recs) < 1: + 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") + + # 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 "")) + + if len(recs) < 2: + return 0 + + 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.") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/diagnose.py b/diagnose-comms.py similarity index 75% rename from diagnose.py rename to diagnose-comms.py index e95cd44..eefcb9f 100644 --- a/diagnose.py +++ b/diagnose-comms.py @@ -6,7 +6,7 @@ broken, so a failure points at a specific wire rather than "it doesn't work": adapter -> I2C bus -> sensor identity -> registers -> analog section -Usage: ./.venv/bin/python diagnose.py +Usage: ./.venv/bin/python diagnose-comms.py """ import sys @@ -19,6 +19,9 @@ 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 + def check(label, ok, detail=""): print(f" [{'OK ' if ok else 'FAIL'}] {label}" + (f" -- {detail}" if detail else "")) @@ -33,7 +36,8 @@ def main(): check("open CH347", False, str(exc)) print("\n-> Run ./setup.sh, or check the adapter is plugged in and in Mode 1.") return 1 - check("open CH347 and set I2C speed", True, "interface 2 claimed, 100 kHz") + check("open CH347 and set I2C speed", True, + f"interface {ch347.INTERFACE} claimed, {BUS_KHZ} kHz") with bus: print("\n2. I2C bus") @@ -77,7 +81,7 @@ def main(): sensor.write_reg(rm3100.REG_BIST, [BIST_RUN]) sensor.write_reg(rm3100.REG_POLL, [rm3100.POLL_ALL_AXES]) time.sleep(0.5) - bist = sensor.read_reg(rm3100.REG_BIST)[0] + bist = sensor.read_reg(rm3100.REG_BIST)[0][0] sensor.write_reg(rm3100.REG_BIST, [0x00]) axes = {"X": (bist >> 4) & 1, "Y": (bist >> 5) & 1, "Z": (bist >> 6) & 1} @@ -98,18 +102,31 @@ def main(): print("\n6. Live measurement") sensor.configure() + # BIST zeroes the measurement registers and leaves DRDY set. With + # DRC1=1 only a results read clears it, so discard one -- otherwise + # wait_for_data() returns instantly on the stale flag and reports zeros. + sensor.read_raw() sensor.set_rate(0x96) sensor.start_cmm() try: ready = sensor.wait_for_data(timeout=2.0) if not check("DRDY asserts", ready, "" if ready else "no data within 2 s"): return 1 - counts, ut = sensor.read_measurements() - magnitude = sum(v * v for v in ut) ** 0.5 - check("field magnitude is plausible", 25 <= magnitude <= 65, - f"{magnitude:.1f} uT (Earth's field is 25-65)") + counts, tesla = sensor.read_measurements() + magnitude = sum(v * v for v in tesla) ** 0.5 + nt = [v * rm3100.NT_PER_TESLA for v in tesla] + plausible = check( + "field magnitude is plausible", 25e-6 <= magnitude <= 65e-6, + f"{magnitude * rm3100.NT_PER_TESLA:,.0f} nT " + "(Earth's field is 25,000-65,000)") print(f" raw {counts}") - print(f" X {ut[0]:+.3f} Y {ut[1]:+.3f} Z {ut[2]:+.3f} uT") + print(f" X {nt[0]:+,.0f} Y {nt[1]:+,.0f} Z {nt[2]:+,.0f} nT") + if not plausible: + print("\n-> Communication works but the numbers are wrong. A reading of\n" + " exactly zero means the results registers were read without a\n" + " fresh measurement behind them; anything else suggests the gain\n" + " or the int24 decode.") + return 1 finally: sensor.stop_cmm() diff --git a/logger.py b/logger.py index 6bc0b02..415e535 100755 --- a/logger.py +++ b/logger.py @@ -3,33 +3,34 @@ Run ./setup.sh first to install the udev rule and create the virtualenv, then: - ./.venv/bin/python logger.py --duration 10 + ./.venv/bin/python logger.py --duration 60 ./.venv/bin/python logger.py --scan-only -The sampling loop is kept clear of everything that is not I2C traffic: raw -counts go onto a queue, and a writer thread does the unit conversion, CSV -formatting and console display. At the rates the sensor can reach (~540 Hz at -cycle count 50) that formatting work is otherwise the bottleneck, not the bus. +The chip samples on its own internal schedule (manual sections 5.7.2, 5.8.2), so +`sample_index` is a grid coordinate, not a count of reads. A late read returns the +*newer* measurement rather than a delayed one, so an unnoticed miss would silently +compress the time axis. Misses are therefore detected and written as explicit +placeholder rows, keeping the index contiguous and the gap visible; the run +continues and reports the total at the end. + +Only irreducible facts are written: the sample count, the host clock, and the raw +counts. Chip time, tesla and magnitude are all reconstructed on load from the +header -- see capture.py. """ import argparse import csv -import math import queue import sys import threading import time -from datetime import datetime, timedelta, timezone +from datetime import datetime import ch347 import rm3100 -CSV_FIELDS = [ - "timestamp_iso", "elapsed_s", - "x_raw", "y_raw", "z_raw", - "x_uT", "y_uT", "z_uT", - "magnitude_uT", -] +CSV_FIELDS = ["sample_index", "system_time_unix", "x_raw", "y_raw", "z_raw", + "warning"] BUS_SPEEDS = { 20: ch347.SPEED_20KHZ, @@ -41,6 +42,35 @@ BUS_SPEEDS = { CONSOLE_REFRESH_S = 0.05 FLUSH_INTERVAL_S = 0.5 +# The sampling loop releases the GIL inside each USB transfer, then has to take +# it back. Python's default 5 ms switch interval means the writer thread can +# hold it for that long, stalling the sampler and inflating the DRDY bracket -- +# measured at ~10 ms against a 6.67 ms budget at 150 Hz, while the same loop +# without a writer thread peaked at 3.1 ms. Handing off 10x more often costs +# negligible throughput (the writer is not CPU-bound) and keeps the stall well +# inside the bracket budget. +GIL_SWITCH_INTERVAL_S = 0.0005 + +# Miss detection uses the DRDY bracket -- the span between the last poll showing +# DRDY clear and the poll showing it set -- rather than the read-to-read +# interval. +# +# The bracket is exact where the interval is only a heuristic. If DRDY reads +# clear at t_c then every earlier measurement has already been read, since a +# results read is what clears it. Measurements complete one period apart, so a +# bracket narrower than one period can contain at most one completion, and DRDY +# going high proves it contained at least one. Exactly one, whatever the host +# was doing beforehand. +# +# The read-to-read interval cannot make that claim: measured on this rig it +# reaches 37 ms against a 28.8 ms period -- a 29% overshoot from host stalls +# alone -- while the bracket stays under 10 ms. Thresholding the interval +# therefore flags healthy captures as lossy. +# +# The bracket is compared against the period measured by calibrate_period(), not +# the table value: the latter is 6-9% out on this unit, which is enough to +# miscount grid points inside a bracket and slip the sample index. + # Bounded so a stalled writer degrades predictably instead of exhausting memory # on a long run. Far above the depth a healthy writer ever reaches. QUEUE_MAX = 200_000 @@ -48,29 +78,111 @@ QUEUE_MAX = 200_000 _SENTINEL = object() +# Rows carry flags in a `warning` column, space separated, empty when fine. +# Keeping them out of the numeric columns means x/y/z stay parseable as numbers +# (blank for a row with no data) and the column generalises to future flags. +# +# MISSED this row is a placeholder -- the measurement was never read, so it +# has no data. The row exists to keep sample_index contiguous, so it +# remains a valid chip-time grid coordinate and the gap is explicit +# rather than silently compressing the timeline. +# AMBIGUOUS the gap ending at this row was of uncertain length, so the index +# may have slipped from here on. +# +# The two are independent. A gap measuring 1.35 periods rounds to one, so no +# placeholder is written -- yet it sits far enough from an integer to distrust, +# and that case carries AMBIGUOUS on the *real* sample terminating the gap, +# which keeps its data. A placeholder inside an uncertain gap carries both. +# +# The flag is a confidence measure, not a claim that 1.35 and 1.65 are equally +# likely: 1.35 probably is one period and 1.65 probably two. Deciding which +# needs neighbouring timestamps and assumptions, and can still be wrong, which +# is exactly why the rows are flagged rather than silently resolved. +WARN_MISSED = "MISSED" +WARN_AMBIGUOUS = "AMBIGUOUS" + +# How near a half-period the rounding may fall before the count is a coin toss. +AMBIGUITY_MARGIN = 0.25 + + def parse_args(): p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) - p.add_argument("--cycle-count", type=int, default=rm3100.DEFAULT_CYCLE_COUNT, - help="cycle count per axis; lower is faster but noisier. " - "50 allows ~540 Hz, 200 allows ~150 Hz (default: %(default)s)") - p.add_argument("--tmrc", type=lambda s: int(s, 0), default=0x96, + knob = p.add_mutually_exclusive_group() + knob.add_argument("--rate", type=float, default=None, + help="target sample rate in Hz. The cycle count and TMRC " + "are derived from it, which is usually what you want " + "-- the cycle count is a continuous rate knob where " + "TMRC offers only factor-of-two steps") + knob.add_argument("--cycle-count", type=int, default=None, + help="cycle count per axis, setting both the rate and the " + "resolution. Lower is faster but more coarsely " + f"quantised (default: {rm3100.DEFAULT_CYCLE_COUNT})") + p.add_argument("--tmrc", type=lambda s: int(s, 0), default=None, help="continuous-mode rate register, 0x92 (fastest) to 0x9D " - "(slowest) (default: 0x96, ~37 Hz)") + "(slowest). Default is the fastest, which lets the cycle " + "count set the rate -- give one only to sample slower " + "than the cycle count allows") p.add_argument("--duration", type=float, default=0.0, help="seconds to log, or 0 to run until Ctrl-C (default: %(default)s)") p.add_argument("--output", default=None, 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=100, - help="I2C bus speed in kHz; 400 measured best for high sample " - "rates (default: %(default)s)") + p.add_argument("--bus-speed", type=int, choices=[20, 100, 400, 750], default=750, + help="I2C bus speed in kHz. The default cycle count runs near " + "534 Hz, where 100 kHz would spend 80%% of each period on " + "the bus; 750 spends 11%% (default: %(default)s)") + p.add_argument("--calibrate", type=float, default=1.0, + help="seconds of loss-free samples used to measure the true " + "measurement period before recording starts; the run " + "aborts if no clean stretch can be found (default: " + "%(default)s)") + p.add_argument("--note", default=None, + help="free-text label recorded in the capture header, e.g. " + "the supply under test. Keeps the configuration with the " + "data instead of only in the filename") p.add_argument("--scan-only", action="store_true", help="scan the bus, report what responded, and exit") return p.parse_args() +def print_plan(cfg, bus_speed, requested_rate): + """Show how the configuration was derived, so it can be checked not trusted.""" + lsb = rm3100.tesla_per_count(cfg.cycle_count) * rm3100.NT_PER_TESLA + dither = rm3100.expected_noise_nt(cfg.cycle_count) / lsb + period = 1.0 / cfg.predicted_hz + bus = i2c_bus_time(bus_speed) + share = bus / period + + print("Configuration") + if requested_rate is not None: + print(f" requested {requested_rate:g} Hz") + print(f" cycle count {cfg.cycle_count:,}" + f"{'':<8}1 / (3 x ({cfg.cycle_count}/" + f"{rm3100.COUNTS_PER_SECOND:,.0f} + " + f"{rm3100.AXIS_OVERHEAD_S * 1e6:.1f} us))") + print(f" sensor ceiling {1 / rm3100.sample_period(cfg.cycle_count):.2f} Hz") + print(f" TMRC 0x{cfg.tmrc:02x} = " + f"{rm3100.TMRC_RATES[cfg.tmrc]:g} Hz -- {cfg.governed_by} governs") + print(f" predicted rate {cfg.predicted_hz:.2f} Hz " + f"+/-{rm3100.RATE_TOLERANCE:.0%} (oscillator tolerance; measured below)") + print(f" resolution {lsb:.2f} nT/LSB dither {dither:.2f} LSB " + "at spec noise") + print(f" duty {cfg.duty:.1%} integration / period") + print(f" bus {bus_speed} kHz {bus * 1e3:.3f} ms/sample " + f"{share:.1%} of the period") + for note in cfg.notes: + print(f" note: {note}") + # Bus speed is validated but never changed silently -- swapping it would + # hide a wiring problem behind a configuration change. + if share > 0.5: + faster = [s for s in sorted(BUS_SPEEDS) if i2c_bus_time(s) / period < 0.25] + fix = f"; {faster[0]} kHz would fit" if faster else "" + print(f" WARNING: the bus needs {share:.0%} of every period{fix}", + file=sys.stderr) + + def find_sensor(bus, address): """Locate the RM3100, or exit with wiring guidance.""" if address is not None: @@ -87,7 +199,7 @@ def find_sensor(bus, address): else: print("Nothing responded anywhere on the bus. Check, in order:\n" " - I2CEN tied high (otherwise the chip stays in SPI mode)\n" - " - the adapter's voltage jumper (the RM3100 wants ~3.3 V)\n" + " - AVDD/VDD powered, not just DVDD\n" " - SDA/SCL not swapped\n" " - SDA and SCL pull-up resistors present", file=sys.stderr) sys.exit(1) @@ -103,7 +215,7 @@ def find_sensor(bus, address): def start_sensor(sensor, cycle_count, tmrc): - """Identify and configure the sensor, printing each step.""" + """Identify and configure the sensor, printing each step. Returns REVID.""" revid = sensor.revid() if revid == rm3100.EXPECTED_REVID: print(f"REVID 0x{revid:02x} -- RM3100 confirmed") @@ -117,25 +229,29 @@ def start_sensor(sensor, cycle_count, tmrc): print(f"ERROR: cycle count read back as {readback}, expected " f"{(cycle_count,) * 3}", file=sys.stderr) sys.exit(1) - gain = rm3100.gain_lsb_per_ut(cycle_count) - print(f"Cycle counts set to {readback} -- gain {gain:.1f} LSB/uT") + lsb_nt = rm3100.tesla_per_count(cycle_count) * rm3100.NT_PER_TESLA + print(f"Cycle counts set to {readback} -- 1 LSB = {lsb_nt:.2f} nT") + # Raises if DRC1 did not take; exactly-once sampling depends on it. sensor.configure() sensor.set_rate(tmrc) - rate = rm3100.TMRC_RATES[tmrc] - print(f"Rate register 0x{tmrc:02x} -- requesting about {rate:g} Hz") + sensor.start_cmm() print("Continuous measurement mode started") - return rate + return revid -def writer_thread(q, path, gain, start_wall, stats): - """Drain raw samples: convert, format, write CSV, drive the console. +def write_header(handle, meta): + for key, value in meta.items(): + handle.write(f"# {key}: {value}\n") - Everything here is deliberately off the sampling thread. - """ + +def writer_thread(q, path, meta, stats): + """Drain samples to CSV and drive the console, off the sampling thread.""" + lsb_nt = float(meta["tesla_per_count"]) * rm3100.NT_PER_TESLA with open(path, "w", newline="") as handle: + write_header(handle, meta) out = csv.writer(handle) out.writerow(CSV_FIELDS) last_print = 0.0 @@ -145,72 +261,261 @@ def writer_thread(q, path, gain, start_wall, stats): item = q.get() if item is _SENTINEL: break - elapsed, cx, cy, cz = item + index, system_time, cx, cy, cz, warning = item - ux, uy, uz = cx / gain, cy / gain, cz / gain - magnitude = math.sqrt(ux * ux + uy * uy + uz * uz) - # Wall-clock is reconstructed from the monotonic offset rather than - # sampled per reading: one fewer syscall in the hot path, and immune - # to NTP steps mid-capture. - stamp = (start_wall + timedelta(seconds=elapsed)).isoformat() - - out.writerow([stamp, f"{elapsed:.5f}", cx, cy, cz, - f"{ux:.4f}", f"{uy:.4f}", f"{uz:.4f}", - f"{magnitude:.4f}"]) + # 6 decimal places is microsecond resolution on an epoch value, + # far finer than the ~1 ms bracket uncertainty on each timestamp. + # Placeholder rows carry zeros rather than blanks: the warning + # column already says the row has no data, so keeping x/y/z strictly + # integer makes the file trivial to parse. + out.writerow([index, f"{system_time:.6f}", + 0 if cx is None else cx, + 0 if cy is None else cy, + 0 if cz is None else cz, warning]) stats["rows"] += 1 now = time.monotonic() if now - last_flush >= FLUSH_INTERVAL_S: handle.flush() last_flush = now - if now - last_print >= CONSOLE_REFRESH_S: + if now - last_print >= CONSOLE_REFRESH_S and cx is not None: + # Placeholders carry no field values, so the display holds the + # last real reading rather than blanking or showing nonsense. last_print = now - print(f"\r{elapsed:8.2f}s X {ux:+9.3f} Y {uy:+9.3f} " - f"Z {uz:+9.3f} |B| {magnitude:8.3f} uT " - f"({stats['rows']} samples, " - f"{stats['rows'] / max(elapsed, 1e-9):.0f} Hz)", + x, y, z = cx * lsb_nt, cy * lsb_nt, cz * lsb_nt + print(f"\r{stats['rows']:,} samples ({stats['missed']} missed)" + f" | X {x:+10.1f} | Y {y:+10.1f} | Z {z:+10.1f}" + f" | T {(x * x + y * y + z * z) ** 0.5:10.1f} nT", end="", flush=True) handle.flush() -def sample_loop(sensor, q, duration, rate, stats): - """Read the sensor as fast as it produces data. I2C and nothing else.""" - sample_timeout = max(2.0, 5.0 / rate) +def i2c_bus_time(bus_speed_khz): + """Seconds of I2C traffic the bus carries per sample. + + A combined register read of n payload bytes clocks n+3 bytes -- write + address, register, read address, payload -- at 9 bit-times each (8 data plus + ACK), and spends roughly a bit-time on each START and STOP. One sample costs + a STATUS poll (1 payload byte) plus a results read (9). + + This is the irreducible traffic, not bus occupancy: the loop polls + continuously while waiting for DRDY, so raw occupancy approaches 100% and + says nothing useful. What matters is this against the sample period. + """ + bits = 0 + for payload in (1, 9): + bits += (payload + 3) * 9 + 3 # 2 STARTs and a STOP + return bits / (bus_speed_khz * 1000.0) + + +def calibrate_period(sensor, dt_nominal, seconds, attempts=3): + """Measure the true measurement period from a clean stretch of samples. + + Miss detection needs the real period, not the table value: the RC oscillator + runs up to 7% off (manual section 5.2.1) and this unit is ~6% slow, which is + enough to miscount grid points inside a bracket and slip the sample index. + + Least-squares slope of read time against index. Over one second at 35 Hz + (~35 points, ~0.5 ms jitter) that pins the period to roughly 0.03%, against + the 6-9% error of the nominal value -- and 0.03% is far finer than the + resolution needed to tell k grid points from k+1. + + Requires a stretch with no miss, since an interval spanning a miss is a + multiple of the period, not one period. Returns None if no clean stretch is + available, which is itself the answer: the rate is not sustainable. + + Also returns the mean host cost of one sample's traffic -- a DRDY poll plus + a results read -- timed here at no extra cost, since those are the very + transactions being performed. + """ + poll_ready, read_raw = sensor.poll_ready, sensor.read_raw + monotonic = time.monotonic + + # Bootstrap the period from observation rather than trusting the TMRC table. + # When the cycle count governs instead of TMRC -- which it does whenever the + # cycle-count ceiling falls below the requested rate (manual section 5.2.1) + # -- the table value is not merely 7% out but wrong by a large factor. At + # cc=400 with TMRC 0x92 the table says 1.67 ms and the truth is 13.5 ms, so + # a threshold built on it rejects every interval and calibration can never + # start. A handful of raw intervals settles it: the median is robust to the + # occasional stall, and only needs to be close enough to seed the real fit. + probe = [] + clear = monotonic() + deadline = clear + max(seconds, 0.25) + while monotonic() < deadline and len(probe) < 24: + ready_now, stamp, _ = poll_ready() + if not ready_now: + continue + read_raw() + probe.append(stamp) + if len(probe) >= 4: + gaps = sorted(b - a for a, b in zip(probe, probe[1:])) + dt_nominal = gaps[len(gaps) // 2] + + for _ in range(attempts): + times = [] + poll_cost = read_cost = 0.0 + polls = 0 + clear = monotonic() + deadline = clear + seconds + clean = True + while monotonic() < deadline: + before = monotonic() + ready_now, after, _ = poll_ready() + poll_cost += monotonic() - before + polls += 1 + if not ready_now: + clear = after + continue + ready = after + if ready - clear >= dt_nominal: # conservative: nominal is short + clean = False + break + read_raw() + read_cost += monotonic() - ready + times.append(ready) + clear = ready + n = len(times) + if not clean or n < 8: + continue + # slope of t against i for i = 0..n-1, where Sxx = n(n^2-1)/12. + mean_i = (n - 1) / 2.0 + mean_t = sum(times) / n + sxy = sum((i - mean_i) * (t - mean_t) for i, t in enumerate(times)) + return (sxy / (n * (n * n - 1) / 12.0), + poll_cost / max(polls, 1) + read_cost / n) + return None, None + + +def sample_loop(sensor, q, duration, dt_nominal, stats): + """Read every measurement, recording placeholders for any that are lost. + + The timestamp recorded per sample is the midpoint of the bracket between the + last poll showing DRDY clear and the poll showing it set: the measurement + completed somewhere in that window, and the midpoint is the best estimate + available without a hardware DRDY line. + """ read_raw = sensor.read_raw - data_ready = sensor.data_ready + poll_ready = sensor.poll_ready monotonic = time.monotonic put = q.put_nowait - start = monotonic() - timeouts = 0 - deadline = start + duration if duration > 0 else float("inf") - while monotonic() < deadline: - if not data_ready(): - if monotonic() - start > sample_timeout and stats["rows"] == 0: - timeouts += 1 - print(f"\nWARNING: no data ready within {sample_timeout:.1f} s", - file=sys.stderr) - if timeouts >= 3: - print("Giving up after 3 consecutive timeouts.", file=sys.stderr) - return + # 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. + if poll_ready()[0]: + read_raw() + + deadline = monotonic() + duration if duration > 0 else float("inf") + index = 0 + # Bracket lower bound: the most recent moment DRDY was seen clear. + clear_mono, clear_wall = monotonic(), time.time() + + # dt_nominal here is the *calibrated* period from calibrate_period(), not + # the table value. Keep refining it from clean intervals so the estimate + # follows the oscillator's thermal drift (~2500 ppm measured over 13 h). + dt_est = dt_nominal + previous_mono = previous_wall = None + DT_SMOOTHING = 0.02 + + while True: + if monotonic() >= deadline: + return + # The stamps come back with the STATUS reply itself, taken the instant + # it landed rather than several frames later up here. + ready, ready_mono, ready_wall = poll_ready() + if not ready: + clear_mono, clear_wall = ready_mono, ready_wall continue + # Monotonic drives detection so an NTP step cannot fake a miss; the + # wall clock is recorded so a real NTP correction stays visible. + bracket = ready_mono - clear_mono + if bracket > stats["max_bracket"]: + stats["max_bracket"] = bracket - cx, cy, cz = read_raw() + # Counting lost measurements in two stages, because the two questions + # have different best answers. + # + # "Did we lose any?" is settled by the bracket, and rigorously: a + # bracket shorter than one period cannot contain two completions, so + # nothing was lost regardless of timing precision. + # + # "How many?" is better answered by the interval since the previous + # accepted sample. Each accepted read sits on a grid point, so that + # interval is very close to an exact multiple of the period, and + # rounding recovers the multiple. The bracket alone cannot -- it only + # spans back to the last poll that saw DRDY clear, which discards where + # the grid actually is. + # + # Precision: the period is calibrated to ~0.07% and each completion is + # located to about half a poll interval (~0.2 ms, under 1% of a period). + # Over a gap of k periods the total error is well under a tenth of a + # period for any small k, so the rounding is unambiguous -- which is + # the case that matters, since a long gap means the run is unusable + # anyway. + lost = 0 + uncertain = False + if bracket >= dt_est and previous_mono is not None: + periods = (ready_mono - previous_mono) / dt_est + lost = max(0, round(periods) - 1) + # Uncertainty is a property of the gap, not of the placeholders: + # a gap rounding to zero losses can still be a coin toss, and then + # there is no placeholder to carry the flag. + uncertain = abs(periods - round(periods)) > AMBIGUITY_MARGIN + if uncertain: + stats["ambiguous"] += 1 + elif previous_mono is not None: + # A clean interval is exactly one period, so it calibrates dt_est. + # Only clean ones qualify: an interval spanning a miss is a multiple. + dt_est += (ready_mono - previous_mono - dt_est) * DT_SMOOTHING + if lost: + stats["missed"] += lost + # Placeholders sit on the grid -- one period after the previous + # accepted sample, and so on -- because that is where the lost + # measurements actually completed. Spreading them across the + # bracket instead would bunch them at the end of the gap and skew + # the rate fit. + flags = (f"{WARN_MISSED} {WARN_AMBIGUOUS}" if uncertain + else WARN_MISSED) + for k in range(lost): + try: + put((index, previous_wall + dt_est * (k + 1), + None, None, None, flags)) + except queue.Full: + stats["dropped"] += 1 + index += 1 + + sample_wall = (clear_wall + ready_wall) / 2.0 + (cx, cy, cz), _, _ = read_raw() try: - put((monotonic() - start, cx, cy, cz)) + # An uncertain gap that produced no placeholder still has to be + # flagged, so it rides on the real sample that ends it. + put((index, sample_wall, cx, cy, cz, + WARN_AMBIGUOUS if uncertain and not lost else "")) except queue.Full: stats["dropped"] += 1 + index += 1 + previous_mono, previous_wall = ready_mono, sample_wall + # The read just cleared DRDY, so this is a known-clear instant. + clear_mono, clear_wall = ready_mono, ready_wall def main(): args = parse_args() + sys.setswitchinterval(GIL_SWITCH_INTERVAL_S) - if args.cycle_count < 1 or args.cycle_count > 0xFFFF: + if args.cycle_count is not None and not 1 <= args.cycle_count <= 0xFFFF: sys.exit(f"--cycle-count {args.cycle_count} outside 1..65535") - if args.tmrc not in rm3100.TMRC_RATES: + if args.rate is not None and args.rate <= 0: + sys.exit(f"--rate {args.rate:g} must be positive") + if args.tmrc is not None and args.tmrc not in rm3100.TMRC_RATES: sys.exit(f"--tmrc 0x{args.tmrc:02x} is not a valid rate register value") + cfg = rm3100.plan(rate_hz=args.rate, cycle_count=args.cycle_count, + tmrc=args.tmrc) + args.cycle_count, args.tmrc = cfg.cycle_count, cfg.tmrc + try: bus = ch347.CH347I2C(BUS_SPEEDS[args.bus_speed]) except IOError as exc: @@ -223,26 +528,84 @@ def main(): found = bus.scan() print("Devices found: " + ", ".join(f"0x{a:02x}" for a in found) if found else "No devices responded on the bus.") - return + return 0 address = find_sensor(bus, args.address) sensor = rm3100.RM3100(bus, address) - rate = start_sensor(sensor, args.cycle_count, args.tmrc) + try: + revid = start_sensor(sensor, args.cycle_count, args.tmrc) + except IOError as exc: + sys.exit(str(exc)) + + print() + print_plan(cfg, args.bus_speed, args.rate) + print() + nominal_rate = cfg.predicted_hz + + # Establish the real period before recording. Everything downstream -- + # miss detection, placeholder counting, the sample-index grid -- depends + # on it, and the table value is 6-9% out on this unit. + dt, host_cost = calibrate_period(sensor, 1.0 / nominal_rate, + args.calibrate) + if dt is None: + sensor.stop_cmm() + sys.exit( + f"Could not find {args.calibrate:g} s of loss-free samples at " + f"{nominal_rate:.1f} Hz.\n" + "Without a clean stretch the true period cannot be measured, so " + "misses cannot be\ncounted reliably and the sample index would " + "not track chip time.\n" + "Lower the rate (higher --tmrc) or raise --bus-speed.") + # Compared against the predicted rate, not TMRC's nominal: when the + # cycle count governs, the TMRC figure is not what was aimed for and the + # error against it is meaningless. + print(f"Calibrated period {dt * 1e3:.4f} ms = {1 / dt:.4f} Hz " + f"({(1 / dt) / cfg.predicted_hz - 1:+.2%} vs predicted)") + # Whether the host is the constraint is a measurement, not a given: it + # costs a near-constant ~0.9 ms per sample regardless of configuration, + # so it dominates at short periods and disappears at long ones. + share = host_cost / dt + margin = ("little margin -- expect misses under load" if share > 0.5 + else "modest margin" if share > 0.25 + else "ample margin") + print(f"Host cost {host_cost * 1e3:.3f} ms/sample " + f"({share:.1%} of the period, {margin})") + + meta = { + "rm3100_capture": 1, + # The rate this configuration is predicted to produce, whichever of + # the two ceilings governs -- not TMRC's table value, which is not + # what was aimed for when the cycle count wins. + "nominal_rate_hz": nominal_rate, + "tmrc_nominal_hz": rm3100.TMRC_RATES[args.tmrc], + "tmrc": f"0x{args.tmrc:02x}", + "cycle_count": args.cycle_count, + # repr() so the constant round-trips through float64 exactly. + "tesla_per_count": repr(rm3100.tesla_per_count(args.cycle_count)), + "i2c_address": f"0x{address:02x}", + "bus_speed_khz": args.bus_speed, + "revid": f"0x{revid:02x}", + # What the logger actually used for miss detection. Not recoverable + # from the data, since capture.py fits the whole run rather than the + # first second. + "calibrated_period_s": repr(dt), + } + if args.note: + # Newlines would break the one-line-per-key header format. + meta["note"] = " ".join(args.note.split()) path = args.output or datetime.now().strftime("rm3100_%Y%m%d_%H%M%S.csv") print(f"Logging to {path} -- Ctrl-C to stop\n") q = queue.Queue(maxsize=QUEUE_MAX) - stats = {"rows": 0, "dropped": 0} - writer = threading.Thread( - target=writer_thread, - args=(q, path, rm3100.gain_lsb_per_ut(args.cycle_count), - datetime.now(timezone.utc), stats), - daemon=True) + stats = {"rows": 0, "dropped": 0, "missed": 0, "ambiguous": 0, + "max_bracket": 0.0} + writer = threading.Thread(target=writer_thread, + args=(q, path, meta, stats), daemon=True) writer.start() try: - sample_loop(sensor, q, args.duration, rate, stats) + sample_loop(sensor, q, args.duration, dt, stats) except KeyboardInterrupt: pass finally: @@ -254,11 +617,38 @@ def main(): print("Continuous measurement mode stopped") 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}") + + 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. + # Quoted against the calibrated period, which is the threshold actually + # enforced -- the nominal one is 6-9% out. + print(f"Worst DRDY bracket {stats['max_bracket'] * 1e3:.2f} ms of " + f"{dt * 1e3:.2f} ms allowed " + f"({stats['max_bracket'] / dt * 100:.0f}% of margin used)") + + if stats["ambiguous"]: + print(f"\nWARNING: {stats['ambiguous']:,} gap(s) could not be " + "counted confidently -- the interval fell near a half-period, " + "so\nthe number of lost measurements is a guess and " + f"sample_index may have slipped.\nThose rows carry " + f"{WARN_AMBIGUOUS} in the warning column. Do not trust this " + "capture's\ntime axis for spectral work.", file=sys.stderr) + + if stats["missed"]: + pct = stats["missed"] / max(stats["rows"], 1) * 100 + print(f"\nWARNING: {stats['missed']:,} measurement(s) were lost " + f"({pct:.3f}% of rows), flagged {WARN_MISSED}.\n" + "The sample index still tracks chip time -- the gaps are " + "explicit, not compressed --\nbut those rows carry no field " + "data. Lower the rate (higher --tmrc) or raise --bus-speed\n" + "to remove them.", file=sys.stderr) + return 1 + return 0 if __name__ == "__main__": - main() + sys.exit(main()) diff --git a/noise_floor_test.sh b/noise_floor_test.sh new file mode 100755 index 0000000..af3e215 --- /dev/null +++ b/noise_floor_test.sh @@ -0,0 +1,94 @@ +#!/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/plot.py b/plot.py index b15edb1..4e28068 100644 --- a/plot.py +++ b/plot.py @@ -11,7 +11,6 @@ independently, and note the per-panel mean/sd annotation for context. """ import argparse -import csv import sys import matplotlib @@ -19,6 +18,8 @@ matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np +import capture + # Light-mode design tokens. SURFACE = "#fcfcfb" TEXT_PRIMARY = "#0b0b0b" @@ -31,33 +32,13 @@ GRID = "#e3e2df" # also keeps the categorical set at the three slots that validate for small # multiples. SERIES = [ - ("x_uT", "X axis", "#2a78d6"), - ("y_uT", "Y axis", "#eb6834"), - ("z_uT", "Z axis", "#1baf7a"), - (None, "Norm |B| = sqrt(X^2 + Y^2 + Z^2)", TEXT_PRIMARY), + ("x", "X axis", "#2a78d6"), + ("y", "Y axis", "#eb6834"), + ("z", "Z axis", "#1baf7a"), + ("total", "Norm |B| = sqrt(X^2 + Y^2 + Z^2)", TEXT_PRIMARY), ] -def load(path): - t, x, y, z, gains = [], [], [], [], [] - with open(path, newline="") as fh: - for row in csv.DictReader(fh): - t.append(float(row["elapsed_s"])) - x.append(float(row["x_uT"])) - y.append(float(row["y_uT"])) - z.append(float(row["z_uT"])) - # Recover the gain from the raw/uT ratio so the caption reports the - # settings actually used rather than an assumption. - for axis in "xyz": - ut = float(row[f"{axis}_uT"]) - if abs(ut) > 1.0: - gains.append(float(row[f"{axis}_raw"]) / ut) - if not t: - sys.exit(f"{path} contains no samples") - gain = float(np.median(gains)) if gains else float("nan") - return np.array(t), np.array(x), np.array(y), np.array(z), gain - - def rolling_mean(v, window): """Centred moving average that stays smooth all the way to both ends. @@ -83,14 +64,16 @@ def main(): "(default: %(default)s)") args = ap.parse_args() - t, x, y, z, gain = load(args.csv) - # Inverse of rm3100.gain_lsb_per_ut(). - cycle_count = (gain - 1.5) / 0.3671 - norm = np.sqrt(x**2 + y**2 + z**2) - series = {"x_uT": x, "y_uT": y, "z_uT": z, None: norm} + try: + cap = capture.load(args.csv) + except (OSError, capture.CaptureError) as exc: + sys.exit(str(exc)) - duration = t[-1] - t[0] - rate = len(t) / duration if duration > 0 else float("nan") + # Drift-corrected time: the nominal grid would be ~6% off real seconds. + t = cap.elapsed + series = dict(cap.axes()) + duration = cap.duration + rate = cap.true_rate_hz window = max(1, int(round(args.smooth * rate))) if args.smooth > 0 else 0 fig, axes = plt.subplots(4, 1, figsize=(12, 9.5), sharex=True, dpi=150) @@ -108,13 +91,13 @@ def main(): ax.plot(t, rolling_mean(v, window), color=color, linewidth=1.6, solid_capstyle="round") - ax.set_ylabel("µT", color=TEXT_SECONDARY, fontsize=10) + ax.set_ylabel("nT", color=TEXT_SECONDARY, fontsize=10) # Direct label instead of a legend: one series per panel, so the title # names it. This is also the relief the palette's contrast WARN requires. ax.set_title(label, color=TEXT_PRIMARY, fontsize=12, loc="left", pad=8, fontweight="bold") - ax.annotate(f"mean {v.mean():.3f} sd {v.std() * 1000:.0f} nT " - f"span {v.max() - v.min():.3f} µT", + ax.annotate(f"mean {v.mean():,.0f} nT sd {v.std():.0f} nT " + f"span {v.max() - v.min():,.0f} nT", xy=(1.0, 1.0), xycoords="axes fraction", xytext=(0, 8), textcoords="offset points", ha="right", va="bottom", @@ -137,8 +120,9 @@ def main(): fig.suptitle("RM3100 magnetometer capture", color=TEXT_PRIMARY, fontsize=15, fontweight="bold", x=0.5, y=0.985) fig.text(0.5, 0.955, - f"{len(t):,} samples over {duration:.1f} s ({rate:.0f} Hz), " - f"cycle count {cycle_count:.0f} ({gain:.1f} LSB/µT)" + f"{len(t):,} samples over {duration:.1f} s at {rate:.2f} Hz " + f"(nominal {cap.nominal_rate_hz:g} Hz, {cap.rate_error * 100:+.1f}%), " + f"cycle count {cap.cycle_count} (1 LSB = {cap.lsb_nt:.1f} nT)" f"{smooth_note}. Panels have independent y-scales.", color=TEXT_SECONDARY, fontsize=10, ha="center") @@ -148,8 +132,8 @@ def main(): print(f"{len(t):,} samples, {duration:.2f} s, {rate:.1f} Hz -> {out}") for key, label, _ in SERIES: v = series[key] - print(f" {label.split()[0]:5s} mean {v.mean():+9.3f} uT " - f"sd {v.std()*1000:6.1f} nT span {v.max()-v.min():6.3f} uT") + print(f" {label.split()[0]:5s} mean {v.mean():+12,.1f} nT " + f"sd {v.std():8.1f} nT span {v.max()-v.min():10,.1f} nT") if __name__ == "__main__": diff --git a/rm3100.py b/rm3100.py index 981d863..c16cbf5 100644 --- a/rm3100.py +++ b/rm3100.py @@ -8,6 +8,7 @@ and read(addr, count). """ import time +from collections import namedtuple # Register addresses (manual Table 5-1). REG_POLL = 0x00 # single measurement trigger @@ -45,16 +46,183 @@ STATUS_DRDY = 0x80 EXPECTED_REVID = 0x22 -DEFAULT_CYCLE_COUNT = 200 +# 50 buys bandwidth almost for free. Against cycle count 200 it costs 4.4% in +# post-filter noise (89% duty against 97%) and nothing else, while giving 3.7x +# the spectrum -- 267 Hz of Nyquist against 73 Hz. That matters because aliased +# interference cannot be filtered out afterwards at any cycle count, so seeing +# it is worth more than a few percent of noise. +# +# The one thing to verify rather than assume is dither: 50 sits at 0.58 LSB of +# intrinsic noise, which simulation puts safely in the region where averaging +# still recovers sub-LSB resolution (it fails below ~0.2). characterize.py +# prints sd/LSB, which answers it from the first capture. Fall back to 100 or +# 200 if that comes back low. +DEFAULT_CYCLE_COUNT = 50 + +# The fastest rate register. Used by default so the cycle count, not TMRC, sets +# the rate -- TMRC offers only factor-of-two steps and has no effect at all once +# the cycle count governs. +TMRC_FASTEST = 0x92 -def gain_lsb_per_ut(cycle_count): - """Sensitivity in LSB per microtesla for a given cycle count. +# Table 3-1 quotes gain in LSB per microtesla; these are the coefficients of a +# linear fit to its 20, 38 and 75 LSB/uT at cycle counts 50, 100 and 200, which +# reproduces all three to within a count. Working units in this codebase are +# tesla, so the fit is scaled by 1e6 on the way out. +_GAIN_SLOPE_LSB_PER_UT = 0.3671 +_GAIN_OFFSET_LSB_PER_UT = 1.5 +UT_PER_TESLA = 1e6 - Linear fit to manual Table 3-1, which quotes 20, 38 and 75 LSB/uT at cycle - counts of 50, 100 and 200; this reproduces all three to within a count. +# Displays and plots use nanotesla; no computation is done in it. +NT_PER_TESLA = 1e9 + + +def gain_lsb_per_tesla(cycle_count): + """Sensitivity in LSB per tesla for a given cycle count.""" + return (_GAIN_SLOPE_LSB_PER_UT * cycle_count + + _GAIN_OFFSET_LSB_PER_UT) * UT_PER_TESLA + + +def tesla_per_count(cycle_count): + """Calibration constant: multiply a raw count by this to get tesla. + + Stored in capture headers, so a reader needs no knowledge of the gain fit -- + and it is a multiply at the point of use rather than a divide. """ - return 0.3671 * cycle_count + 1.5 + return 1.0 / gain_lsb_per_tesla(cycle_count) + + +# Timing model for a three-axis measurement: +# +# per-axis time = cycle_count / COUNTS_PER_SECOND + AXIS_OVERHEAD_S +# +# COUNTS_PER_SECOND comes from the specification, not from a fit: Table 3-1 +# gives a 180 kHz circuit oscillation frequency, and section 4.1 measures each +# cycle count in *both* bias directions, so one count costs two oscillations. +# Rates measured here across cycle counts 229..29769 agree with that figure to +# within 0.8%, which confirms the spec value rather than improving on it -- the +# residual is this particular unit's oscillator sitting inside ordinary +# component tolerance, and another part would sit somewhere else. +# +# The overhead has no specified value and must be measured; 68.7 us reproduces +# the observed rates to better than 1% when the divisor is held at spec. It is +# why the naive "rate x cycle_count" constant is not constant, drifting from +# 84,429 at cycle count 100 to 89,191 at 1200. +COUNTS_PER_SECOND = 90000.0 # 180 kHz (Table 3-1) / 2 bias directions +AXIS_OVERHEAD_S = 68.7e-6 +AXES = 3 + +# Rates predicted from the model are good to roughly this much on a given unit, +# and no better across units, since the manual quotes no tolerance on the +# oscillator. Anything needing the real number measures it: logger.py calibrates +# the period against the host clock before recording. +RATE_TOLERANCE = 0.02 + +# Table 3-1 quotes 30/20/15 nT at cycle counts 50/100/200, which fits +# K/sqrt(cycle_count). Extrapolation past ~400 is unverified: the manual calls +# that its useful upper limit for noise and gives no data beyond it. +_NOISE_K_NT = 208.0 + + +def sample_period(cycle_count, axes=AXES): + """Seconds between measurements when the cycle count governs the rate.""" + return axes * (cycle_count / COUNTS_PER_SECOND + AXIS_OVERHEAD_S) + + +def integration_time(cycle_count, axes=AXES): + """Seconds per sample actually spent integrating, excluding fixed overhead. + + Only this part reduces noise. Against the sample period it gives the duty + cycle: how much of the wall clock the sensor is doing useful work rather + than idling between measurements or paying per-axis overhead. + """ + return axes * cycle_count / COUNTS_PER_SECOND + + +def cycle_count_for_rate(rate_hz, axes=AXES): + """Cycle count that makes the sensor free-run at about rate_hz. + + Inverse of sample_period(). Use with a TMRC faster than the target so the + cycle count governs -- then this is a continuous rate knob, where TMRC only + offers factor-of-two steps, and the duty cycle is ~100% by construction. + + The achieved rate will sit within about RATE_TOLERANCE of the target, set by + oscillator tolerance rather than by this calculation. Measure it if it + matters. + """ + cc = round((1.0 / (rate_hz * axes) - AXIS_OVERHEAD_S) * COUNTS_PER_SECOND) + return max(1, min(0xFFFF, cc)) + + +# Section 5.1: "quantization issues generally dictate working above a cycle +# count value of ~30". Below that the LSB grows faster than the sensor's own +# noise, so the quantiser stops being dithered and averaging stalls. +MIN_CYCLE_COUNT = 30 +MAX_CYCLE_COUNT = 0xFFFF +MIN_RATE_BY_CYCLE_COUNT = 1.0 / (AXES * (MAX_CYCLE_COUNT / COUNTS_PER_SECOND + + AXIS_OVERHEAD_S)) + +Plan = namedtuple("Plan", "cycle_count tmrc predicted_hz governed_by duty notes") + + +def plan(rate_hz=None, cycle_count=None, tmrc=None): + """Resolve a full configuration from whichever knob the caller specified. + + Two ceilings compete and **the slower one wins** (manual section 5.2.1): the + cycle count sets how long a measurement takes, TMRC sets how often one is + started. Leaving TMRC faster than the cycle-count ceiling means the sensor + free-runs at ~100% duty; setting it slower makes the sensor idle, which + costs sensitivity for nothing unless low power is the goal. + + So the default in every branch is TMRC_FASTEST, and the cycle count is the + rate knob -- it is continuous where TMRC offers only factor-of-two steps. + """ + notes = [] + if cycle_count is None: + if rate_hz is None: + cycle_count = DEFAULT_CYCLE_COUNT + elif rate_hz < MIN_RATE_BY_CYCLE_COUNT: + # The register is 16 bits, so below ~0.46 Hz the cycle count runs + # out of range and TMRC is the only way to go slower. Max the cycle + # count anyway: it costs nothing and buys resolution. + cycle_count = MAX_CYCLE_COUNT + notes.append( + f"{rate_hz:g} Hz is below the {MIN_RATE_BY_CYCLE_COUNT:.3f} Hz " + f"floor of a {MAX_CYCLE_COUNT:,}-count measurement, so TMRC " + "must set the cadence and the sensor will idle") + if tmrc is None: + tmrc = min(TMRC_RATES, key=lambda t: abs(TMRC_RATES[t] - rate_hz)) + else: + cycle_count = cycle_count_for_rate(rate_hz) + + if cycle_count < MIN_CYCLE_COUNT: + notes.append(f"cycle count raised to the {MIN_CYCLE_COUNT} the manual " + "advises as a quantisation floor (section 5.1)") + cycle_count = max(MIN_CYCLE_COUNT, min(MAX_CYCLE_COUNT, cycle_count)) + ceiling = 1.0 / sample_period(cycle_count) + if tmrc is None: + tmrc = TMRC_FASTEST + requested = TMRC_RATES[tmrc] + + if requested <= ceiling: + predicted, governed_by = requested, "TMRC" + idle = 1.0 - requested / ceiling + if idle > 0.05: + notes.append(f"TMRC leaves the sensor idle {idle:.0%} of each " + "period, which costs sensitivity") + else: + predicted, governed_by = ceiling, "cycle count" + + if rate_hz is not None and predicted < rate_hz * 0.98: + notes.append(f"{rate_hz:g} Hz is faster than this configuration can " + f"reach; {predicted:.1f} Hz is the ceiling") + return Plan(cycle_count, tmrc, predicted, governed_by, + integration_time(cycle_count) * predicted, notes) + + +def expected_noise_nt(cycle_count): + """Per-sample noise from the Table 3-1 fit, in nanotesla.""" + return _NOISE_K_NT / cycle_count ** 0.5 def decode_measurements(data): @@ -84,10 +252,17 @@ class RM3100: def read_reg(self, reg, count=1): """Read count bytes starting at reg, using the sensor's auto-increment. - The pointer write is a separate transaction terminated by STOP rather - than a repeated START, which is exactly what the manual's I2C read - diagrams (sections 4.5.2 and 5.8.4) specify. + Returns (data, mono, wall): the bus stamps the clocks when the reply + lands, and the stamp travels with its own data. + + Prefers a single combined transaction (repeated START) where the bus + offers one, since the host round trip is what limits the sample rate. + Falls back to the manual's STOP-then-START form (sections 4.5.2, 5.8.4) + for a bus that cannot do it. """ + combined = getattr(self.bus, "write_read", None) + if combined is not None: + return combined(self.address, [reg], count) self.bus.write(self.address, [reg]) return self.bus.read(self.address, count) @@ -95,7 +270,7 @@ class RM3100: self.bus.write(self.address, bytes([reg]) + bytes(data)) def revid(self): - return self.read_reg(REG_REVID)[0] + return self.read_reg(REG_REVID)[0][0] def set_cycle_counts(self, count): """Set all three axes to the same cycle count.""" @@ -106,7 +281,7 @@ class RM3100: def get_cycle_counts(self): """Read back (ccx, ccy, ccz).""" - data = self.read_reg(REG_CCX, 6) + data, _, _ = self.read_reg(REG_CCX, 6) return tuple( int.from_bytes(data[i:i + 2], "big") for i in (0, 2, 4) ) @@ -120,8 +295,21 @@ class RM3100: self.write_reg(REG_TMRC, [tmrc]) def configure(self): - """Put DRDY into a state where polling STATUS actually works.""" + """Put DRDY into a state where polling STATUS actually works. + + Verified rather than assumed: exactly-once sampling relies on DRC1=1 + clearing DRDY when the measurement registers are read, which is what + makes a second read of the same measurement impossible. + """ self.write_reg(REG_HSHAKE, [HSHAKE_DRDY_ON_READ_ONLY]) + # Bits 4-6 are read-only NACK status, so compare only the writable ones. + readback = self.read_reg(REG_HSHAKE)[0][0] & 0x0F + if readback != HSHAKE_DRDY_ON_READ_ONLY & 0x0F: + raise IOError( + f"HSHAKE did not take: wrote 0x{HSHAKE_DRDY_ON_READ_ONLY:02x}, " + f"read back 0x{readback:02x}. Exactly-once sampling cannot be " + "guaranteed without DRC1=1." + ) def start_cmm(self): self.write_reg(REG_CMM, [CMM_ALL_AXES]) @@ -130,7 +318,17 @@ class RM3100: self.write_reg(REG_CMM, [CMM_OFF]) def data_ready(self): - return bool(self.read_reg(REG_STATUS)[0] & STATUS_DRDY) + """True if a measurement is waiting.""" + return self.poll_ready()[0] + + def poll_ready(self): + """Like data_ready(), but returns (ready, mono, wall). + + The stamp is when the STATUS reply landed, which is the tightest bound + available on when DRDY actually went high. + """ + data, mono, wall = self.read_reg(REG_STATUS) + return bool(data[0] & STATUS_DRDY), mono, wall def wait_for_data(self, timeout=2.0, interval=0.001): """Block until DRDY is set. Returns False if timeout elapses first.""" @@ -143,15 +341,16 @@ class RM3100: time.sleep(interval) def read_raw(self): - """Return (x, y, z) as signed counts -- the fast path. + """Return ((x, y, z) counts, mono, wall) -- the fast path. Callers logging at high rates should use this and defer the microtesla conversion, so the sampling loop does I2C and nothing else. """ - return decode_measurements(self.read_reg(REG_MX, 9)) + data, mono, wall = self.read_reg(REG_MX, 9) + return decode_measurements(data), mono, wall def read_measurements(self): - """Return ((x, y, z) counts, (x, y, z) microtesla).""" - counts = self.read_raw() - gain = gain_lsb_per_ut(self.cycle_count) + """Return ((x, y, z) counts, (x, y, z) tesla).""" + counts, _, _ = self.read_raw() + gain = gain_lsb_per_tesla(self.cycle_count) return counts, tuple(c / gain for c in counts) diff --git a/sweep.py b/sweep.py new file mode 100644 index 0000000..f175374 --- /dev/null +++ b/sweep.py @@ -0,0 +1,144 @@ +#!/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 --from 8 --to 128 + ./.venv/bin/python sweep.py --rates 10,25,50 + +Cycle count is the rate knob, not TMRC. TMRC only offers factor-of-two steps +and, once the cycle count governs, has no effect at all -- measured 73.85 Hz at +TMRC 0x92 against 73.86 Hz at 0x94 for the same cycle count. So each point runs +with TMRC set fast and the cycle count chosen from rm3100.cycle_count_for_rate(), +which makes the rate continuous and the duty cycle ~100% by construction. + +Each row is measured, not predicted: the logger calibrates the true period +against the host clock before recording, and that figure is what is reported. +""" + +import argparse +import re +import subprocess +import sys +from pathlib import Path + +import logger +import rm3100 + +HERE = Path(__file__).resolve().parent +CALIBRATED = re.compile(r"Calibrated period [\d.]+ ms = ([\d.]+) Hz") +MISSED = re.compile(r"WARNING: ([\d,]+) measurement\(s\) were lost") +AMBIGUOUS = re.compile(r"WARNING: ([\d,]+) gap\(s\) could not be counted") + +CALIBRATE_SECONDS = 5.0 +RECORD_SECONDS = 5.0 + +# One width per column, shared by the header rule and every cell, so the table +# cannot drift out of alignment as the cell formats change. +COLUMNS = [("target", 9), ("measured", 11), ("err", 7), ("cycle count", 11), + ("nT/LSB", 10), ("spec noise", 10), ("duty", 6), ("bus use", 8)] + + +def row(cells): + return " ".join(f"{c:>{width}}" for c, (_, width) in zip(cells, COLUMNS)) + + +def measure(rate, bus_speed, tmrc, output, record_seconds): + """Run one configuration. Returns (cycle_count, measured_hz or None, note).""" + cc = rm3100.cycle_count_for_rate(rate) + proc = subprocess.run( + [sys.executable, "logger.py", + "--cycle-count", str(cc), "--tmrc", hex(tmrc), + "--bus-speed", str(bus_speed), + "--calibrate", f"{CALIBRATE_SECONDS:g}", + "--duration", f"{record_seconds:g}", + "--output", str(output)], + capture_output=True, text=True, cwd=HERE) + out = proc.stdout + proc.stderr + found = CALIBRATED.search(out) + if found: + # Both warnings matter and they are independent. Reporting only losses + # would hide the worse case: an ambiguous gap means the sample index + # itself may have slipped, so the timeline is suspect even where no + # measurement was lost. + notes = [] + for pattern, label in ((MISSED, "lost"), (AMBIGUOUS, "ambiguous")): + hit = pattern.search(out) + if hit: + notes.append(f"{hit.group(1)} {label}") + return cc, float(found.group(1)), ", ".join(notes) + if "loss-free" in out: + return cc, None, "host cannot sustain" + return cc, None, "failed" + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--from", dest="low", type=float, default=2.0, + help="lowest target rate in Hz (default: %(default)s)") + ap.add_argument("--to", dest="high", type=float, default=512.0, + 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("--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, + help="seconds to record at each point, after calibration " + "(default: %(default)s)") + ap.add_argument("--output", default="/tmp/sweep_point.csv", + help="scratch capture path, overwritten each point") + args = ap.parse_args() + + if args.rates: + rates = [float(r) for r in args.rates.split(",")] + else: + rates, r = [], args.low + while r <= args.high * 1.001: + rates.append(r) + r *= 2 + + bus_time = logger.i2c_bus_time(args.bus_speed) + print(f"TMRC {hex(args.tmrc)} (held fast), I2C {args.bus_speed} kHz, " + f"{bus_time * 1e3:.3f} ms of bus traffic per sample, " + f"{args.duration:g} s per config\n") + print(row(name for name, _ in COLUMNS)) + print(row("-" * width for _, width in COLUMNS)) + + for rate in rates: + cc, hz, note = measure(rate, args.bus_speed, args.tmrc, args.output, + args.duration) + # One count is one LSB, so nT/LSB is the quantisation step directly. + nt_per_lsb = rm3100.tesla_per_count(cc) * rm3100.NT_PER_TESLA + noise = rm3100.expected_noise_nt(cc) + # Fraction of the period spent integrating: what actually reduces noise. + duty = rm3100.integration_time(cc) * hz if hz else None + # A note always trails the row rather than sitting in a cell: "host + # cannot sustain" is wider than any sensible column and would shove the + # rest of the line out of alignment. + print(row([ + f"{rate:g} Hz", + f"{hz:.2f} Hz" if hz else "-", + f"{hz / rate - 1:+.1%}" if hz else "-", + f"{cc:,}", + f"{nt_per_lsb:.3f} nT", + f"{noise:.2f} nT", + f"{duty:.1%}" if duty else "-", + f"{bus_time * hz:.1%}" if hz else "-", + ]) + (f" {note}" if note else "")) + + print("\nspec noise is Table 3-1's figure for the cycle count, not a\n" + "measurement: 208/sqrt(cc) nT, fitted to its 30/20/15 nT at 50/100/200.\n" + "Past cycle count ~400 the manual gives no data, so those are\n" + "extrapolation. Measure the real floor with characterize.py.\n" + "duty is integration time against the period -- what reduces noise.\n" + "It falls at high rates as the fixed per-axis overhead grows relative\n" + "to the integration, and would fall further if TMRC rather than the\n" + "cycle count governed, leaving the sensor idle between measurements.\n" + "bus use is irreducible I2C traffic against the period -- not occupancy,\n" + "which approaches 100% because the loop polls continuously for DRDY.") + + +if __name__ == "__main__": + main()