// MQ Synthesizer
// Harmonic oscillator bank for sinusoidal synthesis, plus two-pole resonator mode

// Deterministic LCG PRNG
function randFloat(seed, min, max) {
  seed = (1664525 * seed + 1013904223) % 0x100000000;
  return min + (seed / 0x100000000) * (max - min);
}

// Build harmonic list from harmonics config.
// Fundamental (ratio=1.0, ampMult=1.0) is always first.
// Then harmonics at n*freq_mult for n=1,2,... with ampMult=decay^n (added on top).
function buildHarmonics(harmonics) {
  const decay    = Math.min(harmonics.decay ?? 0.0, 0.90);
  const freqMult = harmonics.freq_mult ?? 2.0;
  const result   = [{ ratio: 1.0, ampMult: 1.0 }]; // fundamental always
  if (decay > 0) {
    for (let n = 1; ; ++n) {
      const ampMult = Math.pow(decay, n);
      if (ampMult < 0.001) break;
      result.push({ ratio: n * freqMult, ampMult });
    }
  }
  return result;
}

// Synthesize audio from MQ partials
// partials: array of {freqCurve (with a0-a3 for amp), harmonics?, resonator?}
// harmonics: {decay, freq_mult, jitter, spread}
// resonator: {enabled, r, gainComp}  — two-pole resonator mode per partial
// integratePhase: true = accumulate 2π*f/SR per sample (correct for varying freq)
//                 false = 2π*f*t (simpler, only correct for constant freq)
// options.k1: LP coefficient in (0,1] — omit to bypass
// options.k2: HP coefficient in (0,1] — omit to bypass
function synthesizeMQ(partials, sampleRate, duration, integratePhase = true, options = {}) {
  const numSamples = Math.floor(sampleRate * duration);
  const pcm = new Float32Array(numSamples);

  const defaultHarmonics = {
    decay:     0.0,
    freq_mult: 1.0,
    jitter:    0.05,
    spread:    0.02
  };

  // Pre-build per-partial configs with fixed spread/jitter and phase accumulators
  const configs = [];
  for (let p = 0; p < partials.length; ++p) {
    const partial = partials[p];
    const fc = partial.freqCurve;

    if ((partial.resonator && partial.resonator.enabled) || options.forceResonator) {
      // --- Two-pole resonator mode ---
      // Driven by band-limited noise scaled by amp curve.
      // r controls pole radius (bandwidth): r→1 = narrow, r→0 = wide.
      // gainNorm = sqrt(1 - r²) normalises steady-state output power to ~A.
      const res      = partial.resonator || {};
      const r        = options.forceRGain ? clamp(options.globalR, 0, 0.9999)
                     : (res.r        != null ? clamp(res.r, 0, 0.9999) : 0.995);
      const gainComp = options.forceRGain ? options.globalGain
                     : (res.gainComp != null ? res.gainComp : 1.0);
      const gainNorm = Math.sqrt(Math.max(0, 1.0 - r * r));

      // Build harmonic list (jitter/spread not applied to resonator)
      const harm = partial.harmonics || defaultHarmonics;
      const harmonicList = buildHarmonics(harm);

      configs.push({
        mode: 'resonator',
        fc,
        r, gainComp, gainNorm,
        harmonicList,
        y1: new Float64Array(harmonicList.length),
        y2: new Float64Array(harmonicList.length),
        noiseSeed: ((p * 1664525 + 1013904223) & 0xFFFFFFFF) >>> 0
      });
    } else {
      // --- Sinusoidal (harmonic) mode ---
      const harm = partial.harmonics || defaultHarmonics;
      const spread = harm.spread ?? 0.0;
      const jitter = harm.jitter ?? 0.0;

      const harmonicList = buildHarmonics(harm);
      const replicaData = [];
      for (let h = 0; h < harmonicList.length; ++h) {
        const hc        = harmonicList[h];
        const spreadVal = randFloat(p * 67890 + h * 999, -spread, spread);
        const initPhase = randFloat(p * 67890 + h, 0.0, 1.0) * jitter * 2.0 * Math.PI;
        replicaData.push({ ratio: hc.ratio, ampMult: hc.ampMult, spread: spreadVal, phase: initPhase });
      }
      configs.push({ mode: 'sinusoid', fc, replicaData });
    }
  }

  for (let i = 0; i < numSamples; ++i) {
    const t = i / sampleRate;
    let sample = 0.0;

    for (let p = 0; p < configs.length; ++p) {
      const cfg = configs[p];
      const {fc} = cfg;

      if (cfg.mode === 'resonator') {
        if (t < fc.t0 || t > fc.t3) {
          cfg.y1.fill(0.0); cfg.y2.fill(0.0); continue;
        }

        const f0 = evalBezier(fc, t);
        const A  = evalBezierAmp(fc, t);

        // LCG noise excitation (deterministic per-partial, shared across harmonics)
        cfg.noiseSeed = (Math.imul(1664525, cfg.noiseSeed) + 1013904223) >>> 0;
        const noise   = cfg.noiseSeed / 0x100000000 * 2.0 - 1.0;

        for (let h = 0; h < cfg.harmonicList.length; ++h) {
          const hc    = cfg.harmonicList[h];
          const fh    = f0 * hc.ratio;
          const omega = 2.0 * Math.PI * fh / sampleRate;
          const b1    = 2.0 * cfg.r * Math.cos(omega);

          const x    = A * cfg.gainNorm * noise * hc.ampMult;
          const y    = b1 * cfg.y1[h] - cfg.r * cfg.r * cfg.y2[h] + x;
          cfg.y2[h]  = cfg.y1[h];
          cfg.y1[h]  = y;
          sample    += y * cfg.gainComp;
        }

      } else {
        if (t < fc.t0 || t > fc.t3) continue;

        const f0 = evalBezier(fc, t);
        const A0 = evalBezierAmp(fc, t);

        for (let h = 0; h < cfg.replicaData.length; ++h) {
          const rep = cfg.replicaData[h];
          const f   = f0 * rep.ratio * (1.0 + rep.spread);
          const A   = A0 * rep.ampMult;

          let phase;
          if (integratePhase) {
            rep.phase += 2.0 * Math.PI * f / sampleRate;
            phase = rep.phase;
          } else {
            phase = 2.0 * Math.PI * f * t + rep.phase;
          }

          sample += A * Math.sin(phase);
        }
      }
    }

    pcm[i] = sample;
  }

  // Post-synthesis filters (applied before normalization)
  // LP: y[n] = k1*x[n] + (1-k1)*y[n-1]  — options.k1 in (0,1], omit to bypass
  // HP: y[n] = k2*(y[n-1] + x[n] - x[n-1]) — options.k2 in (0,1], omit to bypass
  if (options.k1 != null) {
    const k1 = clamp(options.k1, 0, 1);
    let y = 0.0;
    for (let i = 0; i < numSamples; ++i) {
      y = k1 * pcm[i] + (1.0 - k1) * y;
      pcm[i] = y;
    }
  }
  if (options.k2 != null) {
    const k2 = clamp(options.k2, 0, 1);
    let y = 0.0, xp = 0.0;
    for (let i = 0; i < numSamples; ++i) {
      const x = pcm[i];
      y = k2 * (y + x - xp);
      xp = x;
      pcm[i] = y;
    }
  }

  // Normalize
  let maxAbs = 0;
  for (let i = 0; i < numSamples; ++i) maxAbs = Math.max(maxAbs, Math.abs(pcm[i]));
  if (maxAbs > 1.0) {
    for (let i = 0; i < numSamples; ++i) pcm[i] /= maxAbs;
  }

  return pcm;
}