/**
 * Phase maps: the angle phi_k applied to band k, as a function of the band
 * center frequency. The maps must vary smoothly from one band to the next,
 * otherwise neighbouring bands partly cancel where they overlap and the
 * magnitude response acquires notches.
 */

export type PhaseMapSpec =
  | {
      kind: 'random'
      /** Gaussian smoothing width across bands */
      smoothBands: number
      seed: number
      /** multiplies the map; 1 gives the spread of uniform random phases */
      strength: number
    }
  | {
      kind: 'delay'
      /** phi = -2 pi f tau: a pure delay, which must be inaudible */
      delayMs: number
    }
  | {
      kind: 'chirp'
      /** group delay grows linearly with frequency, reaching this value at the Nyquist frequency */
      delayAtNyquistMs: number
    }

export const defaultPhaseMap: PhaseMapSpec = { kind: 'random', smoothBands: 8, seed: 1, strength: 1 }

/** Small seeded generator (mulberry32) so a map can be reproduced from its seed. */
export function mulberry32(seed: number): () => number {
  let a = seed >>> 0
  return () => {
    a = (a + 0x6d2b79f5) >>> 0
    let t = a
    t = Math.imul(t ^ (t >>> 15), t | 1)
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61)
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296
  }
}

const UNIFORM_STD = Math.PI / Math.sqrt(3)

function smoothedRandom(n: number, smoothBands: number, seed: number, strength: number): Float64Array {
  const rnd = mulberry32(seed)
  const raw = new Float64Array(n)
  for (let i = 0; i < n; i++) raw[i] = (rnd() * 2 - 1) * Math.PI
  if (smoothBands <= 0) return raw.map((v) => v * strength)
  const radius = Math.ceil(4 * smoothBands)
  const kern = new Float64Array(2 * radius + 1)
  for (let j = -radius; j <= radius; j++) kern[j + radius] = Math.exp(-0.5 * (j / smoothBands) ** 2)
  const out = new Float64Array(n)
  for (let i = 0; i < n; i++) {
    let s = 0
    let wsum = 0
    for (let j = -radius; j <= radius; j++) {
      const m = i + j
      if (m < 0 || m >= n) continue
      const w = kern[j + radius]
      s += w * raw[m]
      wsum += w
    }
    out[i] = s / wsum
  }
  // Smoothing shrinks the spread; rescale so the map keeps the std of uniform phases.
  let mean = 0
  for (let i = 0; i < n; i++) mean += out[i]
  mean /= n
  let v = 0
  for (let i = 0; i < n; i++) v += (out[i] - mean) ** 2
  const std = Math.sqrt(v / Math.max(1, n - 1))
  const scale = std > 0 ? (strength * UNIFORM_STD) / std : 0
  for (let i = 0; i < n; i++) out[i] = (out[i] - mean) * scale
  return out
}

/** phi_k for every band center. */
export function phaseMap(spec: PhaseMapSpec, centers: Float64Array, fs: number): Float64Array {
  const n = centers.length
  switch (spec.kind) {
    case 'random':
      return smoothedRandom(n, spec.smoothBands, spec.seed, spec.strength)
    case 'delay': {
      const tau = spec.delayMs / 1000
      return centers.map((f) => -2 * Math.PI * f * tau)
    }
    case 'chirp': {
      // tau(f) = tauN f / fN  =>  phi(f) = -2 pi int tau df = -pi tauN f^2 / fN
      const tauN = spec.delayAtNyquistMs / 1000
      const fN = fs / 2
      return centers.map((f) => (-Math.PI * tauN * f * f) / fN)
    }
  }
}

/** Largest band-to-band slope of the map, as a group delay in seconds. Used to size the FFT padding. */
export function maxGroupDelaySec(phi: Float64Array, df: number): number {
  let m = 0
  for (let i = 1; i < phi.length; i++) {
    let d = Math.abs(phi[i] - phi[i - 1])
    d = Math.min(d, 2 * Math.PI) // a jump larger than 2 pi wraps to something smaller anyway
    if (d > m) m = d
  }
  return m / (2 * Math.PI * df)
}

export const wrapAngle = (x: number): number => {
  let y = x % (2 * Math.PI)
  if (y > Math.PI) y -= 2 * Math.PI
  if (y <= -Math.PI) y += 2 * Math.PI
  return y
}