/**
 * Phase maps: the angle phi_k applied to band k, as a function of the band
 * center frequency and, for the random map, optionally of time. Maps must
 * vary smoothly from one band to the next, otherwise neighbouring bands
 * partly cancel where they overlap and the magnitude response acquires
 * notches. A map that varies in time is frequency modulation of each band:
 * the instantaneous frequency shifts by (d phi / dt) / 2 pi.
 */

export type PhaseMapSpec =
  | {
      kind: 'random'
      /** Gaussian smoothing width across bands */
      smoothBands: number
      /** Gaussian smoothing width in time; null means the map does not change in time */
      timeScaleMs: number | null
      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, timeScaleMs: null, seed: 1, strength: 1 }

export const isTimeVarying = (spec: PhaseMapSpec): boolean => spec.kind === 'random' && spec.timeScaleMs !== null

/** 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
  }
}

/** Standard deviation of a phase drawn uniformly from (-pi, pi]. */
export const UNIFORM_STD = Math.PI / Math.sqrt(3)

/**
 * A random phase field on a coarse grid of (band, time) cells, smoothed with
 * Gaussians in both directions and interpolated bilinearly. Row r holds band
 * index r * dK and column c holds time c * dT.
 */
export type PhaseField = {
  values: Float32Array
  nK: number
  nT: number
  dK: number
  dT: number
  /** largest |d phi / dt| / 2 pi over the grid, in Hz: the largest frequency shift the field imposes */
  maxRateHz: number
}

/** Gaussian smoothing along one axis of a row-major nK x nT grid, with the kernel renormalized at the edges. */
function smoothAxis(v: Float32Array, nK: number, nT: number, alongTime: boolean, sigmaCells: number): void {
  if (sigmaCells < 0.3) return
  const radius = Math.ceil(4 * sigmaCells)
  const kern = new Float64Array(2 * radius + 1)
  for (let j = -radius; j <= radius; j++) kern[j + radius] = Math.exp(-0.5 * (j / sigmaCells) ** 2)
  const len = alongTime ? nT : nK
  const lines = alongTime ? nK : nT
  const stride = alongTime ? 1 : nT
  const lineStride = alongTime ? nT : 1
  const tmp = new Float64Array(len)
  for (let l = 0; l < lines; l++) {
    const base = l * lineStride
    for (let i = 0; i < len; i++) {
      let s = 0
      let wsum = 0
      const j0 = Math.max(-radius, -i)
      const j1 = Math.min(radius, len - 1 - i)
      for (let j = j0; j <= j1; j++) {
        const w = kern[j + radius]
        s += w * v[base + (i + j) * stride]
        wsum += w
      }
      tmp[i] = s / wsum
    }
    for (let i = 0; i < len; i++) v[base + i * stride] = tmp[i]
  }
}

export function randomField(nBands: number, durationSec: number, smoothBands: number, timeScaleSec: number | null, seed: number, strength: number): PhaseField {
  const sb = Math.max(0, smoothBands)
  // grid spacing: about a third of each smoothing scale, bounded so the grid stays small
  const dK = Math.max(1, Math.floor(sb / 3), Math.ceil(nBands / 1024))
  const nK = Math.ceil(Math.max(0, nBands - 1) / dK) + 1
  let nT = 1
  let dT = durationSec
  if (timeScaleSec !== null) {
    dT = Math.max(timeScaleSec / 3, durationSec / 4096)
    nT = Math.ceil(durationSec / dT) + 1
  }
  const v = new Float32Array(nK * nT)
  const rnd = mulberry32(seed)
  for (let i = 0; i < v.length; i++) v[i] = (rnd() * 2 - 1) * Math.PI
  smoothAxis(v, nK, nT, false, sb / dK)
  if (timeScaleSec !== null) smoothAxis(v, nK, nT, true, timeScaleSec / dT)
  // rescale so the field keeps the spread of uniform phases, times the strength
  let mean = 0
  for (let i = 0; i < v.length; i++) mean += v[i]
  mean /= v.length
  let s2 = 0
  for (let i = 0; i < v.length; i++) s2 += (v[i] - mean) ** 2
  const std = Math.sqrt(s2 / Math.max(1, v.length - 1))
  const scale = std > 0 ? (strength * UNIFORM_STD) / std : 0
  for (let i = 0; i < v.length; i++) v[i] = (v[i] - mean) * scale
  let maxRate = 0
  if (nT > 1) {
    for (let r = 0; r < nK; r++) {
      for (let c = 0; c + 1 < nT; c++) {
        const d = Math.abs(v[r * nT + c + 1] - v[r * nT + c])
        if (d > maxRate) maxRate = d
      }
    }
    maxRate /= 2 * Math.PI * dT
  }
  return { values: v, nK, nT, dK, dT, maxRateHz: maxRate }
}

/** The field for band k as a function of the grid time index (interpolated between grid rows). */
export function fieldRow(field: PhaseField, k: number): Float64Array {
  const { values, nK, nT, dK } = field
  const r = Math.min(nK - 1, k / dK)
  const r0 = Math.floor(r)
  const r1 = Math.min(nK - 1, r0 + 1)
  const w = r - r0
  const out = new Float64Array(nT)
  for (let c = 0; c < nT; c++) out[c] = (1 - w) * values[r0 * nT + c] + w * values[r1 * nT + c]
  return out
}

/** Linear interpolation of a grid row at arbitrary times (seconds). */
export function sampleRow(row: Float64Array, dT: number, times: ArrayLike<number>, out: Float64Array): Float64Array {
  const nT = row.length
  for (let i = 0; i < times.length; i++) {
    if (nT === 1) {
      out[i] = row[0]
      continue
    }
    const c = Math.min(nT - 1, Math.max(0, times[i] / dT))
    const c0 = Math.floor(c)
    const c1 = Math.min(nT - 1, c0 + 1)
    const w = c - c0
    out[i] = (1 - w) * row[c0] + w * row[c1]
  }
  return out
}

export function fieldAt(field: PhaseField, k: number, tSec: number): number {
  const row = fieldRow(field, k)
  return sampleRow(row, field.dT, [tSec], new Float64Array(1))[0]
}

/** phi_k for every band center, for maps that do not vary in time. */
export function phaseMap(spec: PhaseMapSpec, centers: Float64Array, fs: number): Float64Array {
  const n = centers.length
  switch (spec.kind) {
    case 'random': {
      const field = randomField(n, 1, spec.smoothBands, null, spec.seed, spec.strength)
      const out = new Float64Array(n)
      for (let k = 0; k < n; k++) out[k] = fieldAt(field, k, 0)
      return out
    }
    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: ArrayLike<number>, 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
}