// In-place radix-2 complex FFT and a Hann window. The twiddle factors and
// the bit-reversal permutation are cached per transform length, so calling
// fft repeatedly at one length (as the spectrogram does) costs only the
// butterflies.

interface Tables {
  cos: Float32Array;
  sin: Float32Array;
  rev: Uint32Array;
}

const cache = new Map<number, Tables>();

function tables(n: number): Tables {
  let t = cache.get(n);
  if (t) return t;
  const half = n >> 1;
  const cos = new Float32Array(half);
  const sin = new Float32Array(half);
  for (let k = 0; k < half; k++) {
    const a = (-2 * Math.PI * k) / n;
    cos[k] = Math.cos(a);
    sin[k] = Math.sin(a);
  }
  const rev = new Uint32Array(n);
  const bits = 31 - Math.clz32(n);
  for (let i = 0; i < n; i++) {
    let r = 0;
    for (let b = 0; b < bits; b++) r |= ((i >> b) & 1) << (bits - 1 - b);
    rev[i] = r;
  }
  t = { cos, sin, rev };
  cache.set(n, t);
  return t;
}

/**
 * Forward complex FFT of (re, im) in place, X[k] = sum_n x[n] exp(-2 pi i k n / N).
 * The length must be a power of two and the two arrays must have equal length.
 */
export function fft(re: Float32Array, im: Float32Array): void {
  const n = re.length;
  if (n < 1 || (n & (n - 1)) !== 0) throw new Error(`fft: length ${n} is not a power of two`);
  if (im.length !== n) throw new Error(`fft: re (${n}) and im (${im.length}) differ in length`);
  if (n === 1) return;
  const { cos, sin, rev } = tables(n);
  for (let i = 0; i < n; i++) {
    const r = rev[i];
    if (r > i) {
      const tr = re[i];
      re[i] = re[r];
      re[r] = tr;
      const ti = im[i];
      im[i] = im[r];
      im[r] = ti;
    }
  }
  for (let size = 2; size <= n; size <<= 1) {
    const half = size >> 1;
    const step = n / size;
    for (let start = 0; start < n; start += size) {
      for (let k = 0, w = 0; k < half; k++, w += step) {
        const c = cos[w];
        const s = sin[w];
        const a = start + k;
        const b = a + half;
        const xr = re[b] * c - im[b] * s;
        const xi = re[b] * s + im[b] * c;
        re[b] = re[a] - xr;
        im[b] = im[a] - xi;
        re[a] += xr;
        im[a] += xi;
      }
    }
  }
}

/** Periodic Hann window of length n, w[k] = 0.5 (1 - cos(2 pi k / n)). */
export function hann(n: number): Float32Array {
  const w = new Float32Array(n);
  for (let k = 0; k < n; k++) w[k] = 0.5 - 0.5 * Math.cos((2 * Math.PI * k) / n);
  return w;
}