path: root/tools/mq_editor/mq_extract.js

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

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

  // Get mono channel (mix to mono if stereo)
  const signal = getMono(audioBuffer);
  const numFrames = Math.floor((signal.length - fftSize) / hopSize);

  // Analyze frames
  const frames = [];
  for (let i = 0; i < numFrames; ++i) {
    const offset = i * hopSize;
    const frame = signal.slice(offset, offset + fftSize);
    const peaks = detectPeaks(frame, fftSize, sampleRate, threshold);
    const time = offset / sampleRate;
    frames.push({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;
}

// Get mono signal
function getMono(audioBuffer) {
  const data = audioBuffer.getChannelData(0);
  if (audioBuffer.numberOfChannels === 1) {
    return data;
  }

  // Mix to mono
  const left = audioBuffer.getChannelData(0);
  const right = audioBuffer.getChannelData(1);
  const mono = new Float32Array(left.length);
  for (let i = 0; i < left.length; ++i) {
    mono[i] = (left[i] + right[i]) * 0.5;
  }
  return mono;
}

// Detect peaks in FFT frame
function detectPeaks(frame, fftSize, sampleRate, thresholdDB) {
  // Apply Hann window
  const windowed = new Float32Array(fftSize);
  for (let i = 0; i < fftSize; ++i) {
    const w = 0.5 - 0.5 * Math.cos(2 * Math.PI * i / fftSize);
    windowed[i] = frame[i] * w;
  }

  // FFT (using built-in)
  const spectrum = realFFT(windowed);

  // Convert to magnitude dB
  const mag = new Float32Array(fftSize / 2);
  for (let i = 0; i < fftSize / 2; ++i) {
    const re = spectrum[i * 2];
    const im = spectrum[i * 2 + 1];
    const magLin = Math.sqrt(re * re + im * im);
    mag[i] = 20 * Math.log10(Math.max(magLin, 1e-10));
  }

  // 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};
}