path: root/tools/mq_editor/mq_extract.js

// MQ Extraction Algorithm
// McAulay-Quatieri sinusoidal analysis

// Extract partials from audio buffer
function extractPartials(params, stftCache) {
  const {fftSize, threshold, sampleRate} = params;

  // Analyze frames from cache
  const frames = [];
  const numFrames = stftCache.getNumFrames();
  for (let i = 0; i < numFrames; ++i) {
    const cachedFrame = stftCache.getFrameAtIndex(i);
    const squaredAmp = stftCache.getSquaredAmplitude(cachedFrame.time);
    const peaks = detectPeaks(squaredAmp, fftSize, sampleRate, threshold);
    frames.push({time: cachedFrame.time, peaks});
  }

  // Track trajectories
  const partials = trackPartials(frames, sampleRate);

  // Fit bezier curves
  for (const partial of partials) {
    partial.freqCurve = fitBezier(partial.times, partial.freqs);
    partial.ampCurve = fitBezier(partial.times, partial.amps);
  }

  return {partials, frames};
}

// Detect peaks in FFT frame (squaredAmp is pre-computed cached re*re+im*im)
function detectPeaks(squaredAmp, fftSize, sampleRate, thresholdDB) {
  // Convert squared amplitude to dB (10*log10 == 20*log10 of magnitude)
  const mag = new Float32Array(fftSize / 2);
  for (let i = 0; i < fftSize / 2; ++i) {
    mag[i] = 10 * Math.log10(Math.max(squaredAmp[i], 1e-20));
  }

  // Find local maxima above threshold
  const peaks = [];
  for (let i = 2; i < mag.length - 2; ++i) {
    if (mag[i] > thresholdDB &&
        mag[i] > mag[i-1] && mag[i] > mag[i-2] &&
        mag[i] > mag[i+1] && mag[i] > mag[i+2]) {

      // Parabolic interpolation for sub-bin accuracy
      const alpha = mag[i-1];
      const beta = mag[i];
      const gamma = mag[i+1];
      const p = 0.5 * (alpha - gamma) / (alpha - 2*beta + gamma);

      const binFreq = (i + p) * sampleRate / fftSize;
      const ampDB = beta - 0.25 * (alpha - gamma) * p;
      const ampLin = Math.pow(10, ampDB / 20);

      peaks.push({freq: binFreq, amp: ampLin});
    }
  }

  return peaks;
}

// Track partials across frames (birth/death/continuation)
function trackPartials(frames, sampleRate) {
  const partials = [];
  const activePartials = [];
  const candidatePartials = []; // Pre-birth candidates
  const trackingThresholdRatio = 0.05; // 5% frequency tolerance
  const minTrackingHz = 20; // Minimum 20 Hz
  const birthPersistence = 3; // Require 3 consecutive frames to birth
  const deathAge = 5; // Allow 5 frame gap before death
  const minPartialLength = 10; // Minimum 10 frames for valid partial

  for (const frame of frames) {
    const matched = new Set();

    // Match peaks to existing partials
    for (const partial of activePartials) {
      const lastFreq = partial.freqs[partial.freqs.length - 1];
      const threshold = Math.max(lastFreq * trackingThresholdRatio, minTrackingHz);

      let bestPeak = null;
      let bestDist = Infinity;

      for (let i = 0; i < frame.peaks.length; ++i) {
        if (matched.has(i)) continue;

        const peak = frame.peaks[i];
        const dist = Math.abs(peak.freq - lastFreq);

        if (dist < threshold && dist < bestDist) {
          bestPeak = peak;
          bestDist = dist;
          partial.matchIdx = i;
        }
      }

      if (bestPeak) {
        // Continuation
        partial.times.push(frame.time);
        partial.freqs.push(bestPeak.freq);
        partial.amps.push(bestPeak.amp);
        partial.age = 0;
        matched.add(partial.matchIdx);
      } else {
        // No match
        partial.age++;
      }
    }

    // Update candidate partials (pre-birth)
    for (let i = candidatePartials.length - 1; i >= 0; --i) {
      const candidate = candidatePartials[i];
      const lastFreq = candidate.freqs[candidate.freqs.length - 1];
      const threshold = Math.max(lastFreq * trackingThresholdRatio, minTrackingHz);

      let bestPeak = null;
      let bestDist = Infinity;

      for (let i = 0; i < frame.peaks.length; ++i) {
        if (matched.has(i)) continue;

        const peak = frame.peaks[i];
        const dist = Math.abs(peak.freq - lastFreq);

        if (dist < threshold && dist < bestDist) {
          bestPeak = peak;
          bestDist = dist;
          candidate.matchIdx = i;
        }
      }

      if (bestPeak) {
        candidate.times.push(frame.time);
        candidate.freqs.push(bestPeak.freq);
        candidate.amps.push(bestPeak.amp);
        matched.add(candidate.matchIdx);

        // Birth if persistent enough
        if (candidate.times.length >= birthPersistence) {
          activePartials.push(candidate);
          candidatePartials.splice(i, 1);
        }
      } else {
        // Candidate died, remove
        candidatePartials.splice(i, 1);
      }
    }

    // Create new candidate partials from unmatched peaks
    for (let i = 0; i < frame.peaks.length; ++i) {
      if (matched.has(i)) continue;

      const peak = frame.peaks[i];
      candidatePartials.push({
        times: [frame.time],
        freqs: [peak.freq],
        amps: [peak.amp],
        age: 0,
        matchIdx: -1
      });
    }

    // Death old partials
    for (let i = activePartials.length - 1; i >= 0; --i) {
      if (activePartials[i].age > deathAge) {
        // Move to finished if long enough
        if (activePartials[i].times.length >= minPartialLength) {
          partials.push(activePartials[i]);
        }
        activePartials.splice(i, 1);
      }
    }
  }

  // Finish remaining active partials
  for (const partial of activePartials) {
    if (partial.times.length >= minPartialLength) {
      partials.push(partial);
    }
  }

  return partials;
}

// Fit cubic bezier curve to trajectory using Catmull-Rom tangents.
// Estimates end tangents via forward/backward differences, then converts
// Hermite form to Bezier: B1 = P0 + m0*dt/3, B2 = P3 - m3*dt/3.
// This guarantees the curve passes through both endpoints.
function fitBezier(times, values) {
  const n = times.length - 1;
  const t0 = times[0],  v0 = values[0];
  const t3 = times[n],  v3 = values[n];
  const dt = t3 - t0;

  if (dt <= 0 || n === 0) {
    return {t0, v0, t1: t0, v1: v0, t2: t3, v2: v3, t3, v3};
  }

  // Catmull-Rom endpoint tangents (forward diff at start, backward at end)
  const dt0 = times[1] - times[0];
  const m0  = dt0 > 0 ? (values[1] - values[0]) / dt0 : 0;

  const dtn = times[n] - times[n - 1];
  const m3  = dtn > 0 ? (values[n] - values[n - 1]) / dtn : 0;

  // Hermite -> Bezier control points
  const v1 = v0 + m0 * dt / 3;
  const v2 = v3 - m3 * dt / 3;

  return {t0, v0, t1: t0 + dt / 3, v1, t2: t0 + 2 * dt / 3, v2, t3, v3};
}