"""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)