Compare commits
17
Commits
88c0e2b765
...
8265e32ad5
@@ -4,6 +4,84 @@ 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 (series-3 *and* series-4).** Analyses Blastware/Thor
|
||||
derive from event data, reverse-engineered against BE12844 (MiniMate Plus) and
|
||||
UM (Thor) events and reproduced in seismo-relay: the compliance chart and the
|
||||
sensor-check strip render on the event-report PDF, and the FFT reproduces
|
||||
Blastware's FFT Report.
|
||||
|
||||
The FFT and the compliance chart are additive and read from the existing `.h5`
|
||||
samples — no `.h5` or DB change for those. The **sensor self-check** now lives
|
||||
in the standardized `.h5` (schema **v2**, a new `/sensor_check` group) so SFM
|
||||
serves it device-agnostically rather than decoding at report time.
|
||||
|
||||
⚠ **The sensor-check needs a backfill.** Existing `.h5` files are schema v1 and
|
||||
carry no `/sensor_check` group, so their reports show no sensor-check strip
|
||||
until regenerated. `TOOL_VERSION` is bumped to **0.31.0**, so the standard
|
||||
backfill regenerates every event and picks up the traces with **no `--force`**:
|
||||
`scripts/backfill_thor_events.py` for series-4 (it already owed a v0.30.0 Thor
|
||||
backfill — this rides along) and the series-3 sidecar/shape backfill for
|
||||
MiniMate events. Purely additive — no decoded value changes, and v1 `.h5` files
|
||||
read fine until then (empty strip). DB backup first, as always.
|
||||
|
||||
### 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 — both series.** The little
|
||||
"Sensor Check" traces (geophone ring-downs — the transducer's damped impulse
|
||||
response — plus a MicL pulse train, the mic's known-signal gain check) are the
|
||||
unit's proof its sensors were healthy when it recorded the event.
|
||||
- **Series-3** (`minimateplus.sensor_check`): four records (`0x3c`–`0x3f`) in
|
||||
the binary's trailing block, same delta-block codec as the main waveform.
|
||||
Verified against all 7 BE12844 reports (mic zero-crossing = 20.1 Hz exact;
|
||||
geophone ring-downs ~7.5 Hz, overswing ~3.5).
|
||||
- **Series-4** (`micromate.sensor_check`): the same self-test in the Thor IDFW
|
||||
fixed header — four `01 0e 3c/3d/3e/3f` records (same channel ids) storing
|
||||
raw int16 traces; three-channel (mic-disabled) units carry only the three
|
||||
geophones. Validated by shape + cross-event consistency.
|
||||
- **Standardized into the `.h5`** (`/sensor_check`, schema v2): each series'
|
||||
decoder attaches the traces to the event at decode, the writer persists
|
||||
them, and `gather_report_data` reads them back — so SFM renders the strip
|
||||
(flush against the waveform panel) plus the **Sensor Check → Frequency /
|
||||
Overswing Ratio** sub-rows without knowing the source instrument.
|
||||
- Tests: `tests/test_sensor_check.py`, `tests/test_sensor_check_idf.py`,
|
||||
`tests/test_event_hdf5_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
|
||||
|
||||
@@ -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 —
|
||||
<https://www.ecfr.gov/current/title-30/chapter-VII/subchapter-K/part-816/section-816.67>
|
||||
@@ -0,0 +1,89 @@
|
||||
r"""Decode the Thor / Micromate (series-4) sensor self-check waveforms from an
|
||||
IDFW event binary.
|
||||
|
||||
Reverse-engineered 2026-09-15 against 4 UM (Thor) oracle events. The IDFW
|
||||
binary carries the sensor self-check in its fixed-header region (before the
|
||||
waveform body), as up to four records tagged ``01 0e 3c/3d/3e/3f`` — the SAME
|
||||
channel ids as the series-3 MiniMate Plus (Tran / Vert / Long / MicL), which is
|
||||
the physical self-test:
|
||||
|
||||
* 3c / 3d / 3e = Tran / Vert / Long geophone ring-downs (a damped impulse
|
||||
response — resonant frequency + damping).
|
||||
* 3f = MicL pulse train (the mic's known-signal gain check). Absent
|
||||
on three-channel (mic-disabled) units.
|
||||
|
||||
Record framing (per record)::
|
||||
|
||||
01 0e [id:1] [flags:3] [count:2 BE] [pad:10] [int16-BE samples × count]
|
||||
\___ 18-byte header ___/
|
||||
|
||||
Unlike series-3's delta-coded trailing block, series-4 stores each trace as a
|
||||
raw int16 big-endian array. ``count`` (the 2-byte field at header offset +8)
|
||||
is the sample count; the record is padded to a fixed stride after that.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import struct
|
||||
from typing import Dict, List
|
||||
|
||||
# Record id → channel. Same ids/order as series-3 (minimateplus.sensor_check).
|
||||
_ID_TO_CHANNEL = {0x3C: "Tran", 0x3D: "Vert", 0x3E: "Long", 0x3F: "MicL"}
|
||||
_CHAIN_IDS = (0x3C, 0x3D, 0x3E, 0x3F)
|
||||
|
||||
_MARKER = b"\x01\x0e" # precedes the 1-byte channel id
|
||||
_HEADER_LEN = 18 # bytes from the marker start to the first sample
|
||||
_COUNT_OFF = 8 # 2-byte BE sample count, from the marker start
|
||||
_MAX_COUNT = 4000 # sanity cap (traces are ~70-200 samples)
|
||||
|
||||
|
||||
def _find_chain(raw: bytes):
|
||||
"""Locate the sensor-check record chain. Returns a list of
|
||||
``(offset, id, count)`` for the first run of markers whose ids run
|
||||
3c, 3d, 3e[, 3f] in order, or ``[]``.
|
||||
|
||||
Records are padded to a fixed stride, so the next marker is not at
|
||||
``header + count*2``; instead collect every ``01 0e [id]`` marker with a
|
||||
sane count and take the first id-ordered run. Validating the id sequence
|
||||
(not a lone ``01 0e 3c``) keeps a stray marker in the waveform body from
|
||||
matching — the real chain sits in the fixed header, ahead of the body.
|
||||
"""
|
||||
n = len(raw)
|
||||
markers = []
|
||||
for p in range(n - _HEADER_LEN):
|
||||
if raw[p:p + 2] == _MARKER and raw[p + 2] in _ID_TO_CHANNEL:
|
||||
count = int.from_bytes(raw[p + _COUNT_OFF:p + _COUNT_OFF + 2], "big")
|
||||
if 0 < count <= _MAX_COUNT:
|
||||
markers.append((p, raw[p + 2], count))
|
||||
|
||||
for i, (off, rid, _c) in enumerate(markers):
|
||||
if rid != 0x3C:
|
||||
continue
|
||||
run = [markers[i]]
|
||||
for m in markers[i + 1:]:
|
||||
if len(run) < len(_CHAIN_IDS) and m[1] == _CHAIN_IDS[len(run)]:
|
||||
run.append(m)
|
||||
else:
|
||||
break
|
||||
if len(run) >= 3: # 3-channel (mic-disabled) units are valid
|
||||
return run
|
||||
return []
|
||||
|
||||
|
||||
def decode_idf_sensor_check(raw: bytes) -> Dict[str, List[int]]:
|
||||
"""Decode the sensor self-check traces from a Thor/Micromate IDFW binary.
|
||||
|
||||
Returns ``{"Tran": [...], "Vert": [...], "Long": [...], "MicL": [...]}`` in
|
||||
raw int16 ADC counts (MicL omitted on 3-channel units), or ``{}`` if the
|
||||
binary carries no sensor-check chain (a non-IDF file, or an IDFH histogram).
|
||||
"""
|
||||
chain = _find_chain(raw)
|
||||
if not chain:
|
||||
return {}
|
||||
out: Dict[str, List[int]] = {}
|
||||
for off, rid, count in chain:
|
||||
start = off + _HEADER_LEN
|
||||
blob = raw[start:start + count * 2]
|
||||
if len(blob) < count * 2:
|
||||
continue
|
||||
out[_ID_TO_CHANNEL[rid]] = list(struct.unpack(">%dh" % count, blob))
|
||||
return out
|
||||
@@ -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)
|
||||
@@ -50,7 +50,7 @@ SIDECAR_KIND = "sfm.event"
|
||||
# bumped without a `pip install` re-run — leading to confusing stale
|
||||
# version stamps in sidecars. Bump this constant and CHANGELOG.md
|
||||
# together at release time.
|
||||
TOOL_VERSION = "0.30.0"
|
||||
TOOL_VERSION = "0.31.0" # +/sensor_check group (schema v2); gates the backfill regen
|
||||
|
||||
try:
|
||||
# Best-effort: prefer the installed metadata when it's NEWER than the
|
||||
@@ -960,6 +960,11 @@ def read_blastware_file(path: Union[str, Path]) -> Event:
|
||||
project=project, client=client, operator=user, sensor_location=seisloc,
|
||||
)
|
||||
ev.raw_samples = samples
|
||||
# Sensor self-check traces from the binary's trailing block (waveform
|
||||
# events only; returns {} for histograms / when absent). Carried on the
|
||||
# Event so the .h5 writer persists them device-agnostically.
|
||||
from minimateplus.sensor_check import decode_sensor_check
|
||||
ev.sensor_check = decode_sensor_check(raw) or None
|
||||
# Only compute peaks from samples when we actually have samples.
|
||||
# For events the codec couldn't decode (histogram-mode bodies, until
|
||||
# the §7.6.2 histogram codec is wired in), samples is an empty dict
|
||||
|
||||
@@ -544,6 +544,15 @@ class Event:
|
||||
pretrig_samples: Optional[int] = None # from STRT record: pre-trigger sample count
|
||||
rectime_seconds: Optional[int] = None # from STRT record: record duration (seconds)
|
||||
|
||||
# Sensor self-check traces keyed by channel label — the short diagnostic
|
||||
# waveforms the unit records when it pulses each sensor before monitoring
|
||||
# (geophone ring-downs + a mic pulse train). Decoded from the binary by
|
||||
# the per-series decoder (minimateplus.sensor_check / micromate.sensor_check)
|
||||
# and carried here so the .h5 writer can persist them device-agnostically.
|
||||
# Raw ADC counts; the source series' scale differs but the trace is a
|
||||
# shape diagnostic (rendered fit-to-box). None when absent.
|
||||
sensor_check: Optional[dict] = None # {"Tran": [...], ..., "MicL": [...]}
|
||||
|
||||
# ── Debug / introspection ─────────────────────────────────────────────────
|
||||
# Raw 210-byte waveform record bytes, set when debug mode is active.
|
||||
# Exposed by the SFM server via ?debug=true so field layouts can be verified.
|
||||
|
||||
@@ -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
|
||||
@@ -305,6 +305,11 @@ def main(argv=None) -> int:
|
||||
default=0,
|
||||
)
|
||||
ev.total_samples = ev.total_samples or n_samp
|
||||
# Sensor self-check traces from the IDFW fixed
|
||||
# header, so regenerated .h5 files gain the v2
|
||||
# /sensor_check group (mirrors save_imported_idf).
|
||||
from micromate.sensor_check import decode_idf_sensor_check
|
||||
ev.sensor_check = decode_idf_sensor_check(binary_bytes) or None
|
||||
|
||||
event_hdf5.write_event_hdf5(
|
||||
hdf5_path, ev,
|
||||
|
||||
@@ -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)
|
||||
|
||||
|
||||
@@ -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")
|
||||
+44
-4
@@ -12,8 +12,11 @@ Layout written to `<filename>.h5`:
|
||||
├─ samples_int16/ (optional)
|
||||
│ ├─ Tran (int16, raw ADC counts) shape: (N,)
|
||||
│ └─ ... per channel (only when present in the source)
|
||||
├─ sensor_check/ (optional, schema v2+)
|
||||
│ ├─ Tran (int32, raw counts) shape: (M,) M ≪ N
|
||||
│ └─ ... per channel present in the source (MicL absent on 3-channel units)
|
||||
└─ root attrs (event metadata):
|
||||
schema_version int = 1
|
||||
schema_version int = 2
|
||||
kind str = "sfm.event.hdf5"
|
||||
serial str
|
||||
waveform_key str (8-hex)
|
||||
@@ -64,7 +67,7 @@ from minimateplus.models import Event
|
||||
|
||||
log = logging.getLogger(__name__)
|
||||
|
||||
SCHEMA_VERSION = 1
|
||||
SCHEMA_VERSION = 2 # v2 adds the optional /sensor_check group
|
||||
HDF5_KIND = "sfm.event.hdf5"
|
||||
|
||||
# Geophone full-scale velocity per range (in/s). Confirmed in CLAUDE.md
|
||||
@@ -270,6 +273,22 @@ def write_event_hdf5(
|
||||
)
|
||||
igrp.attrs["mic_psi_per_count"] = float(mic_factor)
|
||||
|
||||
# /sensor_check — optional short diagnostic self-check traces (schema
|
||||
# v2+). Raw ADC counts (a shape diagnostic; the per-series count scale
|
||||
# differs, and the renderer fits each trace to its box). Only channels
|
||||
# the decoder found are written — 3-channel units carry no MicL.
|
||||
sc = event.sensor_check or {}
|
||||
if sc:
|
||||
scgrp = f.create_group("sensor_check")
|
||||
for ch in ("Tran", "Vert", "Long", "MicL"):
|
||||
vals = sc.get(ch)
|
||||
if vals:
|
||||
scgrp.create_dataset(
|
||||
ch, data=np.asarray(vals, dtype=np.int32),
|
||||
compression="gzip", compression_opts=4, shuffle=True,
|
||||
)
|
||||
scgrp.attrs["units"] = "raw_counts"
|
||||
|
||||
import os
|
||||
os.replace(tmp, path)
|
||||
|
||||
@@ -334,6 +353,16 @@ def read_event_hdf5(path: Union[str, Path]) -> dict:
|
||||
if mic_attr is not None:
|
||||
mic_psi = float(mic_attr)
|
||||
|
||||
# /sensor_check — optional (schema v2+); absent on older files.
|
||||
sensor_check = None
|
||||
scgrp = f.get("sensor_check")
|
||||
if scgrp is not None:
|
||||
sensor_check = {}
|
||||
for ch in ("Tran", "Vert", "Long", "MicL"):
|
||||
ds = scgrp.get(ch)
|
||||
if ds is not None:
|
||||
sensor_check[ch] = np.asarray(ds[()])
|
||||
|
||||
return {
|
||||
"schema_version": sv,
|
||||
"kind": attrs.get("kind"),
|
||||
@@ -341,6 +370,7 @@ def read_event_hdf5(path: Union[str, Path]) -> dict:
|
||||
"samples": samples,
|
||||
"samples_int16": samples_int16,
|
||||
"mic_psi_per_count": mic_psi,
|
||||
"sensor_check": sensor_check,
|
||||
}
|
||||
|
||||
|
||||
@@ -431,11 +461,16 @@ def plot_json_from_hdf5(
|
||||
event_id: Optional[str] = None,
|
||||
index: Optional[int] = None,
|
||||
) -> dict:
|
||||
"""Build a `sfm.plot.v1` JSON dict from a stored .h5 file."""
|
||||
"""Build a `sfm.plot.v1` JSON dict from a stored .h5 file.
|
||||
|
||||
The dict also carries a top-level ``sensor_check`` key (the raw self-check
|
||||
traces as ``{ch: [int]}``, or None) beyond the plot schema, so report
|
||||
generation can read the traces from the same single .h5 load.
|
||||
"""
|
||||
data = read_event_hdf5(path)
|
||||
a = data["attrs"]
|
||||
s = data["samples"]
|
||||
return _build_plot_dict(
|
||||
out = _build_plot_dict(
|
||||
n_samples=len(s["Tran"]) if "Tran" in s else 0,
|
||||
sample_rate=int(a.get("sample_rate", 1024) or 1024),
|
||||
pretrig_samples=int(a.get("pretrig_samples", 0) or 0),
|
||||
@@ -463,6 +498,11 @@ def plot_json_from_hdf5(
|
||||
event_id=event_id,
|
||||
index=index,
|
||||
)
|
||||
scd = data.get("sensor_check")
|
||||
out["sensor_check"] = (
|
||||
{ch: v.tolist() for ch, v in scd.items()} if scd else None
|
||||
)
|
||||
return out
|
||||
|
||||
|
||||
def _build_plot_dict(
|
||||
|
||||
+162
-35
@@ -121,6 +121,13 @@ class ReportData:
|
||||
t0_ms: Optional[float] = None
|
||||
dt_ms: Optional[float] = None
|
||||
|
||||
# Sensor self-check traces — {ch: [samples]} in raw counts, read from the
|
||||
# standardized .h5 (/sensor_check group, schema v2+) where the per-series
|
||||
# decoder stored them at ingest. The little diagnostic waveforms BW draws
|
||||
# in its "Sensor Check" strip. Empty when absent (pre-v2 .h5, histogram,
|
||||
# or 3-channel unit's MicL).
|
||||
sensor_check_waveforms: dict = field(default_factory=dict)
|
||||
|
||||
# Record-type discriminator
|
||||
record_type: Optional[str] = None
|
||||
is_histogram: bool = False
|
||||
@@ -246,6 +253,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,
|
||||
})
|
||||
@@ -287,6 +296,12 @@ def gather_report_data(
|
||||
rd.pretrig_samples = ta.get("pretrig_samples")
|
||||
rd.t0_ms = ta.get("t0_ms")
|
||||
rd.dt_ms = ta.get("dt_ms")
|
||||
# Sensor self-check traces — read from the standardized .h5 (schema
|
||||
# v2+). Device-agnostic: whichever decoder produced the event
|
||||
# stored them at ingest, so SFM reads them here without knowing or
|
||||
# caring about the source instrument series. Empty on pre-v2 files
|
||||
# (until backfilled) and on 3-channel / histogram events.
|
||||
rd.sensor_check_waveforms = wf.get("sensor_check") or {}
|
||||
except Exception as exc:
|
||||
log.warning("gather_report_data: hdf5 read failed: %s", exc)
|
||||
|
||||
@@ -396,9 +411,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 +517,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 +584,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 +617,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 +674,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 +746,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 +794,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 +820,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 +840,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 +877,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 +900,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 +920,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
|
||||
|
||||
@@ -662,6 +662,11 @@ class WaveformStore:
|
||||
ev.raw_samples = idf_samples
|
||||
n_samples = max((len(idf_samples.get(ch, [])) for ch in ("Tran", "Vert", "Long", "MicL")), default=0)
|
||||
ev.total_samples = ev.total_samples or n_samples
|
||||
# Sensor self-check traces from the IDFW fixed header (waveform
|
||||
# events only; {} on histograms / when absent). Carried on the
|
||||
# bridged Event so the .h5 writer persists them like series-3.
|
||||
from micromate.sensor_check import decode_idf_sensor_check
|
||||
ev.sensor_check = decode_idf_sensor_check(idf_bytes) or None
|
||||
|
||||
# For IDFH histograms there are no per-sample waveform arrays — the
|
||||
# device stores one peak ADC count per interval per channel. Synthesise
|
||||
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -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
|
||||
@@ -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
|
||||
@@ -0,0 +1,71 @@
|
||||
"""The event .h5 carries the sensor self-check traces (schema v2).
|
||||
|
||||
The sensor check is decoded by the per-series decoder and attached to the
|
||||
standardized Event, so the .h5 writer persists it device-agnostically and SFM
|
||||
reads it back without knowing which instrument produced it. Old v1 files (no
|
||||
sensor_check group) must still read cleanly.
|
||||
"""
|
||||
import tempfile
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
|
||||
from minimateplus.models import Event
|
||||
from minimateplus.event_file_io import read_blastware_file
|
||||
from sfm import event_hdf5
|
||||
|
||||
S3_FIX = Path(__file__).parent / "fixtures" / "fft-oracle-2026-09-14" / "N844LQHB.ZT0W"
|
||||
|
||||
|
||||
def _write(ev, **kw):
|
||||
d = Path(tempfile.mkdtemp())
|
||||
p = d / "e.h5"
|
||||
event_hdf5.write_event_hdf5(p, ev, serial="BE12844", **kw)
|
||||
return p
|
||||
|
||||
|
||||
def test_sensor_check_roundtrips_through_hdf5():
|
||||
ev = Event(index=0)
|
||||
ev.raw_samples = {"Tran": [1, 2, -3], "Vert": [0, 1], "Long": [2], "MicL": [5, -5]}
|
||||
ev.sample_rate = 1024
|
||||
sc = {"Tran": [0, -990, -500, -100], "Vert": [0, -980, -480],
|
||||
"Long": [0, -986, -470], "MicL": [0, -1800, 1800, -1800]}
|
||||
ev.sensor_check = sc
|
||||
|
||||
r = event_hdf5.read_event_hdf5(_write(ev))
|
||||
assert r["schema_version"] == 2
|
||||
assert set(r["sensor_check"]) == {"Tran", "Vert", "Long", "MicL"}
|
||||
for ch, vals in sc.items():
|
||||
assert r["sensor_check"][ch].tolist() == vals
|
||||
|
||||
|
||||
def test_plot_json_carries_sensor_check():
|
||||
ev = Event(index=0)
|
||||
ev.raw_samples = {"Tran": [1, 2, 3]}
|
||||
ev.sample_rate = 1024
|
||||
ev.sensor_check = {"Tran": [0, -990, -500], "Vert": [0, -980],
|
||||
"Long": [0, -986]} # 3-channel: no MicL
|
||||
pj = event_hdf5.plot_json_from_hdf5(_write(ev))
|
||||
assert pj["sensor_check"] is not None
|
||||
assert "MicL" not in pj["sensor_check"]
|
||||
assert pj["sensor_check"]["Tran"] == [0, -990, -500]
|
||||
|
||||
|
||||
def test_event_without_sensor_check_still_reads_as_v2():
|
||||
ev = Event(index=0)
|
||||
ev.raw_samples = {"Tran": [1, 2, 3]}
|
||||
ev.sample_rate = 1024
|
||||
r = event_hdf5.read_event_hdf5(_write(ev))
|
||||
assert r["schema_version"] == 2
|
||||
assert r["sensor_check"] is None
|
||||
assert event_hdf5.plot_json_from_hdf5(_write(ev))["sensor_check"] is None
|
||||
|
||||
|
||||
def test_series3_decode_populates_event_sensor_check():
|
||||
# The real series-3 decoder attaches the traces to the Event, so the
|
||||
# ingest/backfill .h5 write picks them up with no extra plumbing.
|
||||
ev = read_blastware_file(S3_FIX)
|
||||
assert ev.sensor_check is not None
|
||||
assert set(ev.sensor_check) == {"Tran", "Vert", "Long", "MicL"}
|
||||
tran = np.asarray(ev.sensor_check["Tran"], dtype=float)
|
||||
assert tran.min() < -800 # the geophone ring-down deflection
|
||||
@@ -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"") == {}
|
||||
@@ -0,0 +1,67 @@
|
||||
"""Series-4 (Thor / Micromate IDFW) sensor self-check waveform decode.
|
||||
|
||||
Reverse-engineered 2026-09-15 against 4 UM (Thor) oracle events. The IDFW
|
||||
binary carries the sensor self-check in its fixed-header region (before the
|
||||
waveform body) as up to four records tagged ``01 0e 3c/3d/3e/3f`` — the SAME
|
||||
channel ids as series-3 (Tran/Vert/Long/MicL). Unlike series-3's delta-coded
|
||||
trailing block, series-4 stores each trace as a raw int16-BE array after an
|
||||
18-byte record header whose sample count is a 2-byte field at offset +8.
|
||||
|
||||
Three-channel (mic-disabled) Thor units carry only 3c/3d/3e — no MicL record.
|
||||
|
||||
Validated by shape (geophone ring-down / mic pulse train) and cross-event
|
||||
consistency, since there's no Thor Event-Report strip to exact-match against.
|
||||
"""
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
|
||||
from micromate.sensor_check import decode_idf_sensor_check
|
||||
|
||||
FIXDIR = Path(__file__).parent / "fixtures" / "thor-idf-sc"
|
||||
EVENTS = sorted(p.name for p in FIXDIR.glob("*.IDFW"))
|
||||
|
||||
|
||||
def _decode(name):
|
||||
return decode_idf_sensor_check((FIXDIR / name).read_bytes())
|
||||
|
||||
|
||||
def test_geo_channels_present_and_ringdown_shaped():
|
||||
# Every IDFW event has the three geophone self-checks; each is a large
|
||||
# one-sided deflection (~15000 raw counts) that rings back — the geophone's
|
||||
# damped impulse response.
|
||||
for name in EVENTS:
|
||||
sc = _decode(name)
|
||||
for ch in ("Tran", "Vert", "Long"):
|
||||
assert ch in sc, f"{name} missing {ch}"
|
||||
tr = np.asarray(sc[ch], dtype=float)
|
||||
tr = tr - tr[:4].mean() # reference to the pre-trigger baseline
|
||||
assert 40 <= len(tr) <= 300, f"{name}:{ch} n={len(tr)}"
|
||||
assert tr.min() < -8000, f"{name}:{ch} min {tr.min()}"
|
||||
# deflects one way and rings back toward / past the baseline
|
||||
assert tr.max() < abs(tr.min()), f"{name}:{ch} not one-sided"
|
||||
|
||||
|
||||
def test_mic_present_only_on_four_channel_units():
|
||||
# UM11719 / UM12947 record a mic; UM13981 / UM20147 are 3-channel
|
||||
# (mic-disabled) units and carry no MicL self-check.
|
||||
got = {name: ("MicL" in _decode(name)) for name in EVENTS}
|
||||
assert any(got.values()), "expected at least one 4-channel unit"
|
||||
assert not all(got.values()), "expected at least one 3-channel unit"
|
||||
for name, has_mic in got.items():
|
||||
if has_mic:
|
||||
tr = np.asarray(_decode(name)["MicL"], dtype=float)
|
||||
tr = tr - tr[:4].mean()
|
||||
# mic self-check is a bipolar pulse train — swings both ways, wide range
|
||||
assert tr.max() > 5000 and tr.min() < -5000, f"{name} mic not bipolar"
|
||||
|
||||
|
||||
def test_channel_ids_and_order():
|
||||
# ids decode to the canonical channel names, geo always in Tran/Vert/Long order
|
||||
sc = _decode(EVENTS[0])
|
||||
assert [c for c in ("Tran", "Vert", "Long") if c in sc] == ["Tran", "Vert", "Long"]
|
||||
|
||||
|
||||
def test_returns_empty_on_non_idf_input():
|
||||
assert decode_idf_sensor_check(b"not an IDF file") == {}
|
||||
assert decode_idf_sensor_check(b"") == {}
|
||||
@@ -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)
|
||||
@@ -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])
|
||||
Reference in New Issue
Block a user