import { timeStep } from './solver';
import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';

// Anti-aliased pickup. The solvers record the pickup once per output sample (point
// sampling), so motion above fs/2 would fold back into the audible band: a 10 mm pluck of
// the 5-bead string puts about 28% of its pickup power above 24 kHz. Instead, a run is
// solved with an internal sample rate R fs, and the pickup is low-pass filtered (a Kaiser-
// windowed sinc, passband to 0.45 fs, stopband from 0.5 fs, about 86 dB down) and kept
// every R-th sample. R is chosen so that the internal Nyquist frequency is above twice the
// system's highest linear frequency plus fs/2 (room for nonlinear combination tones),
// capped at the solver's steps per sample, so the time step does not shrink.

export function oversampling(sys: BeadSystem, exc: Excitation, opt: SimOptions): number {
  const { omegaMax, stepsPerSample } = timeStep(sys, exc.mallet, opt.fs, opt.cfl);
  const fmax = omegaMax / (2 * Math.PI);
  return Math.max(1, Math.min(stepsPerSample, Math.ceil((2 * fmax + opt.fs / 2) / opt.fs)));
}

// the options the solver runs with: R times the sample rate, delays in internal samples
export function oversampledOptions(opt: SimOptions, R: number): SimOptions {
  return { ...opt, fs: opt.fs * R, pickupDelay: opt.pickupDelay?.map((d) => d * R) };
}

const BETA = 8.6; // Kaiser window: about 86 dB stopband

// filter taps at the internal rate, symmetric about the middle, summing to 1
const kernels = new Map<number, Float64Array>();
export function lowpass(R: number): Float64Array {
  let h = kernels.get(R);
  if (h) return h;
  const df = 0.05 / R; // transition width, cycles per internal sample (0.45 to 0.5 of fs)
  const fc = 0.475 / R; // cutoff, cycles per internal sample
  let N = Math.ceil((86 - 7.95) / (14.36 * df)) + 1;
  if (N % 2 === 0) N++;
  const M = (N - 1) / 2;
  h = new Float64Array(N);
  const i0b = besselI0(BETA);
  let sum = 0;
  for (let i = 0; i < N; i++) {
    const x = i - M;
    const sinc = x === 0 ? 2 * fc : Math.sin(2 * Math.PI * fc * x) / (Math.PI * x);
    const r = x / M;
    h[i] = sinc * (besselI0(BETA * Math.sqrt(Math.max(0, 1 - r * r))) / i0b);
    sum += h[i];
  }
  for (let i = 0; i < N; i++) h[i] /= sum;
  kernels.set(R, h);
  return h;
}

function besselI0(x: number): number {
  let s = 1;
  let term = 1;
  for (let k = 1; k < 50; k++) {
    term *= (x / (2 * k)) ** 2;
    s += term;
    if (term < 1e-12 * s) break;
  }
  return s;
}

// Filter and decimate an internal-rate result to fs. complete: the whole run (the signal
// is zero after its end); otherwise a partial result, of which only the output samples
// whose filter window is covered are kept.
export function decimateResult(r: SimResult, R: number, complete: boolean): SimResult {
  if (R === 1) return { ...r, oversampling: 1 };
  const h = lowpass(R);
  const M = (h.length - 1) / 2;
  const x = r.audio;
  const nOut = complete ? Math.ceil(x.length / R) : Math.max(0, Math.floor((x.length - 1 - M) / R) + 1);
  const audio = new Float32Array(nOut);
  for (let i = 0; i < nOut; i++) {
    const c = i * R; // centre, in internal samples; the signal is zero before t = 0
    const a = Math.max(0, c - M);
    const b = Math.min(x.length - 1, c + M);
    let s = 0;
    for (let j = a; j <= b; j++) s += h[j - c + M] * x[j];
    audio[i] = s;
  }
  let contactForce: Float32Array | undefined;
  if (r.contactForce) {
    contactForce = new Float32Array(nOut);
    for (let i = 0; i < nOut; i++) contactForce[i] = r.contactForce[Math.min(i * R, r.contactForce.length - 1)];
  }
  return { ...r, audio, contactForce, fs: r.fs / R, stepsPerSample: r.stepsPerSample * R, oversampling: R };
}