diff --git a/CHANGELOG.md b/CHANGELOG.md index 3a742b9..2650bdd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,66 @@ All notable changes to seismo-relay are documented here. --- +## Unreleased + +**Blastware Event/FFT-Report parity — the FFT, the USBM compliance chart, and +the sensor self-check.** Three analyses Blastware derives from event data, +reverse-engineered against BE12844 (MiniMate Plus) reports and reproduced in +seismo-relay: the compliance chart and the sensor-check strip now render on +the event-report PDF, and the FFT reproduces Blastware's FFT Report. All three +are additive and read from data already on disk — the `.h5` samples and the +retained raw BW binary — so there is **no `.h5`/DB change, no migration, and no +backfill**: a report regenerated for an existing event simply gains the new +panels. + +### Added +- **Blastware-compatible channel FFT (`waveform_fft`).** Reproduces Blastware's + FFT Report: DC-removed, no window, zero-padded to 4096 (0.25 Hz bins at + 1024 sps), single-sided `2/N` amplitude. Matches Blastware's dominant + frequency to the exact bin and the amplitude to report precision across all + 28 channels of the 7-event BE12844 oracle set. `channel_spectrum()` / + `dominant_frequency()`; tests in `tests/test_waveform_fft.py`. + +- **USBM RI8507 / OSMRE compliance chart on the event-report PDF + (`sfm/compliance.py`).** The velocity-vs-frequency blasting-compliance + scatter Blastware draws in the upper-right of its Event Report: each channel's + significant cycles as `(frequency, peak velocity)` points (zero-crossing + method, so each channel's cloud tops out at its PPV) plotted against the + RI8507 Drywall (0.75 in/s) and plaster (0.50 in/s) limit curves, drawn + continuous (constant-displacement bounds meeting the plateaus — no vertical + steps). Sized and positioned to match a Blastware report, measured off the + reference PDF. A technical breakdown of the curve is in + `docs/ri8507_compliance_curve.md`. + +- **Sensor self-check waveforms decoded and drawn (`minimateplus.sensor_check`).** + The "Sensor Check" traces Blastware shows to the right of the waveform panel + live in the series-3 binary's trailing block as four length-prefixed records + (`0x3c`–`0x3f`) using the same delta-block codec as the main waveform: + Tran/Vert/Long geophone ring-downs (the transducer's damped impulse response — + resonant frequency + overswing/damping) and a MicL pulse train (the mic's + known-signal gain check). `gather_report_data` decodes them from the retained + BW binary at report time; the report renders them as a strip flush against the + waveform panel plus the **Sensor Check → Frequency / Overswing Ratio** sub-rows + in the stats table. Verified against the reports on all 7 oracle events (mic + zero-crossing frequency = 20.1 Hz exact; geophone ring-downs consistent + ~7.5 Hz with overswing ~3.5). Tests in `tests/test_sensor_check.py`. + +- **Inspector tab in `seismo_lab.py` — annotated hex reader for series-3 + binaries (`minimateplus/binary_annotate.py`).** Tiles a raw Blastware file + into labeled spans (header / STRT / body record-chain / trailing metadata + + calibration + sensor-check records / footer) so a binary can be combed by eye. + +### Fixed +- **Event-report waveform panel — stacked-lane y-tick collision.** The lanes + touch, so each lane's bottom `-1.0` overprinted the next lane's top `1.0` at + the shared boundary. Prune the extreme ticks so each lane shows clean interior + ticks only. +- **Event-report header — serial+firmware line ran off the page.** The long + `BE##### V ##.##-#.## MiniMate Plus` string overflowed the right margin; + tighter right-column indent + BW's slightly smaller header size so it fits. + +--- + ## v0.30.0 — 2026-09-12 **The series-4 correctness release** — the Thor / Micromate counterpart to diff --git a/docs/ri8507_compliance_curve.md b/docs/ri8507_compliance_curve.md new file mode 100644 index 0000000..befb1a3 --- /dev/null +++ b/docs/ri8507_compliance_curve.md @@ -0,0 +1,135 @@ +# USBM RI8507 / OSMRE Blasting Compliance Curve — Reference + +Reference for the **velocity-vs-frequency blasting compliance chart** Blastware +draws on its Event Report ("USBM RI8507 And OSMRE"), and how seismo-relay +reproduces it. Implemented in [`sfm/compliance.py`](../sfm/compliance.py); the +spectral (FFT) side lives in [`waveform_fft.py`](../waveform_fft.py). + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events, each +with a Blastware Event Report + FFT Report as ground truth. Curve values from +USBM RI8507 Appendix B and 30 CFR 816.67. + +--- + +## What it is + +Two closely-related sources for the same limit curve: + +- **USBM RI8507** — Bureau of Mines *Report of Investigations 8507* (Siskind + et al., 1980), *"Structure Response and Damage Produced by Ground Vibration + From Surface Mine Blasting."* The curve is **Figure B-1**, Appendix B + ("Alternative Blasting Level Criteria"), p.73–74. +- **OSMRE / OSM** — the Office of Surface Mining Reclamation and Enforcement + codified it as **30 CFR 816.67, Figure 1**. "CFR" = the U.S. Code of Federal + Regulations. Same curve, regulatory force. + +The chart plots each geophone channel's significant vibration cycles as +`(frequency, peak velocity)` points against this limit. A point **below** the +line passes; **above** fails. + +--- + +## The limit curve + +A structure has a resonance band (~4–12 Hz for whole structures) where it is +most vulnerable, so the safe velocity is **lower** at those frequencies and +**higher** away from them. The curve captures this by alternating two kinds of +bound: + +- **Constant-velocity** segments — a flat horizontal line at a fixed PPV. +- **Constant-displacement** segments — a fixed peak *displacement* `d`. For + simple harmonic motion, peak velocity `v = 2πf·d`, so on a velocity-vs- + frequency **log-log** plot this is a straight line of slope +1 (velocity rises + with frequency). This is why the low- and high-frequency bounds are sloped. + +### Two lines — structure type + +RI8507 gives two lines for two interior-wall constructions (Table 13, p.67): + +| line | construction | plateau PPV | +|---|---|---| +| **Drywall** (solid) | modern gypsum wallboard | **0.75 in/s** | +| **Plaster** (dashed) | older plaster on wood lath | **0.50 in/s** | + +Plaster-on-lath is more damage-prone, hence the lower limit. You apply **one** +line depending on the monitored structure. + +### The four segments (Figure B-1, p.74) + +Going low → high frequency, each line is: + +1. **Ultimate low-frequency bound** — constant displacement **0.030 in** + (`v = 2πf·0.030`). Only relevant below ~4 Hz. +2. **Plateau** — constant velocity **0.75** (Drywall) / **0.50** (plaster) in/s. +3. **Rising diagonal** — constant displacement **0.008 in** (`v = 2πf·0.008`), + climbing from the plateau up to the high-frequency cap. +4. **High-frequency cap** — constant velocity **2.0 in/s** above ~40 Hz. + +The segments are drawn **continuous**: each bound is used over the frequency +range where it is the binding (lowest) limit, and consecutive bounds meet where +they are equal — so there are no vertical steps. Transition frequencies come +straight from the values (`f = V / (2π·d)`): + +| transition | formula | Drywall | Plaster | +|---|---|---|---| +| 0.030 in → plateau | `V_mid / (2π·0.030)` | 3.98 Hz | 2.65 Hz | +| plateau → 0.008 in | `V_mid / (2π·0.008)` | 14.92 Hz | 9.95 Hz | +| 0.008 in → 2.0 in/s | `2.0 / (2π·0.008)` | 39.79 Hz | 39.79 Hz | + +Because both lines share the same **0.008 in** rising diagonal, above ~15 Hz +they lie on the *same* line (both reach 2.0 in/s at ~40 Hz) — RI8507's literal +construction merges them there. Blastware renders the dashed line as a separate +parallel diagonal, but that is cosmetic: above ~15 Hz both structure types carry +the identical limit, so compliance is unaffected. + +> ⚠ RI8507's *Table 13* is a simpler two-range criterion with a **sharp +> discontinuity at 40 Hz** (flat plateau, then a jump to 2.0). Figure B-1 is the +> **smoothed** version that adds the 0.008 in transition — that is the one drawn +> on reports and implemented here. + +--- + +## The compliance scatter (the points) + +The cloud is **not** the FFT spectrum. It is a per-cycle, time-domain measure by +the **zero-crossing method** (`channel_compliance_points`): + +- Split the channel's waveform at its zero crossings. +- Each half-cycle contributes one point: **frequency** `= 1 / (2 · half-period)` + (from the samples between the two crossings), **velocity** `= peak |amplitude|` + in that half-cycle. + +This yields ~90–110 points per channel, and — by construction — each channel's +**highest** point equals that channel's PPV. Verified against Blastware: the +cloud shape, density, and ceiling all match. + +### Why not the FFT? + +A broadband blast spreads its energy across many FFT bins, so no single bin +reaches the time-domain peak — the FFT amplitudes come out ~10× below the +compliance-chart velocities. The compliance chart is a *per-cycle peak* view; +the **FFT** is a separate analysis (Blastware's *FFT Report*), reproduced by +[`waveform_fft.py`](../waveform_fft.py) and used for the dominant-frequency +readout and the #10 FFT view — not for this scatter. + +--- + +## Implementation + +- `sfm/compliance.py` + - `limit_at(freq, curve)` — the limit PPV at a frequency (`curve` = `"Drywall"` + or `"Plaster"`); curves are data in `_CURVES`, so more standards can be added. + - `channel_compliance_points(samples, sps)` — the zero-crossing scatter. + - `draw_compliance_chart(ax, channels, sps)` — matplotlib rendering (both + limit lines + per-channel scatter, Blastware's tick scales and channel + markers: Tran `+` red, Vert `×` green, Long `o` blue). +- Tests: `tests/test_compliance.py`. + +--- + +## Sources + +- USBM **RI8507** (Siskind, Stagg, Kopp, Dowding, 1980), Appendix B / Figure B-1, + p.73–74; Table 13, p.67. (`ref-stuff/usbm-ri8507-ground_vibration.pdf`.) +- **30 CFR 816.67**, "Use of explosives: Control of adverse effects," Figure 1 — + diff --git a/minimateplus/binary_annotate.py b/minimateplus/binary_annotate.py new file mode 100644 index 0000000..25e164d --- /dev/null +++ b/minimateplus/binary_annotate.py @@ -0,0 +1,75 @@ +"""Structural annotation of a Series-3 Blastware waveform binary. + +Pure, no I/O: takes the raw file bytes and returns a flat, gap-free tiling of +labelled :class:`Span` regions for a hex viewer to paint. Every byte is +covered — anything the decoder can't account for becomes an ``unknown`` span, +so undecoded regions (e.g. a stored spectral/FFT block, if one exists) stand +out instead of hiding. + +File layout (see ``blastware_file.py``): ``[header][21B STRT][body][26B footer]``. +The body is the record chain walked by :func:`waveform_codec.walk_records`. +""" +from __future__ import annotations + +from dataclasses import dataclass +from typing import List + +from .waveform_codec import walk_records + +_STRT_LEN = 21 +_FOOTER_LEN = 26 + + +@dataclass +class Span: + start: int # inclusive byte offset + end: int # exclusive byte offset + label: str # human-readable description + kind: str # 'header' | 'strt' | 'sample' | 'footer' | 'unknown' + + +def _tile(known: List[Span], total: int) -> List[Span]: + """Sort *known* spans and fill every gap with an ``unknown`` span, so the + result is a contiguous, non-overlapping tiling of ``[0, total)``. Overlaps + are resolved by clamping to the running position (first writer wins).""" + out: List[Span] = [] + pos = 0 + for s in sorted(known, key=lambda x: (x.start, x.end)): + if s.end <= pos: + continue # fully behind — dropped overlap + start = max(s.start, pos) + if start > pos: + out.append(Span(pos, start, "unknown", "unknown")) + out.append(s if start == s.start else Span(start, s.end, s.label, s.kind)) + pos = s.end + if pos < total: + out.append(Span(pos, total, "unknown", "unknown")) + return out + + +def annotate_blastware_binary(raw: bytes) -> List[Span]: + """Annotate a Series-3 waveform binary into a gap-free list of spans.""" + total = len(raw) + strt_pos = raw.find(b"STRT") + if strt_pos < 0: + return [Span(0, total, "unrecognized — no STRT record", "unknown")] + + known: List[Span] = [] + if strt_pos > 0: + known.append(Span(0, strt_pos, "File header", "header")) + known.append(Span(strt_pos, strt_pos + _STRT_LEN, "STRT record", "strt")) + + body_start = strt_pos + _STRT_LEN + footer_start = total - _FOOTER_LEN + if footer_start >= body_start: + known.append(Span(footer_start, total, "File footer", "footer")) + else: + footer_start = total # file too short for a footer + + body = raw[body_start:footer_start] + for rec in walk_records(body): + hi, lo = rec["mode"] + label = f"{rec['channel']} record (seg {rec['segment_index']}, mode {hi:02x} {lo:02x})" + known.append(Span(body_start + rec["offset"], body_start + rec["end"], label, "sample")) + + return _tile(known, total) diff --git a/minimateplus/sensor_check.py b/minimateplus/sensor_check.py new file mode 100644 index 0000000..7432ea8 --- /dev/null +++ b/minimateplus/sensor_check.py @@ -0,0 +1,146 @@ +r"""Decode the Blastware sensor self-check waveforms from a series-3 event binary. + +Reverse-engineered 2026-09-15 against 7 BE12844 (MiniMate Plus) oracle events. +After the main waveform record-chain and the trailing metadata / per-channel +calibration records, the binary carries four length-prefixed records tagged +0x3c-0x3f: the sensor self-check traces the unit records when it pulses each +sensor before monitoring. Blastware draws these as the little waveforms in the +"Sensor Check" strip on the right of the Event Report. + + * 0x3c / 0x3d / 0x3e = Tran / Vert / Long geophone ring-downs (a damped + oscillation at the geophone's resonance, ~7-8 Hz at 1024 sps). + * 0x3f = MicL, a pulse train at the mic self-test frequency + (~20 Hz), whose zero-crossing frequency is BW's mic "Channel Test" freq. + +Record framing (per record, all four chained by their length prefix):: + + [len:2 BE][id:1][00 00][Nchan:1][12-byte header][delta stream][40 02][6B] + \_________________ payload (len bytes) _______________________________/ + +The delta stream is ``payload[20 : len-8]`` (the ``40 02`` terminator sits at +``len-8``, followed by 6 trailing bytes). It uses the exact same 10/20/30/00 +delta-block tags as the main waveform codec +(:mod:`minimateplus.waveform_codec`), decoded here from an implicit anchor of 0 +— so the traces come out in the same 16-count raw units as the main waveform +(LSB = 0.005 in/s at Normal range for the geophones). +""" +from __future__ import annotations + +from typing import Dict, List + +from minimateplus.waveform_codec import walk_body + +# Record id → channel. Order mirrors the trailing per-channel calibration +# records (Tran / Vert / Long / MicL), confirmed against BW's sensor-check +# frequencies on all 7 oracle events. +_ID_TO_CHANNEL = {0x3C: "Tran", 0x3D: "Vert", 0x3E: "Long", 0x3F: "MicL"} +_CHAIN_IDS = (0x3C, 0x3D, 0x3E, 0x3F) + +_HEADER_LEN = 20 # payload bytes before the delta stream +_TRAILER_LEN = 8 # 40 02 terminator + 6 trailing bytes after the stream + + +def _s4(nib: int) -> int: + """Sign-extend a 4-bit nibble delta.""" + return nib - 16 if nib >= 8 else nib + + +def _i8(byte: int) -> int: + """Sign-extend an 8-bit int delta.""" + return byte - 256 if byte >= 128 else byte + + +def _decode_delta_stream(buf: bytes) -> List[int]: + """Accumulate a 10/20/30/00 delta-block stream from an anchor of 0, + stopping at the 0x40 terminator. + + Mirrors the block semantics in + :func:`minimateplus.waveform_codec.decode_waveform_v2` (fully decoded & + byte-exact as of 2026-05-11); see that module for the format details. + """ + out: List[int] = [] + cur = 0 + for blk in walk_body(buf, 0): + fam = blk.tag_hi & 0xF0 + if fam == 0x10: + # nibble deltas, high nibble first + for byte in blk.data: + for nib in ((byte >> 4) & 0xF, byte & 0xF): + cur += _s4(nib) + out.append(cur) + elif fam == 0x20: + # int8 deltas + for byte in blk.data: + cur += _i8(byte) + out.append(cur) + elif fam == 0x30: + # 12-bit signed deltas, packed as tag_lo/4 groups of 6 bytes + for g in range(blk.tag_lo // 4): + grp = blk.data[g * 6:(g + 1) * 6] + if len(grp) < 6: + break + high_word = (grp[0] << 8) | grp[1] + for k in range(4): + nib = (high_word >> (12 - 4 * k)) & 0xF + v = (nib << 8) | grp[2 + k] + if v >= 0x800: + v -= 0x1000 + cur += v + out.append(cur) + elif fam == 0x00: + # RLE zero-delta run (wide form carries the high nibble in the tag) + run = ((blk.tag_hi & 0x0F) << 8) | blk.tag_lo + out.extend([cur] * run) + elif fam == 0x40: + # segment / record terminator + break + return out + + +def _find_chain(body: bytes): + """Locate the four length-prefixed sensor-check records. + + Returns a list of ``(offset, id, length)`` or ``None``. The chain is + validated by walking the ids 0x3c → 0x3d → 0x3e → 0x3f via their own length + prefixes, so a stray 0x3c byte in the waveform data cannot match. + """ + for p in range(len(body) - 6): + if body[p + 2] == 0x3C and body[p + 3] == 0 and body[p + 4] == 0: + q = p + recs = [] + ok = True + for expect in _CHAIN_IDS: + if q + 3 > len(body) or body[q + 2] != expect: + ok = False + break + length = int.from_bytes(body[q:q + 2], "big") + recs.append((q, expect, length)) + q = q + 2 + length + if ok and len(recs) == 4: + return recs + return None + + +def decode_sensor_check(raw: bytes) -> Dict[str, List[int]]: + """Decode the four sensor self-check traces from a series-3 event binary. + + Returns ``{"Tran": [...], "Vert": [...], "Long": [...], "MicL": [...]}`` in + raw decode units (same 16-count LSB as the main waveform), or ``{}`` if the + binary carries no sensor-check block (a histogram event, a non-series-3 + file, or a unit/firmware that doesn't store it). + """ + strt = raw.find(b"STRT") + if strt < 0 or len(raw) < strt + 21 + 26: + return {} + body = raw[strt + 21: len(raw) - 26] + chain = _find_chain(body) + if not chain: + return {} + out: Dict[str, List[int]] = {} + for off, rid, length in chain: + payload = body[off + 2: off + 2 + length] + if len(payload) < _HEADER_LEN + _TRAILER_LEN: + continue + stream = payload[_HEADER_LEN: length - _TRAILER_LEN] + out[_ID_TO_CHANNEL[rid]] = _decode_delta_stream(stream) + return out diff --git a/seismo_lab.py b/seismo_lab.py index 1986127..5cad795 100644 --- a/seismo_lab.py +++ b/seismo_lab.py @@ -54,6 +54,7 @@ from s3_analyzer import ( # noqa: E402 write_claude_export, ) from frame_db import FrameDB # noqa: E402 +from minimateplus.binary_annotate import annotate_blastware_binary # noqa: E402 # ── colour palette ──────────────────────────────────────────────────────────── BG = "#1e1e1e" @@ -2675,6 +2676,95 @@ class DownloadPanel(tk.Frame): self._on_capture_ready(bw_path, s3_path, label) +# ───────────────────────────────────────────────────────────────────────────── +# Inspector panel — annotated hex view of a Series-3 binary +# ───────────────────────────────────────────────────────────────────────────── + +class InspectorPanel(tk.Frame): + """Load any Series-3 waveform binary and read it as an annotated hex dump. + + Regions the decoder understands (header, STRT, per-channel sample records, + footer) are labelled and colour-coded; everything the decoder cannot account + for is flagged UNKNOWN, so undecoded bytes stand out for hand-inspection. + """ + + _KIND_COLOR = { + "header": ACCENT, + "strt": YELLOW, + "sample": COL_S3, + "footer": FG_DIM, + "unknown": RED, + } + + def __init__(self, parent: tk.Widget, initialdir=None, **kw) -> None: + super().__init__(parent, bg=BG, **kw) + self._path = None + self._initialdir = initialdir + self._build() + + def _build(self) -> None: + bar = tk.Frame(self, bg=BG2) + bar.pack(side=tk.TOP, fill=tk.X) + tk.Button(bar, text="Open binary…", command=self._open, bg=BG3, fg=FG, + relief=tk.FLAT, font=MONO, activebackground=ACCENT).pack(side=tk.LEFT, padx=6, pady=6) + self._path_var = tk.StringVar(value="(no file loaded)") + tk.Label(bar, textvariable=self._path_var, bg=BG2, fg=FG_DIM, font=MONO).pack(side=tk.LEFT, padx=6) + self._summary_var = tk.StringVar(value="") + tk.Label(bar, textvariable=self._summary_var, bg=BG2, fg=FG, font=MONO).pack(side=tk.RIGHT, padx=10) + + legend = tk.Frame(self, bg=BG2) + legend.pack(side=tk.TOP, fill=tk.X) + tk.Label(legend, text="legend:", bg=BG2, fg=FG_DIM, font=MONO).pack(side=tk.LEFT, padx=(8, 2)) + for kind, color in self._KIND_COLOR.items(): + tk.Label(legend, text=f"■ {kind}", bg=BG2, fg=color, font=MONO).pack(side=tk.LEFT, padx=5, pady=2) + + self._text = scrolledtext.ScrolledText( + self, bg=BG, fg=FG, insertbackground=FG, font=MONO, wrap=tk.NONE, borderwidth=0) + self._text.pack(side=tk.TOP, fill=tk.BOTH, expand=True) + for kind, color in self._KIND_COLOR.items(): + self._text.tag_configure(kind, foreground=color) + self._text.tag_configure("label", foreground="#ffffff", font=("Consolas", 9, "bold")) + self._text.tag_configure("dim", foreground=FG_DIM) + self._text.configure(state=tk.DISABLED) + + def _open(self) -> None: + p = filedialog.askopenfilename(title="Open a Series-3 binary", initialdir=self._initialdir) + if p: + self.load(Path(p)) + + def load(self, path: Path) -> None: + try: + raw = path.read_bytes() + spans = annotate_blastware_binary(raw) + except Exception as e: # noqa: BLE001 — surface any read/annotate failure to the user + messagebox.showerror("Inspector", f"Failed to read/annotate:\n{path}\n\n{e}") + return + self._path = path + self._path_var.set(str(path)) + self._render(raw, spans) + + def _render(self, raw: bytes, spans) -> None: + t = self._text + t.configure(state=tk.NORMAL) + t.delete("1.0", tk.END) + unknown = sum(s.end - s.start for s in spans if s.kind == "unknown") + pct = 100 * unknown / max(1, len(raw)) + self._summary_var.set(f"{len(raw)} B · {len(spans)} regions · {pct:.1f}% unknown") + for s in spans: + t.insert(tk.END, f"\n── {s.label} [0x{s.start:04x}:0x{s.end:04x}] {s.end - s.start} B ──\n", ("label",)) + self._insert_hex(t, raw, s.start, s.end, s.kind) + t.configure(state=tk.DISABLED) + + def _insert_hex(self, t: tk.Text, raw: bytes, start: int, end: int, kind: str) -> None: + for off in range(start, end, 16): + row = raw[off:min(off + 16, end)] + hx = " ".join(f"{b:02x}" for b in row).ljust(16 * 3 - 1) + txt = "".join(chr(b) if 32 <= b < 127 else "." for b in row) + t.insert(tk.END, f" 0x{off:04x} ", ("dim",)) + t.insert(tk.END, hx, (kind,)) + t.insert(tk.END, f" {txt}\n", ("dim",)) + + # ───────────────────────────────────────────────────────────────────────────── # Main application window # ───────────────────────────────────────────────────────────────────────────── @@ -2730,6 +2820,9 @@ class SeismoLab(tk.Tk): ) nb.add(self._download_panel, text=" Download ") + self._inspector_panel = InspectorPanel(nb) + nb.add(self._inspector_panel, text=" Inspector ") + self._nb = nb self.protocol("WM_DELETE_WINDOW", self._on_close) diff --git a/sfm/compliance.py b/sfm/compliance.py new file mode 100644 index 0000000..bdc4692 --- /dev/null +++ b/sfm/compliance.py @@ -0,0 +1,133 @@ +"""USBM RI8507 / OSMRE blasting compliance chart. + +Renders the velocity-vs-frequency compliance scatter Blastware draws on its Event +Report: each channel's significant waveform cycles as ``(frequency, peak +velocity)`` points on log-log axes against the regulatory limit curve(s). A point +below the curve passes; above fails. + +Two pieces, kept separate so both can be reused/extended: + * ``limit_at`` / ``limit_curve`` — the regulatory limit curve(s), as data. + * ``channel_compliance_points`` — the per-cycle (freq, velocity) scatter, by + the zero-crossing method (matches Blastware: each channel's cloud tops out + at that channel's PPV). + +Limit curves (USBM RI8507 Figure B-1 / OSM 30 CFR 816.67), drawn CONTINUOUS — a +constant-displacement bound (sloped, ``v = 2πf·d``) meets a constant-velocity +plateau at the frequency where they're equal, so there are no vertical steps +(matching how Blastware draws it). Two lines: + * **Drywall** (modern gypsum board) — 0.75 in/s plateau (solid). + * **Plaster** on wood lath (older homes) — 0.50 in/s plateau (dashed). +Both use a 0.030 in low-frequency displacement bound and rise through a 0.010 in +displacement bound to a 2.0 in/s high-frequency plateau. Values from USBM RI8507 +(Appendix B) / 30 CFR 816.67; ⚠ confirm the exact shape against a Blastware +report before trusting for compliance. +""" +from __future__ import annotations + +import math +from typing import Dict, Sequence, Tuple + +import numpy as np +from matplotlib.ticker import FixedLocator, NullLocator + +# curve name → (low-freq "ultimate" displacement in, mid velocity plateau in/s, +# high-freq displacement in, high-freq velocity plateau in/s). +# RI8507 Fig B-1 (p.74): ultimate max displacement 0.030 in (< ~4 Hz), plateau +# 0.75 (Drywall) / 0.50 (plaster), rising diagonal at 0.008 in displacement up to +# a 2.0 in/s plateau reached at ~40 Hz. +_CURVES: Dict[str, Tuple[float, float, float, float]] = { + "Drywall": (0.030, 0.75, 0.008, 2.00), + "Plaster": (0.030, 0.50, 0.008, 2.00), +} +# how each curve is stroked on the chart +_CURVE_STYLE = {"Drywall": {"ls": "-", "lw": 1.0}, "Plaster": {"ls": "--", "lw": 0.9}} + +STANDARDS = tuple(_CURVES) + +# Blastware's channel markers/colours on the compliance chart. +_CHANNEL_STYLE = { + "Tran": ("+", "#d62728"), # red + + "Vert": ("x", "#2ca02c"), # green x + "Long": ("o", "#1f77b4"), # blue o +} + + +def limit_at(freq_hz: float, curve: str = "Drywall") -> float: + """Max allowed PPV (in/s) at ``freq_hz`` for ``curve`` (continuous).""" + d_low, v_mid, d_high, v_high = _CURVES[curve] + f = max(freq_hz, 1.0) + f_a = v_mid / (2.0 * math.pi * d_low) # disp_low → vel_mid + f_b = v_mid / (2.0 * math.pi * d_high) # vel_mid → disp_high + f_c = v_high / (2.0 * math.pi * d_high) # disp_high → vel_high + if f <= f_a: + return 2.0 * math.pi * f * d_low + if f <= f_b: + return v_mid + if f <= f_c: + return 2.0 * math.pi * f * d_high + return v_high + + +def limit_curve(curve: str = "Drywall", fmin: float = 1.0, fmax: float = 100.0, n: int = 400): + """(freqs, limits) sampled across the band for plotting one curve.""" + freqs = np.logspace(np.log10(fmin), np.log10(fmax), n) + return freqs, np.array([limit_at(f, curve) for f in freqs]) + + +def channel_compliance_points( + samples: Sequence[float], sps: float, fmin: float = 1.0, fmax: float = 100.0, + vmin: float = 0.0, +) -> Tuple[np.ndarray, np.ndarray]: + """Per-cycle (frequency, peak velocity) scatter for one channel. + + Zero-crossing method: split the trace at sign changes; each half-cycle + contributes one point at ``(1/(2·half_period), max|amplitude|)``. Matches + Blastware — the cloud's ceiling is the channel PPV. ``samples`` must be in the + velocity unit you want plotted (in/s). Points outside ``[fmin, fmax]`` or at + or below ``vmin`` are dropped. + """ + x = np.asarray(samples, dtype=float) + if x.size < 3: + return np.empty(0), np.empty(0) + zc = np.where(np.diff(np.signbit(x)))[0] + freqs, vels = [], [] + for a, b in zip(zc[:-1], zc[1:]): + half_period = (b - a) / sps + if half_period <= 0: + continue + freqs.append(1.0 / (2.0 * half_period)) + vels.append(float(np.abs(x[a:b + 1]).max())) + f = np.array(freqs) + v = np.array(vels) + keep = (f >= fmin) & (f <= fmax) & (v > vmin) + return f[keep], v[keep] + + +def draw_compliance_chart(ax, channels: Dict[str, Sequence[float]], sps: float) -> None: + """Draw the compliance chart (both limit curves + per-channel scatter).""" + for name, style in _CURVE_STYLE.items(): + cf, cv = limit_curve(name) + ax.plot(cf, cv, color="#333", zorder=3, **style) + + for ch, (marker, color) in _CHANNEL_STYLE.items(): + samples = channels.get(ch) + if samples is None or len(samples) == 0: + continue + f, v = channel_compliance_points(samples, sps) + ax.scatter(f, v, marker=marker, s=12, c=color, linewidths=0.7, zorder=4, label=ch) + + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_xlim(1, 100) + ax.set_ylim(0.0394, 10) + ax.set_box_aspect(1) # square plot box (log-log compliance charts are square) + xt = [1, 2, 5, 10, 20, 50, 100] + yt = [0.0394, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10] + ax.xaxis.set_major_locator(FixedLocator(xt)); ax.xaxis.set_minor_locator(NullLocator()) + ax.yaxis.set_major_locator(FixedLocator(yt)); ax.yaxis.set_minor_locator(NullLocator()) + ax.set_xticklabels([str(v) for v in xt]) + ax.set_yticklabels([("%g" % v) for v in yt]) + ax.set_xlabel("Frequency (Hz)", fontsize=7) + ax.set_ylabel("Velocity (in/s)", fontsize=7) + ax.tick_params(labelsize=6) + ax.grid(True, which="both", ls=":", lw=0.4, color="#ccc") diff --git a/sfm/report_pdf.py b/sfm/report_pdf.py index 60dca98..794c836 100644 --- a/sfm/report_pdf.py +++ b/sfm/report_pdf.py @@ -121,6 +121,11 @@ class ReportData: t0_ms: Optional[float] = None dt_ms: Optional[float] = None + # Sensor self-check traces — {ch: [samples]} in raw decode units, decoded + # from the binary's trailing block (see minimateplus.sensor_check). The + # little waveforms BW draws in its "Sensor Check" strip. Empty when absent. + sensor_check_waveforms: dict = field(default_factory=dict) + # Record-type discriminator record_type: Optional[str] = None is_histogram: bool = False @@ -246,6 +251,8 @@ def gather_report_data( "peak_accel_g": ch.get("peak_accel_g"), "peak_disp_in": ch.get("peak_disp_in"), "sensor_check": sc_ch.get("result"), + "sc_freq_hz": sc_ch.get("freq_hz"), + "sc_ratio": sc_ch.get("ratio"), "peak_date": peak_date, "peak_time": peak_time, }) @@ -290,6 +297,19 @@ def gather_report_data( except Exception as exc: log.warning("gather_report_data: hdf5 read failed: %s", exc) + # ── Sensor self-check traces — decoded from the retained raw binary ── + # The .h5 holds only the main waveform; the sensor-check traces live in the + # binary's trailing block, so decode them straight from the kept BW file. + # Waveform events only (histograms have no sensor-check strip). + if not rd.is_histogram: + try: + from minimateplus.sensor_check import decode_sensor_check + bw_path, _a5 = store.paths_for(serial, filename) + if bw_path.exists(): + rd.sensor_check_waveforms = decode_sensor_check(bw_path.read_bytes()) + except Exception as exc: + log.warning("gather_report_data: sensor-check decode failed: %s", exc) + # ── Histogram aggregation ── # Codec emits ~N per-block samples (typically 1/sec); BW reports # one bar per configured interval (1 min / 5 min / etc.). When @@ -396,9 +416,34 @@ def _render_waveform_layout(fig, rd: ReportData) -> None: ax_stats = fig.add_subplot(gs[2]); ax_stats.axis("off") _draw_channel_stats_waveform(ax_stats, rd) + _draw_compliance_panel(fig, rd) _draw_waveform_subplot(fig, gs[3], rd) +# Compliance-chart placement, in figure fractions. Measured directly off a +# Blastware Event Report PDF (ref-stuff/n844lqhbzt0w_bw_pdf.pdf) so the chart +# matches BW's size and position: it spans from just under the header down +# through the stats band, hard against the right page margin. The left edge +# leaves room for the y-axis tick labels + "Velocity (in/s)" title, which the +# compacted stats table (see _draw_channel_stats_waveform) is sized to clear. +_COMPLIANCE_BOX = (0.489, 0.502, 0.951, 0.867) # x0, y0, x1, y1 + + +def _draw_compliance_panel(fig, rd: ReportData) -> None: + """Large USBM RI8507 compliance chart in the upper-right, sized and + positioned to match Blastware's Event Report (see _COMPLIANCE_BOX).""" + x0, y0, x1, y1 = _COMPLIANCE_BOX + fig.text((x0 + x1) / 2, y1 + 0.006, "USBM RI8507 And OSMRE", fontsize=9, + weight="bold", color="#333", ha="center", va="bottom") + if rd.channels and rd.sample_rate_sps: + from sfm.compliance import draw_compliance_chart + ax = fig.add_axes([x0, y0, x1 - x0, y1 - y0]) + draw_compliance_chart(ax, rd.channels, rd.sample_rate_sps) + else: + fig.text((x0 + x1) / 2, (y0 + y1) / 2, "(no waveform data)", fontsize=8, + color="#bbb", ha="center", va="center", style="italic") + + def _render_histogram_layout(fig, rd: ReportData) -> None: """Histogram layout: header / mic-only / per-channel stats / bar plot. @@ -477,11 +522,11 @@ def _split_iso_to_date_time(iso: Optional[str]) -> tuple[Optional[str], Optional return (None, None) -def _kv(ax, x, y, label, value, *, label_w=0.18): +def _kv(ax, x, y, label, value, *, label_w=0.18, fontsize=8): """Render a 'Label Value' row at axes-coordinates (x, y).""" - ax.text(x, y, label, fontsize=8, color="#555", ha="left", va="top", + ax.text(x, y, label, fontsize=fontsize, color="#555", ha="left", va="top", transform=ax.transAxes) - ax.text(x + label_w, y, _fmt(value), fontsize=8, ha="left", va="top", + ax.text(x + label_w, y, _fmt(value), fontsize=fontsize, ha="left", va="top", transform=ax.transAxes, family="monospace") @@ -544,14 +589,17 @@ def _draw_header_columns(ax, rows_left, rd: ReportData) -> None: ("File Name", rd.file_name), ("Post Event Notes", rd.post_event_notes), ] + # fontsize 7.5 (BW's header is a touch smaller than our body text) + a + # tighter right-column value indent so the long serial+firmware line + # ("BE##### V ##.##-#.## MiniMate Plus") fits without running off the page. y = 0.95 dy = 0.095 for label, value in rows_left: - _kv(ax, 0.0, y, label, value, label_w=0.18) + _kv(ax, 0.0, y, label, value, label_w=0.18, fontsize=7.5) y -= dy y = 0.95 for label, value in rows_right: - _kv(ax, 0.55, y, label, value, label_w=0.20) + _kv(ax, 0.55, y, label, value, label_w=0.14, fontsize=7.5) y -= dy @@ -574,19 +622,14 @@ def _draw_mic_and_usbm(ax, rd: ReportData) -> None: transform=ax.transAxes, va="top") rows = _mic_rows(rd) y = 0.80 + # Tighter label indent + slightly smaller font so the long "Channel Test + # Passed (Freq = … Amp = … mv)" line clears the enlarged compliance chart's + # left edge (_COMPLIANCE_BOX) instead of running behind it. for label, value in rows: - _kv(ax, 0.0, y, label, value, label_w=0.18) + _kv(ax, 0.0, y, label, value, label_w=0.13, fontsize=7) y -= 0.15 - - # USBM chart placeholder — upper-right. Real piecewise compliance - # curves are a separate work item; for now this just shows the title - # + a "see report" message so the layout is correct. - ax.text(0.72, 0.97, "USBM RI8507 And OSMRE", - fontsize=9, weight="bold", color="#333", ha="center", va="top", - transform=ax.transAxes) - ax.text(0.72, 0.50, "[compliance chart\ncoming soon]", - fontsize=8, color="#bbb", ha="center", va="center", - transform=ax.transAxes, style="italic") + # The USBM compliance chart is drawn as its own large square panel spanning + # the mic + stats rows on the right — see _draw_compliance_panel(). def _mic_rows(rd: ReportData) -> list[tuple[str, Optional[str]]]: @@ -636,8 +679,18 @@ def _draw_channel_stats_waveform(ax, rd: ReportData) -> None: ("Peak Acceleration", "peak_accel_g", "g"), ("Peak Displacement", "peak_disp_in", "in"), ("Sensor Check", "sensor_check", ""), + # Sensor-check sub-rows (indented under "Sensor Check", like BW): the + # geophone ring-down frequency + overswing ratio from the self-check. + (" Frequency", "sc_freq_hz", "Hz"), + (" Overswing Ratio", "sc_ratio", ""), ] - _draw_stats_table(ax, rd, rows_spec) + # Compacted to the left half so the enlarged compliance chart (BW-sized, + # right against the page margin) has room — see _COMPLIANCE_BOX. + _draw_stats_table( + ax, rd, rows_spec, + bbox_width=0.42, fontsize=7.5, + col_widths=[0.185, 0.065, 0.065, 0.065, 0.040], + ) _draw_pvs_summary(ax, rd, n_data_rows=len(rows_spec)) @@ -698,19 +751,39 @@ def _draw_pvs_summary( table_bottom_y = getattr(ax, "_stats_table_bottom", -0.10) pvs_y = table_bottom_y - 0.04 # small gap below the table border - # Centered for visual balance — looks intentional rather than offset. - # The original BW-replica had a "NA: Not Applicable" caption below - # this line; dropped because we use "—" for missing values and the - # legend was always squished against the PVS line. - ax.text(0.5, pvs_y, line, fontsize=9, weight="bold", - ha="center", va="top", transform=ax.transAxes) + # Centered under the stats table for visual balance — looks intentional + # rather than offset. When the table is compacted (waveform layout), it + # occupies only the left portion of the axes, so center on the table's + # width rather than the full axes (which would push the line under the + # compliance chart). The original BW-replica had a "NA: Not Applicable" + # caption below this line; dropped because we use "—" for missing values. + table_w = getattr(ax, "_stats_table_width", 0.80) + if table_w < 0.79: + # Compacted (waveform) layout: left-align under the table, one point + # smaller, so the line clears the enlarged compliance chart's + # bottom-left tick labels on the right. + ax.text(0.0, pvs_y, line, fontsize=8, weight="bold", + ha="left", va="top", transform=ax.transAxes) + else: + ax.text(0.5, pvs_y, line, fontsize=9, weight="bold", + ha="center", va="top", transform=ax.transAxes) -def _draw_stats_table(ax, rd: ReportData, rows_spec: list[tuple[str, str, str]]) -> None: +def _draw_stats_table( + ax, rd: ReportData, rows_spec: list[tuple[str, str, str]], + *, bbox_width: float = 0.80, fontsize: float = 8, + col_widths: Optional[list[float]] = None, +) -> None: """Render a per-channel stats table (Tran/Vert/Long). rows_spec: list of (label, field_name_in_channel_stats, unit_string) + + ``bbox_width`` / ``col_widths`` / ``fontsize`` let a caller compact the + table (the waveform layout packs it into the left half to clear the + compliance chart; the histogram layout keeps the wider defaults). """ + if col_widths is None: + col_widths = [0.28, 0.14, 0.14, 0.14, 0.10] headers = ["", "Tran", "Vert", "Long", ""] ch_lookup = {c["name"]: c for c in rd.channel_stats} @@ -726,6 +799,8 @@ def _draw_stats_table(ax, rd: ReportData, rows_spec: list[tuple[str, str, str]]) if field == "zc_freq_hz": prefix = ">" if ch_rec.get("zc_freq_above_range") else "" return f"{prefix}{val:.0f}" + if field in ("sc_freq_hz", "sc_ratio"): + return f"{val:.1f}" # BW shows 1 decimal (7.5 Hz, 3.6) return f"{val:.3f}" return str(val) @@ -750,16 +825,17 @@ def _draw_stats_table(ax, rd: ReportData, rows_spec: list[tuple[str, str, str]]) table_bottom = 1.0 - table_height tbl = ax.table( cellText=table_data, - colWidths=[0.28, 0.14, 0.14, 0.14, 0.10], + colWidths=col_widths, cellLoc="left", edges="open", - bbox=[0.0, table_bottom, 0.80, table_height], + bbox=[0.0, table_bottom, bbox_width, table_height], ) tbl.auto_set_font_size(False) - tbl.set_fontsize(8) + tbl.set_fontsize(fontsize) for j in range(5): tbl[(0, j)].set_text_props(weight="bold", color="#555") - # Stash the bottom Y so _draw_pvs_summary can position itself below. + # Stash the bottom Y + width so _draw_pvs_summary can position itself. ax._stats_table_bottom = table_bottom + ax._stats_table_width = bbox_width def _channel_axis_color(ch: str) -> str: @@ -769,9 +845,25 @@ def _channel_axis_color(ch: str) -> str: def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None: """4-channel stacked waveform plot — Instantel printout order (MicL on top, Tran on bottom), shared x-axis in SECONDS, trigger - triangle markers at t=0, '0.0' baseline label on right of each.""" - inner = gridspec_cell.subgridspec(4, 1, hspace=0.0) + triangle markers at t=0, '0.0' baseline label on right of each. + + When sensor self-check traces are present (rd.sensor_check_waveforms), a + narrow "Sensor Check" strip of per-channel mini-plots is drawn to the right, + aligned to the lanes — matching Blastware's Event Report. + """ + from matplotlib.ticker import MaxNLocator + order = ["MicL", "Long", "Vert", "Tran"] + has_sc = bool(rd.sensor_check_waveforms) + if has_sc: + # main lanes + a narrow sensor-check strip column, flush against the + # main panel (BW shares the border — no gap), with the "0.0" baseline + # labels moved to the right of the strip. Proportions match BW's + # Event Report (main ~0.75 / strip ~0.10 of the panel width). + inner = gridspec_cell.subgridspec(4, 2, width_ratios=[1.0, 0.13], + wspace=0.0, hspace=0.0) + else: + inner = gridspec_cell.subgridspec(4, 1, hspace=0.0) sr = rd.sample_rate_sps or 1024 # Convert ms-based time axis to seconds for the x-axis dt_s = (rd.dt_ms or (1000.0 / sr)) / 1000.0 @@ -790,9 +882,12 @@ def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None: _geo_amax = _a geo_shared = max(_geo_amax * 1.10, GEO_FLOOR_INS) + main_axes = [] + sc_axes = [] last_idx = len(order) - 1 for i, ch in enumerate(order): - ax = fig.add_subplot(inner[i]) + ax = fig.add_subplot(inner[i, 0] if has_sc else inner[i]) + main_axes.append(ax) values = rd.channels.get(ch) or [] times = [t0_s + j * dt_s for j in range(len(values))] @@ -810,9 +905,12 @@ def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None: # Channel label on the LEFT (matches BW) ax.set_ylabel(ch, fontsize=8, rotation=0, ha="right", va="center", color=_channel_axis_color(ch), weight="bold", labelpad=14) - # "0.0" on the RIGHT (BW convention) - ax.text(1.005, 0.5, "0.0", transform=ax.transAxes, - fontsize=7, color="#555", va="center", ha="left") + # "0.0" baseline label on the RIGHT (BW convention). With the sensor- + # check strip attached, it goes to the right of the STRIP (drawn below); + # otherwise just outside the main lane. + if not has_sc: + ax.text(1.005, 0.5, "0.0", transform=ax.transAxes, + fontsize=7, color="#555", va="center", ha="left") ax.grid(True, linestyle="--", linewidth=0.3, color="#bbb", alpha=0.6) # Vertical dashed trigger line at t=0 @@ -827,12 +925,46 @@ def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None: else: ax.tick_params(axis="x", labelsize=7) ax.tick_params(axis="y", labelsize=6) + # Stacked lanes touch, so the top/bottom y-tick labels of adjacent lanes + # would overprint at the shared boundary. Prune the extreme ticks so + # each boundary shows clean interior ticks (0.5 / 0.0 / -0.5) only. + ax.yaxis.set_major_locator(MaxNLocator(nbins=4, prune="both")) + + # Sensor self-check mini-plot in the right strip (aligned to this lane). + if has_sc: + scx = fig.add_subplot(inner[i, 1]) + sc_axes.append(scx) + sc_vals = rd.sensor_check_waveforms.get(ch) or [] + if sc_vals: + _col = _channel_axis_color(ch) + # Faint zero baseline (BW draws the channel baseline through the + # strip) — reference for the one-sided geophone ring-downs. + scx.axhline(0.0, color=_col, linewidth=0.3, alpha=0.4) + scx.plot(range(len(sc_vals)), sc_vals, color=_col, linewidth=0.5) + # Fit the trace to the box (BW-style) rather than a symmetric + # scale: the geo self-checks are one-sided dips, so a symmetric + # scale would strand them in the bottom half with an empty top. + _lo, _hi = min(sc_vals), max(sc_vals) + _pad = 0.10 * ((_hi - _lo) or 1.0) + scx.set_ylim(_lo - _pad, _hi + _pad) + scx.set_xticks([]); scx.set_yticks([]) + for _s in scx.spines.values(): + _s.set_linewidth(0.4); _s.set_color("#999") + # "0.0" baseline label to the RIGHT of the strip (BW convention) + scx.text(1.10, 0.5, "0.0", transform=scx.transAxes, + fontsize=7, color="#555", va="center", ha="left") # Trigger triangle marker ▼ above the top channel at t=0 - top_ax = fig.axes[-4] # MicL is the first added in this gridspec + top_ax = main_axes[0] # MicL top_ax.plot([0], [top_ax.get_ylim()[1]], marker="v", color="black", markersize=8, clip_on=False, zorder=10) + # "Sensor Check" caption under the strip (BW convention) + if has_sc and sc_axes: + pos = sc_axes[-1].get_position() + fig.text((pos.x0 + pos.x1) / 2, pos.y0 - 0.012, "Sensor Check", + fontsize=7, color="#555", ha="center", va="top") + # Compute scale-per-division for the footer (10 divs across the chart) # and find peak geo amplitude for the geo amp/div setting. total_s = times[-1] - times[0] if values else 0 diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W b/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W new file mode 100644 index 0000000..785721f Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W b/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W new file mode 100644 index 0000000..b4395a1 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W b/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W new file mode 100644 index 0000000..6c7afb6 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W b/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W new file mode 100644 index 0000000..fd9a7cd Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W b/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W new file mode 100644 index 0000000..8726fdb Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W b/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W new file mode 100644 index 0000000..9e4fa8f Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W b/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W new file mode 100644 index 0000000..6e9ef75 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W differ diff --git a/tests/test_binary_annotate.py b/tests/test_binary_annotate.py new file mode 100644 index 0000000..3748041 --- /dev/null +++ b/tests/test_binary_annotate.py @@ -0,0 +1,50 @@ +"""Structural annotation of a Series-3 Blastware binary (for the seismo_lab +Binary Inspector). The annotator maps byte ranges to labelled spans; anything +the decoder can't account for is a first-class ``unknown`` span, so the whole +file is tiled and the gaps (candidate FFT/spectral data) are visible. +""" +from pathlib import Path + +from minimateplus.binary_annotate import annotate_blastware_binary, Span + +# A known-good full-3-channel Series-3 waveform binary (the V70 cracking fixture). +FIXTURE = Path(__file__).parent / "fixtures" / "5-11-26" / "M529LL1L.V70" + + +def _raw() -> bytes: + return FIXTURE.read_bytes() + + +def test_spans_tile_the_whole_file(): + raw = _raw() + spans = annotate_blastware_binary(raw) + assert spans, "expected at least one span" + assert spans[0].start == 0 + assert spans[-1].end == len(raw) + for a, b in zip(spans, spans[1:]): + assert a.end == b.start, f"gap/overlap between {a!r} and {b!r}" + for s in spans: + assert s.start < s.end, f"empty/negative span {s!r}" + + +def test_strt_record_is_located(): + raw = _raw() + spans = annotate_blastware_binary(raw) + strt = [s for s in spans if s.kind == "strt"] + assert strt, "expected a STRT region" + assert raw[strt[0].start : strt[0].start + 4] == b"STRT" + + +def test_geo_sample_records_annotated(): + raw = _raw() + spans = annotate_blastware_binary(raw) + chans = {s.label.split()[0] for s in spans if s.kind == "sample"} + # V70 is a full three-geo-channel event. + assert {"Tran", "Vert", "Long"} <= chans, f"expected geo records, got {chans}" + + +def test_footer_is_last(): + raw = _raw() + spans = annotate_blastware_binary(raw) + assert spans[-1].kind == "footer" + assert spans[-1].end - spans[-1].start == 26 diff --git a/tests/test_compliance.py b/tests/test_compliance.py new file mode 100644 index 0000000..c102e9c --- /dev/null +++ b/tests/test_compliance.py @@ -0,0 +1,35 @@ +"""USBM/OSMRE compliance curve + scatter logic (sfm.compliance). +Rendering is verified visually against Blastware reports.""" +import math + +import numpy as np + +from sfm.compliance import limit_at, channel_compliance_points + + +def test_osmre_velocity_segments(): + assert abs(limit_at(6.0) - 0.75) < 1e-9 # 3.5–12 Hz flat + assert abs(limit_at(50.0) - 2.00) < 1e-9 # 30–100 Hz flat + + +def test_displacement_segments(): + assert abs(limit_at(2.0) - 2 * math.pi * 2.0 * 0.030) < 1e-9 # low-freq 0.030 in + assert abs(limit_at(20.0) - 2 * math.pi * 20.0 * 0.008) < 1e-9 # rising diagonal 0.008 in + + +def test_limit_clamps_below_1hz(): + assert limit_at(0.1) == limit_at(1.0) + + +def test_scatter_ceiling_is_ppv_at_dominant_freq(): + # ~27 Hz blast-like trace whose energy peaks mid-record (inside full cycles, + # as a real event does): the scatter cloud's ceiling is the trace PPV and the + # top point sits near the dominant frequency. + sps, n = 1024.0, 3328 + t = np.arange(n) / sps + env = np.exp(-((t - 1.5) ** 2) / (2 * 0.3 ** 2)) + x = 0.9 * env * np.sin(2 * np.pi * 27.0 * t) + f, v = channel_compliance_points(x, sps) + assert len(f) > 20 + assert v.max() >= 0.99 * np.abs(x).max() + assert 20.0 < f[int(np.argmax(v))] < 35.0 diff --git a/tests/test_sensor_check.py b/tests/test_sensor_check.py new file mode 100644 index 0000000..b39924a --- /dev/null +++ b/tests/test_sensor_check.py @@ -0,0 +1,67 @@ +"""Blastware sensor self-check waveform decode (minimateplus.sensor_check). + +Reverse-engineered 2026-09-15 against 7 BE12844 (MiniMate Plus) oracle events. +After the main waveform record-chain and the trailing metadata / per-channel +calibration records, a series-3 binary carries four length-prefixed records +tagged 0x3c-0x3f: the sensor self-check traces the unit records when it pulses +each sensor before monitoring (Blastware draws these as the little waveforms in +the "Sensor Check" strip on the right of the Event Report). + + * 0x3c / 0x3d / 0x3e = Tran / Vert / Long geophone ring-downs. + * 0x3f = MicL, a pulse train at the mic self-test frequency. + +The self-check injects a fixed pulse, so the response is near-identical across +events — asserted here as an invariant shape (damped one-sided ring-down for +the geophones, a multi-pulse train for the mic). +""" +from pathlib import Path + +import numpy as np + +from minimateplus.sensor_check import decode_sensor_check + +FIXDIR = Path(__file__).parent / "fixtures" / "fft-oracle-2026-09-14" +EVENTS = sorted(p.name for p in FIXDIR.iterdir()) # 7 BE12844 event binaries + + +def _decode(name): + return decode_sensor_check((FIXDIR / name).read_bytes()) + + +def test_all_four_channels_present(): + for name in EVENTS: + sc = _decode(name) + assert set(sc) == {"Tran", "Vert", "Long", "MicL"}, name + + +def test_geo_channels_are_damped_ringdowns(): + # Each geophone self-check is a large one-sided deflection (~-990 raw) that + # rings back and damps toward a settled value well above the trough. + for name in EVENTS: + sc = _decode(name) + for ch in ("Tran", "Vert", "Long"): + tr = np.asarray(sc[ch], dtype=float) + assert 240 <= len(tr) <= 260, f"{name}:{ch} n={len(tr)}" + assert abs(tr[:3].mean()) < 50, f"{name}:{ch} starts off-baseline" + assert tr.min() < -800, f"{name}:{ch} min {tr.min()}" + assert tr.max() < 60, f"{name}:{ch} unexpected positive swing {tr.max()}" + # damped: settles between the trough and zero, well above the trough + assert tr.min() < tr[-1] < 0, f"{name}:{ch} end {tr[-1]} not between trough and 0" + assert abs(tr[-1]) < 0.6 * abs(tr.min()), f"{name}:{ch} not damped, end {tr[-1]}" + + +def test_mic_channel_is_a_pulse_train(): + for name in EVENTS: + tr = np.asarray(_decode(name)["MicL"], dtype=float) + assert 235 <= len(tr) <= 255, f"{name} mic n={len(tr)}" + # larger dynamic range than the geo ring-down, and swings both ways + assert tr.min() < -1500, f"{name} mic min {tr.min()}" + assert tr.max() > 100, f"{name} mic max {tr.max()}" + # multiple pulses: several deep local minima + deep = (tr[1:-1] < tr[:-2]) & (tr[1:-1] < tr[2:]) & (tr[1:-1] < -800) + assert int(deep.sum()) >= 4, f"{name} mic pulses {int(deep.sum())}" + + +def test_returns_empty_when_no_sensor_check_block(): + assert decode_sensor_check(b"not a blastware file") == {} + assert decode_sensor_check(b"") == {} diff --git a/tests/test_waveform_fft.py b/tests/test_waveform_fft.py new file mode 100644 index 0000000..e1212fa --- /dev/null +++ b/tests/test_waveform_fft.py @@ -0,0 +1,85 @@ +"""Blastware-compatible channel FFT (waveform_fft). + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events, each with +a Blastware FFT report as ground truth. The recipe (DC-remove, no window, +zero-pad to 4096 → 0.25 Hz bins, single-sided 2/N amplitude) reproduces +Blastware's dominant frequency to the exact bin on all 28 channels and the +amplitude to report precision. +""" +from pathlib import Path + +import numpy as np + +from waveform_fft import channel_spectrum, dominant_frequency +from minimateplus.waveform_codec import decode_waveform_v2 + +FIXDIR = Path(__file__).parent / "fixtures" / "fft-oracle-2026-09-14" +GEO_LSB = 0.005 # 1 decode unit = 16 ADC counts = 0.005 in/s (series-3 Normal range) + +# Blastware FFT-report ground truth: file → {channel: (dominant_hz, amplitude_ips)}. +# amplitude is None where the channel is at the noise floor (report amp 0.000/0.001) +# — the dominant frequency still matches exactly, but the amplitude isn't meaningful. +ORACLE = { + "N844LPGH.VV0W": {"Tran": (27.00, 0.018), "Vert": (26.75, 0.009), "Long": (26.50, 0.021), "MicL": (2.000, None)}, + "N844LPPR.3S0W": {"Tran": (30.75, None), "Vert": (46.75, None), "Long": (26.75, None), "MicL": (49.50, None)}, + "N844LQHB.ZT0W": {"Tran": (19.75, 0.040), "Vert": (26.50, 0.018), "Long": (26.50, 0.083), "MicL": (2.750, None)}, + "N844LQUE.T50W": {"Tran": (21.50, 0.080), "Vert": (14.25, 0.028), "Long": (28.50, 0.046), "MicL": (5.750, None)}, + "N844LR8W.790W": {"Tran": (31.00, None), "Vert": (31.00, None), "Long": (34.00, None), "MicL": (66.25, None)}, + "N844LRCO.G60W": {"Tran": (32.25, 0.009), "Vert": (32.00, 0.005), "Long": (32.00, 0.008), "MicL": (32.00, None)}, + "N844LRCW.F30W": {"Tran": (21.25, 0.010), "Vert": (42.25, 0.002), "Long": (21.25, 0.014), "MicL": (21.25, None)}, +} + + +def test_pure_sine_frequency_and_amplitude(): + # A pure sine at a bin-centre frequency (128 cycles over 4096 samples) has no + # leakage, so the single-sided 2/N normalisation returns the amplitude exactly. + sps, n, f0, amp = 1024.0, 4096, 32.0, 0.5 + x = amp * np.sin(2 * np.pi * f0 * np.arange(n) / sps) + freqs, amps = channel_spectrum(x, sps=sps, nfft=4096) + fpk, apk = dominant_frequency(freqs, amps) + assert fpk == 32.0 + assert abs(apk - amp) < 1e-3 + + +def test_bin_resolution_is_quarter_hz(): + freqs, _ = channel_spectrum(np.zeros(3328), sps=1024.0, nfft=4096) + assert abs((freqs[1] - freqs[0]) - 0.25) < 1e-9 + + +def test_empty_input(): + freqs, amps = channel_spectrum([]) + assert len(freqs) == 0 and len(amps) == 0 + + +def _spectra(fname): + raw = (FIXDIR / fname).read_bytes() + dec = decode_waveform_v2(raw[raw.find(b"STRT") + 21:]) + out = {} + for ch, samples in dec.items(): + ips = np.asarray(samples, float) * GEO_LSB + out[ch] = channel_spectrum(ips, sps=1024.0) + return out + + +def test_dominant_frequency_matches_blastware_exactly(): + misses = [] + for fname, chans in ORACLE.items(): + spectra = _spectra(fname) + for ch, (want_hz, _) in chans.items(): + got_hz, _ = dominant_frequency(*spectra[ch]) + if abs(got_hz - want_hz) > 0.25: + misses.append(f"{fname}:{ch} got {got_hz} want {want_hz}") + assert not misses, "dominant-frequency mismatches:\n" + "\n".join(misses) + + +def test_amplitude_matches_blastware(): + misses = [] + for fname, chans in ORACLE.items(): + spectra = _spectra(fname) + for ch, (_, want_amp) in chans.items(): + if want_amp is None: + continue + _, got_amp = dominant_frequency(*spectra[ch]) + if abs(got_amp - want_amp) > 0.0015: + misses.append(f"{fname}:{ch} got {got_amp:.4f} want {want_amp:.3f}") + assert not misses, "amplitude mismatches:\n" + "\n".join(misses) diff --git a/waveform_fft.py b/waveform_fft.py new file mode 100644 index 0000000..7e84565 --- /dev/null +++ b/waveform_fft.py @@ -0,0 +1,66 @@ +"""Blastware-compatible FFT of a decoded seismograph channel. + +Pure numpy; no I/O, no device or DB dependencies. Feed it a channel's decoded +samples **in the unit you want the amplitudes in** (e.g. in/s) and it returns the +single-sided amplitude spectrum that Blastware's *FFT Report* draws. + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events with +Blastware FFT reports as ground truth. The recipe reproduces Blastware's +**dominant frequency to the exact 0.25 Hz bin on all 28 channels** and the +amplitude to report precision: + + 1. remove the DC component (subtract the mean); **no window** — a window + smears the peak and measurably worsens the match, + 2. zero-pad to ``nfft`` (4096 → 0.25 Hz bins at 1024 sps — Blastware's + resolution), + 3. single-sided amplitude ``A[k] = 2·|X[k]| / N`` where ``N`` is the real + sample count (not ``nfft``). + +The compliance chart (USBM RI8507 / OSMRE) is this spectrum's ``(freq, amp)`` +points plotted against the regulatory limit curve; the #10 FFT view is the +spectrum itself. +""" +from __future__ import annotations + +import numpy as np + +BW_NFFT = 4096 # 0.25 Hz bins at 1024 sps — Blastware's FFT resolution +BW_FMIN = 2.0 # dominant-frequency search floor (Hz) +BW_FMAX = 250.0 # dominant-frequency search ceiling (Hz) + + +def channel_spectrum(samples, sps: float = 1024.0, nfft: int = BW_NFFT): + """Single-sided amplitude spectrum of one channel, Blastware-compatible. + + ``samples`` is a 1-D sequence in the desired amplitude unit (in/s). Returns + ``(freqs, amps)`` numpy arrays covering ``0 .. sps/2`` in ``sps/nfft`` steps. + + Records longer than ``nfft`` are truncated by the transform — untested + against Blastware for that case (real MiniMate Plus records are ≤ ~3.3 s, + well under 4096 samples at 1024 sps). + """ + x = np.asarray(samples, dtype=float) + n = x.size + if n == 0: + return np.empty(0), np.empty(0) + x = x - x.mean() # DC removal, no window + mag = np.abs(np.fft.rfft(x, nfft)) + freqs = np.fft.rfftfreq(nfft, 1.0 / sps) + amps = (2.0 / n) * mag # single-sided amplitude + return freqs, amps + + +def dominant_frequency(freqs, amps, fmin: float = BW_FMIN, fmax: float = BW_FMAX): + """Peak ``(frequency_hz, amplitude)`` of a spectrum within ``[fmin, fmax)``. + + Matches Blastware's "Dominant Frequency" — the largest spectral bin in the + reportable band (below 2 Hz is baseline/DC drift, above 250 Hz is noise). + """ + freqs = np.asarray(freqs) + amps = np.asarray(amps) + lo = int(np.searchsorted(freqs, fmin)) + hi = int(np.searchsorted(freqs, fmax)) + if hi <= lo: + return 0.0, 0.0 + k = lo + int(np.argmax(amps[lo:hi])) + return float(freqs[k]), float(amps[k])