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
function fitBezier(times, values) {
  if (times.length < 4) {
    // Not enough points, just use linear segments
    return {
      t0: times[0], v0: values[0],
      t1: times[0], v1: values[0],
      t2: times[times.length-1], v2: values[values.length-1],
      t3: times[times.length-1], v3: values[values.length-1]
    };
  }

  // Fix endpoints
  const t0 = times[0];
  const t3 = times[times.length - 1];
  const v0 = values[0];
  const v3 = values[values.length - 1];

  // Solve for interior control points via least squares
  // Simplification: place at 1/3 and 2/3 positions
  const t1 = t0 + (t3 - t0) / 3;
  const t2 = t0 + 2 * (t3 - t0) / 3;

  // Find v1, v2 by evaluating at nearest data points
  let v1 = v0, v2 = v3;
  let minDist1 = Infinity, minDist2 = Infinity;

  for (let i = 0; i < times.length; ++i) {
    const dist1 = Math.abs(times[i] - t1);
    const dist2 = Math.abs(times[i] - t2);

    if (dist1 < minDist1) {
      minDist1 = dist1;
      v1 = values[i];
    }
    if (dist2 < minDist2) {
      minDist2 = dist2;
      v2 = values[i];
    }
  }

  return {t0, v0, t1, v1, t2, v2, t3, v3};
}