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