diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W b/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W new file mode 100644 index 0000000..785721f Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W b/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W new file mode 100644 index 0000000..b4395a1 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W b/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W new file mode 100644 index 0000000..6c7afb6 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W b/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W new file mode 100644 index 0000000..fd9a7cd Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W b/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W new file mode 100644 index 0000000..8726fdb Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W b/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W new file mode 100644 index 0000000..9e4fa8f Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W b/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W new file mode 100644 index 0000000..6e9ef75 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W differ diff --git a/tests/test_waveform_fft.py b/tests/test_waveform_fft.py new file mode 100644 index 0000000..e1212fa --- /dev/null +++ b/tests/test_waveform_fft.py @@ -0,0 +1,85 @@ +"""Blastware-compatible channel FFT (waveform_fft). + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events, each with +a Blastware FFT report as ground truth. The recipe (DC-remove, no window, +zero-pad to 4096 → 0.25 Hz bins, single-sided 2/N amplitude) reproduces +Blastware's dominant frequency to the exact bin on all 28 channels and the +amplitude to report precision. +""" +from pathlib import Path + +import numpy as np + +from waveform_fft import channel_spectrum, dominant_frequency +from minimateplus.waveform_codec import decode_waveform_v2 + +FIXDIR = Path(__file__).parent / "fixtures" / "fft-oracle-2026-09-14" +GEO_LSB = 0.005 # 1 decode unit = 16 ADC counts = 0.005 in/s (series-3 Normal range) + +# Blastware FFT-report ground truth: file → {channel: (dominant_hz, amplitude_ips)}. +# amplitude is None where the channel is at the noise floor (report amp 0.000/0.001) +# — the dominant frequency still matches exactly, but the amplitude isn't meaningful. +ORACLE = { + "N844LPGH.VV0W": {"Tran": (27.00, 0.018), "Vert": (26.75, 0.009), "Long": (26.50, 0.021), "MicL": (2.000, None)}, + "N844LPPR.3S0W": {"Tran": (30.75, None), "Vert": (46.75, None), "Long": (26.75, None), "MicL": (49.50, None)}, + "N844LQHB.ZT0W": {"Tran": (19.75, 0.040), "Vert": (26.50, 0.018), "Long": (26.50, 0.083), "MicL": (2.750, None)}, + "N844LQUE.T50W": {"Tran": (21.50, 0.080), "Vert": (14.25, 0.028), "Long": (28.50, 0.046), "MicL": (5.750, None)}, + "N844LR8W.790W": {"Tran": (31.00, None), "Vert": (31.00, None), "Long": (34.00, None), "MicL": (66.25, None)}, + "N844LRCO.G60W": {"Tran": (32.25, 0.009), "Vert": (32.00, 0.005), "Long": (32.00, 0.008), "MicL": (32.00, None)}, + "N844LRCW.F30W": {"Tran": (21.25, 0.010), "Vert": (42.25, 0.002), "Long": (21.25, 0.014), "MicL": (21.25, None)}, +} + + +def test_pure_sine_frequency_and_amplitude(): + # A pure sine at a bin-centre frequency (128 cycles over 4096 samples) has no + # leakage, so the single-sided 2/N normalisation returns the amplitude exactly. + sps, n, f0, amp = 1024.0, 4096, 32.0, 0.5 + x = amp * np.sin(2 * np.pi * f0 * np.arange(n) / sps) + freqs, amps = channel_spectrum(x, sps=sps, nfft=4096) + fpk, apk = dominant_frequency(freqs, amps) + assert fpk == 32.0 + assert abs(apk - amp) < 1e-3 + + +def test_bin_resolution_is_quarter_hz(): + freqs, _ = channel_spectrum(np.zeros(3328), sps=1024.0, nfft=4096) + assert abs((freqs[1] - freqs[0]) - 0.25) < 1e-9 + + +def test_empty_input(): + freqs, amps = channel_spectrum([]) + assert len(freqs) == 0 and len(amps) == 0 + + +def _spectra(fname): + raw = (FIXDIR / fname).read_bytes() + dec = decode_waveform_v2(raw[raw.find(b"STRT") + 21:]) + out = {} + for ch, samples in dec.items(): + ips = np.asarray(samples, float) * GEO_LSB + out[ch] = channel_spectrum(ips, sps=1024.0) + return out + + +def test_dominant_frequency_matches_blastware_exactly(): + misses = [] + for fname, chans in ORACLE.items(): + spectra = _spectra(fname) + for ch, (want_hz, _) in chans.items(): + got_hz, _ = dominant_frequency(*spectra[ch]) + if abs(got_hz - want_hz) > 0.25: + misses.append(f"{fname}:{ch} got {got_hz} want {want_hz}") + assert not misses, "dominant-frequency mismatches:\n" + "\n".join(misses) + + +def test_amplitude_matches_blastware(): + misses = [] + for fname, chans in ORACLE.items(): + spectra = _spectra(fname) + for ch, (_, want_amp) in chans.items(): + if want_amp is None: + continue + _, got_amp = dominant_frequency(*spectra[ch]) + if abs(got_amp - want_amp) > 0.0015: + misses.append(f"{fname}:{ch} got {got_amp:.4f} want {want_amp:.3f}") + assert not misses, "amplitude mismatches:\n" + "\n".join(misses) diff --git a/waveform_fft.py b/waveform_fft.py new file mode 100644 index 0000000..7e84565 --- /dev/null +++ b/waveform_fft.py @@ -0,0 +1,66 @@ +"""Blastware-compatible FFT of a decoded seismograph channel. + +Pure numpy; no I/O, no device or DB dependencies. Feed it a channel's decoded +samples **in the unit you want the amplitudes in** (e.g. in/s) and it returns the +single-sided amplitude spectrum that Blastware's *FFT Report* draws. + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events with +Blastware FFT reports as ground truth. The recipe reproduces Blastware's +**dominant frequency to the exact 0.25 Hz bin on all 28 channels** and the +amplitude to report precision: + + 1. remove the DC component (subtract the mean); **no window** — a window + smears the peak and measurably worsens the match, + 2. zero-pad to ``nfft`` (4096 → 0.25 Hz bins at 1024 sps — Blastware's + resolution), + 3. single-sided amplitude ``A[k] = 2·|X[k]| / N`` where ``N`` is the real + sample count (not ``nfft``). + +The compliance chart (USBM RI8507 / OSMRE) is this spectrum's ``(freq, amp)`` +points plotted against the regulatory limit curve; the #10 FFT view is the +spectrum itself. +""" +from __future__ import annotations + +import numpy as np + +BW_NFFT = 4096 # 0.25 Hz bins at 1024 sps — Blastware's FFT resolution +BW_FMIN = 2.0 # dominant-frequency search floor (Hz) +BW_FMAX = 250.0 # dominant-frequency search ceiling (Hz) + + +def channel_spectrum(samples, sps: float = 1024.0, nfft: int = BW_NFFT): + """Single-sided amplitude spectrum of one channel, Blastware-compatible. + + ``samples`` is a 1-D sequence in the desired amplitude unit (in/s). Returns + ``(freqs, amps)`` numpy arrays covering ``0 .. sps/2`` in ``sps/nfft`` steps. + + Records longer than ``nfft`` are truncated by the transform — untested + against Blastware for that case (real MiniMate Plus records are ≤ ~3.3 s, + well under 4096 samples at 1024 sps). + """ + x = np.asarray(samples, dtype=float) + n = x.size + if n == 0: + return np.empty(0), np.empty(0) + x = x - x.mean() # DC removal, no window + mag = np.abs(np.fft.rfft(x, nfft)) + freqs = np.fft.rfftfreq(nfft, 1.0 / sps) + amps = (2.0 / n) * mag # single-sided amplitude + return freqs, amps + + +def dominant_frequency(freqs, amps, fmin: float = BW_FMIN, fmax: float = BW_FMAX): + """Peak ``(frequency_hz, amplitude)`` of a spectrum within ``[fmin, fmax)``. + + Matches Blastware's "Dominant Frequency" — the largest spectral bin in the + reportable band (below 2 Hz is baseline/DC drift, above 250 Hz is noise). + """ + freqs = np.asarray(freqs) + amps = np.asarray(amps) + lo = int(np.searchsorted(freqs, fmin)) + hi = int(np.searchsorted(freqs, fmax)) + if hi <= lo: + return 0.0, 0.0 + k = lo + int(np.argmax(amps[lo:hi])) + return float(freqs[k]), float(amps[k])