Merge feat/fft-series3: Blastware FFT, USBM compliance chart, sensor self-check
Reverse-engineered Blastware Event/FFT-Report parity, all additive (reads the existing .h5 samples + retained raw binary, no DB/.h5 change or backfill): - Blastware-compatible channel FFT (waveform_fft) - USBM RI8507/OSMRE compliance chart on the event-report PDF (sfm/compliance) - sensor self-check strip decoded from the series-3 binary trailing block (minimateplus/sensor_check) + Frequency/Overswing sub-rows - seismo_lab Inspector hex reader (minimateplus/binary_annotate) - report fixes: stacked-lane y-tick collision, header serial fit Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01YDXjZCr4RqT2U3QvMDhgzf
This commit is contained in:
Claude Opus 4.8
@@ -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
|
||||
|
||||
@@ -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,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)
|
||||
@@ -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
|
||||
@@ -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")
|
||||
+167
-35
@@ -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
|
||||
|
||||
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,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,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