tests/test_waveform_fft.py

@@ -0,0 +1,85 @@
+"""Blastware-compatible channel FFT (waveform_fft).
+
+Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events, each with
+a Blastware FFT report as ground truth. The recipe (DC-remove, no window,
+zero-pad to 4096 → 0.25 Hz bins, single-sided 2/N amplitude) reproduces
+Blastware's dominant frequency to the exact bin on all 28 channels and the
+amplitude to report precision.
+"""
+from pathlib import Path
+
+import numpy as np
+
+from waveform_fft import channel_spectrum, dominant_frequency
+from minimateplus.waveform_codec import decode_waveform_v2
+
+FIXDIR = Path(__file__).parent / "fixtures" / "fft-oracle-2026-09-14"
+GEO_LSB = 0.005   # 1 decode unit = 16 ADC counts = 0.005 in/s (series-3 Normal range)
+
+# Blastware FFT-report ground truth: file → {channel: (dominant_hz, amplitude_ips)}.
+# amplitude is None where the channel is at the noise floor (report amp 0.000/0.001)
+# — the dominant frequency still matches exactly, but the amplitude isn't meaningful.
+ORACLE = {
+    "N844LPGH.VV0W": {"Tran": (27.00, 0.018), "Vert": (26.75, 0.009), "Long": (26.50, 0.021), "MicL": (2.000, None)},
+    "N844LPPR.3S0W": {"Tran": (30.75, None),  "Vert": (46.75, None),  "Long": (26.75, None),  "MicL": (49.50, None)},
+    "N844LQHB.ZT0W": {"Tran": (19.75, 0.040), "Vert": (26.50, 0.018), "Long": (26.50, 0.083), "MicL": (2.750, None)},
+    "N844LQUE.T50W": {"Tran": (21.50, 0.080), "Vert": (14.25, 0.028), "Long": (28.50, 0.046), "MicL": (5.750, None)},
+    "N844LR8W.790W": {"Tran": (31.00, None),  "Vert": (31.00, None),  "Long": (34.00, None),  "MicL": (66.25, None)},
+    "N844LRCO.G60W": {"Tran": (32.25, 0.009), "Vert": (32.00, 0.005), "Long": (32.00, 0.008), "MicL": (32.00, None)},
+    "N844LRCW.F30W": {"Tran": (21.25, 0.010), "Vert": (42.25, 0.002), "Long": (21.25, 0.014), "MicL": (21.25, None)},
+}
+
+
+def test_pure_sine_frequency_and_amplitude():
+    # A pure sine at a bin-centre frequency (128 cycles over 4096 samples) has no
+    # leakage, so the single-sided 2/N normalisation returns the amplitude exactly.
+    sps, n, f0, amp = 1024.0, 4096, 32.0, 0.5
+    x = amp * np.sin(2 * np.pi * f0 * np.arange(n) / sps)
+    freqs, amps = channel_spectrum(x, sps=sps, nfft=4096)
+    fpk, apk = dominant_frequency(freqs, amps)
+    assert fpk == 32.0
+    assert abs(apk - amp) < 1e-3
+
+
+def test_bin_resolution_is_quarter_hz():
+    freqs, _ = channel_spectrum(np.zeros(3328), sps=1024.0, nfft=4096)
+    assert abs((freqs[1] - freqs[0]) - 0.25) < 1e-9
+
+
+def test_empty_input():
+    freqs, amps = channel_spectrum([])
+    assert len(freqs) == 0 and len(amps) == 0
+
+
+def _spectra(fname):
+    raw = (FIXDIR / fname).read_bytes()
+    dec = decode_waveform_v2(raw[raw.find(b"STRT") + 21:])
+    out = {}
+    for ch, samples in dec.items():
+        ips = np.asarray(samples, float) * GEO_LSB
+        out[ch] = channel_spectrum(ips, sps=1024.0)
+    return out
+
+
+def test_dominant_frequency_matches_blastware_exactly():
+    misses = []
+    for fname, chans in ORACLE.items():
+        spectra = _spectra(fname)
+        for ch, (want_hz, _) in chans.items():
+            got_hz, _ = dominant_frequency(*spectra[ch])
+            if abs(got_hz - want_hz) > 0.25:
+                misses.append(f"{fname}:{ch} got {got_hz} want {want_hz}")
+    assert not misses, "dominant-frequency mismatches:\n" + "\n".join(misses)
+
+
+def test_amplitude_matches_blastware():
+    misses = []
+    for fname, chans in ORACLE.items():
+        spectra = _spectra(fname)
+        for ch, (_, want_amp) in chans.items():
+            if want_amp is None:
+                continue
+            _, got_amp = dominant_frequency(*spectra[ch])
+            if abs(got_amp - want_amp) > 0.0015:
+                misses.append(f"{fname}:{ch} got {got_amp:.4f} want {want_amp:.3f}")
+    assert not misses, "amplitude mismatches:\n" + "\n".join(misses)

waveform_fft.py

@@ -0,0 +1,66 @@
+"""Blastware-compatible FFT of a decoded seismograph channel.
+
+Pure numpy; no I/O, no device or DB dependencies. Feed it a channel's decoded
+samples **in the unit you want the amplitudes in** (e.g. in/s) and it returns the
+single-sided amplitude spectrum that Blastware's *FFT Report* draws.
+
+Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events with
+Blastware FFT reports as ground truth. The recipe reproduces Blastware's
+**dominant frequency to the exact 0.25 Hz bin on all 28 channels** and the
+amplitude to report precision:
+
+  1. remove the DC component (subtract the mean); **no window** — a window
+     smears the peak and measurably worsens the match,
+  2. zero-pad to ``nfft`` (4096 → 0.25 Hz bins at 1024 sps — Blastware's
+     resolution),
+  3. single-sided amplitude ``A[k] = 2·|X[k]| / N`` where ``N`` is the real
+     sample count (not ``nfft``).
+
+The compliance chart (USBM RI8507 / OSMRE) is this spectrum's ``(freq, amp)``
+points plotted against the regulatory limit curve; the #10 FFT view is the
+spectrum itself.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+BW_NFFT = 4096              # 0.25 Hz bins at 1024 sps — Blastware's FFT resolution
+BW_FMIN = 2.0              # dominant-frequency search floor (Hz)
+BW_FMAX = 250.0           # dominant-frequency search ceiling (Hz)
+
+
+def channel_spectrum(samples, sps: float = 1024.0, nfft: int = BW_NFFT):
+    """Single-sided amplitude spectrum of one channel, Blastware-compatible.
+
+    ``samples`` is a 1-D sequence in the desired amplitude unit (in/s). Returns
+    ``(freqs, amps)`` numpy arrays covering ``0 .. sps/2`` in ``sps/nfft`` steps.
+
+    Records longer than ``nfft`` are truncated by the transform — untested
+    against Blastware for that case (real MiniMate Plus records are ≤ ~3.3 s,
+    well under 4096 samples at 1024 sps).
+    """
+    x = np.asarray(samples, dtype=float)
+    n = x.size
+    if n == 0:
+        return np.empty(0), np.empty(0)
+    x = x - x.mean()                          # DC removal, no window
+    mag = np.abs(np.fft.rfft(x, nfft))
+    freqs = np.fft.rfftfreq(nfft, 1.0 / sps)
+    amps = (2.0 / n) * mag                     # single-sided amplitude
+    return freqs, amps
+
+
+def dominant_frequency(freqs, amps, fmin: float = BW_FMIN, fmax: float = BW_FMAX):
+    """Peak ``(frequency_hz, amplitude)`` of a spectrum within ``[fmin, fmax)``.
+
+    Matches Blastware's "Dominant Frequency" — the largest spectral bin in the
+    reportable band (below 2 Hz is baseline/DC drift, above 250 Hz is noise).
+    """
+    freqs = np.asarray(freqs)
+    amps = np.asarray(amps)
+    lo = int(np.searchsorted(freqs, fmin))
+    hi = int(np.searchsorted(freqs, fmax))
+    if hi <= lo:
+        return 0.0, 0.0
+    k = lo + int(np.argmax(amps[lo:hi]))
+    return float(freqs[k]), float(amps[k])