feat(mq_editor): replace replicas with harmonics model

- Fundamental f0 always synthesized; harmonics added at n*freq_mult
- decay^n amplitude rolloff per harmonic (capped at 0.90)
- Resonator mode also expanded across harmonics (per-harmonic y1/y2 state)
- UI: h.decay, h.freq (default 2.0), jitter, spread↑/↓ params
- Viewer: faint dotted harmonic bands with spread visualization
- Default freq_mult=2.0 (natural harmonic series)

handoff(Gemini): harmonics model complete, ready for next task

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

1 files changed, 68 insertions, 41 deletions

diff --git a/tools/mq_editor/mq_synth.js b/tools/mq_editor/mq_synth.js
index a9f387c..e5f7e1a 100644
--- a/tools/mq_editor/mq_synth.js
+++ b/tools/mq_editor/mq_synth.js
@@ -1,5 +1,5 @@
 // MQ Synthesizer
-// Replica oscillator bank for sinusoidal synthesis, plus two-pole resonator mode
+// Harmonic oscillator bank for sinusoidal synthesis, plus two-pole resonator mode
 
 // Deterministic LCG PRNG
 function randFloat(seed, min, max) {
@@ -7,9 +7,26 @@ function randFloat(seed, min, max) {
   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), replicas?, resonator?}
-// replicas:  {offsets, decay_alpha, jitter, spread_above, spread_below}
+// partials: array of {freqCurve (with a0-a3 for amp), harmonics?, resonator?}
+// harmonics: {decay, freq_mult, jitter, spread_above, spread_below}
 // 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)
@@ -19,15 +36,12 @@ function synthesizeMQ(partials, sampleRate, duration, integratePhase = true, opt
   const numSamples = Math.floor(sampleRate * duration);
   const pcm = new Float32Array(numSamples);
 
-  const jitterMult = options.disableJitter ? 0 : 1;
-  const spreadMult = options.disableSpread ? 0 : 1;
-
-  const defaultReplicas = {
-    offsets:       [1.0],
-    decay_alpha:   0.1,
-    jitter:        0.05,
-    spread_above:  0.02,
-    spread_below:  0.02
+  const defaultHarmonics = {
+    decay:        0.0,
+    freq_mult:    1.0,
+    jitter:       0.05,
+    spread_above: 0.02,
+    spread_below: 0.02
   };
 
   // Pre-build per-partial configs with fixed spread/jitter and phase accumulators
@@ -47,29 +61,36 @@ function synthesizeMQ(partials, sampleRate, duration, integratePhase = true, opt
       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,
-        y1: 0.0, y2: 0.0,
+        harmonicList,
+        y1: new Float64Array(harmonicList.length),
+        y2: new Float64Array(harmonicList.length),
         noiseSeed: ((p * 1664525 + 1013904223) & 0xFFFFFFFF) >>> 0
       });
     } else {
-      // --- Sinusoidal (replica) mode ---
-      const rep = partial.replicas != null ? partial.replicas : defaultReplicas;
-      const offsets      = rep.offsets      != null ? rep.offsets      : [1.0];
-      const decay_alpha  = rep.decay_alpha  != null ? rep.decay_alpha  : 0.0;
-      const jitter       = rep.jitter       != null ? rep.jitter       : 0.0;
-      const spread_above = rep.spread_above != null ? rep.spread_above : 0.0;
-      const spread_below = rep.spread_below != null ? rep.spread_below : 0.0;
+      // --- Sinusoidal (harmonic) mode ---
+      const harm = partial.harmonics || defaultHarmonics;
+      const spread_above = harm.spread_above ?? 0.0;
+      const spread_below = harm.spread_below ?? 0.0;
+      const jitter       = harm.jitter       ?? 0.0;
 
+      const harmonicList = buildHarmonics(harm);
       const replicaData = [];
-      for (let r = 0; r < offsets.length; ++r) {
-        const spread    = spreadMult * randFloat(p * 67890 + r * 999, -spread_below, spread_above);
-        const initPhase = randFloat(p * 67890 + r, 0.0, 1.0) * (jitter * jitterMult) * 2.0 * Math.PI;
-        replicaData.push({ratio: offsets[r], spread, phase: initPhase});
+      for (let h = 0; h < harmonicList.length; ++h) {
+        const hc        = harmonicList[h];
+        const spread    = randFloat(p * 67890 + h * 999, -spread_below, spread_above);
+        const initPhase = randFloat(p * 67890 + h, 0.0, 1.0) * jitter * 2.0 * Math.PI;
+        replicaData.push({ ratio: hc.ratio, ampMult: hc.ampMult, spread, phase: initPhase });
       }
-      configs.push({ mode: 'sinusoid', fc, decay_alpha, replicaData });
+      configs.push({ mode: 'sinusoid', fc, replicaData });
     }
   }
 
@@ -82,34 +103,40 @@ function synthesizeMQ(partials, sampleRate, duration, integratePhase = true, opt
       const {fc} = cfg;
 
       if (cfg.mode === 'resonator') {
-        if (t < fc.t0 || t > fc.t3) { cfg.y1 = 0.0; cfg.y2 = 0.0; continue; }
+        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);
-        const omega = 2.0 * Math.PI * f0 / sampleRate;
-        const b1    = 2.0 * cfg.r * Math.cos(omega);
+        const f0 = evalBezier(fc, t);
+        const A  = evalBezierAmp(fc, t);
 
-        // LCG noise excitation (deterministic per-partial)
+        // 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;
 
-        const x  = A * cfg.gainNorm * noise;
-        const y  = b1 * cfg.y1 - cfg.r * cfg.r * cfg.y2 + x;
-        cfg.y2   = cfg.y1;
-        cfg.y1   = y;
-        sample  += y * cfg.gainComp;
+        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);
-        const {decay_alpha, replicaData} = cfg;
 
-        for (let r = 0; r < replicaData.length; ++r) {
-          const rep = replicaData[r];
-          const f = f0 * rep.ratio * (1.0 + rep.spread);
-          const A = A0 * Math.exp(-decay_alpha * Math.abs(f - f0));
+        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) {