/**
 * A plain radix-2 FFT, enough for the spectrum view: Hann-windowed
 * magnitude of the first power-of-two chunk of the trace, in dB below the
 * peak. Nothing here is performance-critical; a quarter-million-point
 * transform is well under a frame.
 */

/** In-place complex FFT, re/im interleaved as two arrays, n a power of 2. */
export function fft(re: Float64Array, im: Float64Array): void {
  const n = re.length;
  for (let i = 1, j = 0; i < n; i++) {
    let bit = n >> 1;
    for (; j & bit; bit >>= 1) j ^= bit;
    j ^= bit;
    if (i < j) {
      [re[i], re[j]] = [re[j], re[i]];
      [im[i], im[j]] = [im[j], im[i]];
    }
  }
  for (let len = 2; len <= n; len <<= 1) {
    const ang = (-2 * Math.PI) / len;
    const wr = Math.cos(ang);
    const wi = Math.sin(ang);
    for (let i = 0; i < n; i += len) {
      let cr = 1;
      let ci = 0;
      for (let k = 0; k < len / 2; k++) {
        const ur = re[i + k];
        const ui = im[i + k];
        const vr = re[i + k + len / 2] * cr - im[i + k + len / 2] * ci;
        const vi = re[i + k + len / 2] * ci + im[i + k + len / 2] * cr;
        re[i + k] = ur + vr;
        im[i + k] = ui + vi;
        re[i + k + len / 2] = ur - vr;
        im[i + k + len / 2] = ui - vi;
        const nr = cr * wr - ci * wi;
        ci = cr * wi + ci * wr;
        cr = nr;
      }
    }
  }
}

export interface Spectrum {
  /** Magnitude in dB relative to the strongest bin (all <= 0). */
  db: Float32Array;
  /** Hz per bin. */
  binHz: number;
}

/** The spectrum of a trace sampled at 1/dt: Hann window over the first
 *  power-of-two chunk (at most 2^18 samples). */
export function spectrumOf(trace: Float32Array, dt: number): Spectrum {
  let n = 1;
  while (n * 2 <= Math.min(trace.length, 1 << 18)) n *= 2;
  const re = new Float64Array(n);
  const im = new Float64Array(n);
  for (let i = 0; i < n; i++) {
    re[i] = trace[i] * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (n - 1)));
  }
  fft(re, im);
  const half = n / 2;
  const mag = new Float64Array(half);
  let peak = 0;
  for (let i = 0; i < half; i++) {
    mag[i] = Math.hypot(re[i], im[i]);
    if (mag[i] > peak) peak = mag[i];
  }
  const db = new Float32Array(half);
  for (let i = 0; i < half; i++) {
    db[i] = peak > 0 ? 20 * Math.log10(Math.max(mag[i] / peak, 1e-8)) : -160;
  }
  return { db, binHz: 1 / (n * dt) };
}