// Radix-2 Cooley–Tukey FFT, in-place, iterative, bit-reversal ordered.
// Used by the RTL-SDR spectrum sweep. Length must be a power of two.

export type ComplexArray = { re: Float64Array; im: Float64Array }

export function nextPow2(n: number): number {
  let p = 1
  while (p < n) p <<= 1
  return p
}

/**
 * In-place iterative radix-2 FFT. `re`/`im` must have identical length
 * and that length must be a power of two.
 */
export function fftInPlace(re: Float64Array, im: Float64Array): void {
  const n = re.length
  if (n !== im.length) throw new Error('re/im length mismatch')
  if ((n & (n - 1)) !== 0) throw new Error('FFT length must be a power of two')
  if (n <= 1) return

  // Bit-reversal permutation.
  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) {
      const tr = re[i]
      re[i] = re[j]
      re[j] = tr
      const ti = im[i]
      im[i] = im[j]
      im[j] = ti
    }
  }

  for (let len = 2; len <= n; len <<= 1) {
    const half = 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 cwr = 1
      let cwi = 0
      for (let j = 0; j < half; j++) {
        const uR = re[i + j]
        const uI = im[i + j]
        const vR = re[i + j + half] * cwr - im[i + j + half] * cwi
        const vI = re[i + j + half] * cwi + im[i + j + half] * cwr
        re[i + j] = uR + vR
        im[i + j] = uI + vI
        re[i + j + half] = uR - vR
        im[i + j + half] = uI - vI
        const nwr = cwr * wr - cwi * wi
        cwi = cwr * wi + cwi * wr
        cwr = nwr
      }
    }
  }
}

/** Convenience wrapper returning magnitude spectrum (first n/2 bins). */
export function magnitudeSpectrum(re: Float64Array, im: Float64Array): Float64Array {
  const r = Float64Array.from(re)
  const i = Float64Array.from(im)
  fftInPlace(r, i)
  const half = r.length >> 1
  const mags = new Float64Array(half)
  for (let k = 0; k < half; k++) {
    mags[k] = Math.hypot(r[k], i[k])
  }
  return mags
}

/**
 * Power spectrum in dBFS-ish units from interleaved I/Q samples.
 * `iq` layout: [i0, q0, i1, q1, ...]. Returns n/2 dB values.
 */
export function powerSpectrumDb(iq: Float64Array): Float64Array {
  const n = nextPow2(iq.length >> 1)
  const re = new Float64Array(n)
  const im = new Float64Array(n)
  for (let k = 0; k < n && 2 * k + 1 < iq.length; k++) {
    re[k] = iq[2 * k]
    im[k] = iq[2 * k + 1]
  }
  fftInPlace(re, im)
  const half = n >> 1
  const out = new Float64Array(half)
  for (let k = 0; k < half; k++) {
    const power = (re[k] * re[k] + im[k] * im[k]) / (n * n)
    out[k] = 10 * Math.log10(power + 1e-12)
  }
  return out
}