summaryrefslogtreecommitdiff
path: root/tools/mq_editor/mq_extract.js
diff options
context:
space:
mode:
Diffstat (limited to 'tools/mq_editor/mq_extract.js')
-rw-r--r--tools/mq_editor/mq_extract.js360
1 files changed, 299 insertions, 61 deletions
diff --git a/tools/mq_editor/mq_extract.js b/tools/mq_editor/mq_extract.js
index 97fbb00..18897fb 100644
--- a/tools/mq_editor/mq_extract.js
+++ b/tools/mq_editor/mq_extract.js
@@ -9,28 +9,48 @@ function extractPartials(params, stftCache) {
const frames = [];
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, freqWeight, prominence);
+ const squaredAmp = cachedFrame.squaredAmplitude;
+ const phase = cachedFrame.phase;
+ const peaks = detectPeaks(squaredAmp, phase, fftSize, sampleRate, threshold, freqWeight, prominence);
frames.push({time: cachedFrame.time, peaks});
}
- const partials = trackPartials(frames);
+ const partials = trackPartials(frames, params);
// Second pass: extend partials leftward to recover onset frames
expandPartialsLeft(partials, frames);
for (const partial of partials) {
partial.freqCurve = fitBezier(partial.times, partial.freqs);
- partial.ampCurve = fitBezier(partial.times, partial.amps);
+ const ac = fitBezier(partial.times, partial.amps);
+ partial.freqCurve.a0 = ac.v0; partial.freqCurve.a1 = ac.v1;
+ partial.freqCurve.a2 = ac.v2; partial.freqCurve.a3 = ac.v3;
}
return {partials, frames};
}
+// Helper to interpolate phase via quadratic formula on unwrapped neighbors.
+// This provides a more accurate phase estimate at the sub-bin peak location.
+function phaseInterp(p_minus, p_center, p_plus, p_frac) {
+ // unwrap neighbors relative to center
+ let dp_minus = p_minus - p_center;
+ while (dp_minus > Math.PI) dp_minus -= 2 * Math.PI;
+ while (dp_minus < -Math.PI) dp_minus += 2 * Math.PI;
+
+ let dp_plus = p_plus - p_center;
+ while (dp_plus > Math.PI) dp_plus -= 2 * Math.PI;
+ while (dp_plus < -Math.PI) dp_plus += 2 * Math.PI;
+
+ const p_interp = p_center + (dp_plus - dp_minus) * p_frac * 0.5 + (dp_plus + dp_minus) * 0.5 * p_frac * p_frac;
+ return p_interp;
+}
+
// Detect spectral peaks via local maxima + parabolic interpolation
// squaredAmp: pre-computed re*re+im*im per bin
+// phase: pre-computed atan2(im,re) per bin
// freqWeight: if true, weight by f before peak detection (f * Power(f))
-function detectPeaks(squaredAmp, fftSize, sampleRate, thresholdDB, freqWeight, prominenceDB = 0) {
+function detectPeaks(squaredAmp, phase, fftSize, sampleRate, thresholdDB, freqWeight, prominenceDB = 0) {
const mag = new Float32Array(fftSize / 2);
const binHz = sampleRate / fftSize;
for (let i = 0; i < fftSize / 2; ++i) {
@@ -62,56 +82,84 @@ function detectPeaks(squaredAmp, fftSize, sampleRate, thresholdDB, freqWeight, p
if (mag[i] - valley < prominenceDB) continue;
}
- // Parabolic interpolation for sub-bin accuracy
+ // Parabolic interpolation for sub-bin accuracy on frequency, amplitude, and phase
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 p_phase = phaseInterp(phase[i-1], phase[i], phase[i+1], p);
const freq = (i + p) * sampleRate / fftSize;
const ampDB = beta - 0.25 * (alpha - gamma) * p;
- peaks.push({freq, amp: Math.pow(10, ampDB / 20)});
+ peaks.push({freq, amp: Math.pow(10, ampDB / 20), phase: p_phase});
}
}
return peaks;
}
-// Track partials across frames (birth/death/continuation)
-function trackPartials(frames) {
+// Helper to compute shortest angle difference (e.g., between -pi and pi)
+function normalizeAngle(angle) {
+ return angle - 2 * Math.PI * Math.floor((angle + Math.PI) / (2 * Math.PI));
+}
+
+// Find best matching peak for a predicted freq/phase. Returns {bestIdx, bestCost}.
+function findBestPeak(peaks, matched, predictedFreq, predictedPhase, tol, phaseErrorWeight) {
+ let bestIdx = -1, bestCost = Infinity;
+ for (let i = 0; i < peaks.length; ++i) {
+ if (matched.has(i)) continue;
+ const pk = peaks[i];
+ const freqError = Math.abs(pk.freq - predictedFreq);
+ if (freqError > tol) continue;
+ const phaseError = Math.abs(normalizeAngle(pk.phase - predictedPhase));
+ const cost = freqError + phaseErrorWeight * phaseError * predictedFreq;
+ if (cost < bestCost) { bestCost = cost; bestIdx = i; }
+ }
+ return { bestIdx, bestCost };
+}
+
+// Track partials across frames using phase coherence for robust matching.
+function trackPartials(frames, params) {
+ const {
+ sampleRate, hopSize,
+ birthPersistence = 3,
+ deathAge = 5,
+ minLength = 10,
+ phaseErrorWeight = 2.0
+ } = params;
const partials = [];
const activePartials = [];
const candidates = []; // pre-birth
const trackingRatio = 0.05; // 5% frequency tolerance
const minTrackingHz = 20;
- const birthPersistence = 3; // frames before partial is born
- const deathAge = 5; // frames without match before death
- const minLength = 10; // frames required to keep partial
for (const frame of frames) {
const matched = new Set();
- // Continue active partials
+ // --- Continue active partials ---
for (const partial of activePartials) {
const lastFreq = partial.freqs[partial.freqs.length - 1];
+ const lastPhase = partial.phases[partial.phases.length - 1];
const velocity = partial.velocity || 0;
- const predicted = lastFreq + velocity;
+ const predictedFreq = lastFreq + velocity;
- const tol = Math.max(lastFreq * trackingRatio, minTrackingHz);
- let bestIdx = -1, bestDist = Infinity;
-
- for (let i = 0; i < frame.peaks.length; ++i) {
- if (matched.has(i)) continue;
- const dist = Math.abs(frame.peaks[i].freq - predicted);
- if (dist < tol && dist < bestDist) { bestDist = dist; bestIdx = i; }
- }
+ // Predict phase for the current frame based on the last frame's frequency.
+ // Multiply by (age+1) to account for frames missed during a gap.
+ const phaseAdvance = 2 * Math.PI * lastFreq * (partial.age + 1) * hopSize / sampleRate;
+ const predictedPhase = lastPhase + phaseAdvance;
+
+ const tol = Math.max(predictedFreq * trackingRatio, minTrackingHz);
+ // Find the peak in the new frame with the lowest cost (freq + phase error).
+ const { bestIdx } = findBestPeak(frame.peaks, matched, predictedFreq, predictedPhase, tol, phaseErrorWeight);
if (bestIdx >= 0) {
const pk = frame.peaks[bestIdx];
partial.times.push(frame.time);
partial.freqs.push(pk.freq);
partial.amps.push(pk.amp);
+ partial.phases.push(pk.phase);
partial.age = 0;
partial.velocity = pk.freq - lastFreq;
matched.add(bestIdx);
@@ -120,52 +168,54 @@ function trackPartials(frames) {
}
}
- // Advance candidates
+ // --- Advance candidates ---
for (let i = candidates.length - 1; i >= 0; --i) {
const cand = candidates[i];
const lastFreq = cand.freqs[cand.freqs.length - 1];
+ const lastPhase = cand.phases[cand.phases.length - 1];
const velocity = cand.velocity || 0;
- const predicted = lastFreq + velocity;
+ const predictedFreq = lastFreq + velocity;
- const tol = Math.max(lastFreq * trackingRatio, minTrackingHz);
- let bestIdx = -1, bestDist = Infinity;
+ // Candidates die on first miss so age is always 0 here, but kept consistent.
+ const phaseAdvance = 2 * Math.PI * lastFreq * hopSize / sampleRate;
+ const predictedPhase = lastPhase + phaseAdvance;
- for (let j = 0; j < frame.peaks.length; ++j) {
- if (matched.has(j)) continue;
- const dist = Math.abs(frame.peaks[j].freq - predicted);
- if (dist < tol && dist < bestDist) { bestDist = dist; bestIdx = j; }
- }
+ const tol = Math.max(predictedFreq * trackingRatio, minTrackingHz);
+ const { bestIdx } = findBestPeak(frame.peaks, matched, predictedFreq, predictedPhase, tol, phaseErrorWeight);
if (bestIdx >= 0) {
const pk = frame.peaks[bestIdx];
cand.times.push(frame.time);
cand.freqs.push(pk.freq);
cand.amps.push(pk.amp);
+ cand.phases.push(pk.phase);
cand.velocity = pk.freq - lastFreq;
matched.add(bestIdx);
+ // "graduate" a candidate to a full partial
if (cand.times.length >= birthPersistence) {
activePartials.push(cand);
candidates.splice(i, 1);
}
} else {
- candidates.splice(i, 1);
+ candidates.splice(i, 1); // kill candidate
}
}
- // Spawn candidates from unmatched peaks
+ // --- Spawn new candidates from unmatched peaks ---
for (let i = 0; i < frame.peaks.length; ++i) {
if (matched.has(i)) continue;
const pk = frame.peaks[i];
candidates.push({
times: [frame.time],
freqs: [pk.freq],
- amps: [pk.amp],
+ amps: [pk.amp],
+ phases: [pk.phase],
age: 0,
velocity: 0
});
}
- // Kill aged-out partials
+ // --- Kill aged-out partials ---
for (let i = activePartials.length - 1; i >= 0; --i) {
if (activePartials[i].age > deathAge) {
if (activePartials[i].times.length >= minLength) partials.push(activePartials[i]);
@@ -174,7 +224,7 @@ function trackPartials(frames) {
}
}
- // Collect remaining active partials
+ // --- Collect remaining active partials ---
for (const partial of activePartials) {
if (partial.times.length >= minLength) partials.push(partial);
}
@@ -193,6 +243,8 @@ function expandPartialsLeft(partials, frames) {
for (let i = 0; i < frames.length; ++i) timeToIdx.set(frames[i].time, i);
for (const partial of partials) {
+ if (!partial.phases) partial.phases = []; // Ensure old partials have phase array
+
let startIdx = timeToIdx.get(partial.times[0]);
if (startIdx == null || startIdx === 0) continue;
@@ -213,6 +265,7 @@ function expandPartialsLeft(partials, frames) {
partial.times.unshift(frame.time);
partial.freqs.unshift(pk.freq);
partial.amps.unshift(pk.amp);
+ partial.phases.unshift(pk.phase);
}
}
}
@@ -283,7 +336,200 @@ function autodetectSpread(partial, stftCache, fftSize, sampleRate) {
return {spread_above, spread_below};
}
-// Fit cubic bezier to trajectory using least-squares for inner control points
+// Track a single partial starting from a (time, freq) seed position.
+// Snaps to nearest spectral peak, then tracks forward and backward.
+// Returns a partial object (with freqCurve), or null if no peak found near seed.
+function trackFromSeed(frames, seedTime, seedFreq, params) {
+ if (!frames || frames.length === 0) return null;
+
+ // Find nearest frame to seedTime
+ let seedFrameIdx = 0;
+ let bestDt = Infinity;
+ for (let i = 0; i < frames.length; ++i) {
+ const dt = Math.abs(frames[i].time - seedTime);
+ if (dt < bestDt) { bestDt = dt; seedFrameIdx = i; }
+ }
+
+ // Snap to nearest spectral peak within 10% freq tolerance
+ const seedFrame = frames[seedFrameIdx];
+ const snapTol = Math.max(seedFreq * 0.10, 50);
+ let seedPeak = null, bestDist = snapTol;
+ for (const pk of seedFrame.peaks) {
+ const d = Math.abs(pk.freq - seedFreq);
+ if (d < bestDist) { bestDist = d; seedPeak = pk; }
+ }
+ if (!seedPeak) return null;
+
+ const { hopSize, sampleRate, deathAge = 5, phaseErrorWeight = 2.0 } = params;
+ const trackingRatio = 0.05;
+ const minTrackingHz = 20;
+
+ // Forward pass from seed frame
+ const times = [seedFrame.time];
+ const freqs = [seedPeak.freq];
+ const amps = [seedPeak.amp];
+ const phases = [seedPeak.phase];
+
+ let fwdFreq = seedPeak.freq, fwdPhase = seedPeak.phase, fwdVel = 0, fwdAge = 0;
+ for (let i = seedFrameIdx + 1; i < frames.length; ++i) {
+ const predicted = fwdFreq + fwdVel;
+ const predPhase = fwdPhase + 2 * Math.PI * fwdFreq * (fwdAge + 1) * hopSize / sampleRate;
+ const tol = Math.max(predicted * trackingRatio, minTrackingHz);
+ const { bestIdx } = findBestPeak(frames[i].peaks, new Set(), predicted, predPhase, tol, phaseErrorWeight);
+ if (bestIdx >= 0) {
+ const pk = frames[i].peaks[bestIdx];
+ times.push(frames[i].time);
+ freqs.push(pk.freq);
+ amps.push(pk.amp);
+ phases.push(pk.phase);
+ fwdVel = pk.freq - fwdFreq;
+ fwdFreq = pk.freq; fwdPhase = pk.phase; fwdAge = 0;
+ } else {
+ fwdAge++;
+ if (fwdAge > deathAge) break;
+ }
+ }
+
+ // Backward pass from seed frame
+ const bwdTimes = [], bwdFreqs = [], bwdAmps = [], bwdPhases = [];
+ let bwdFreq = seedPeak.freq, bwdAge = 0;
+ for (let i = seedFrameIdx - 1; i >= 0; --i) {
+ const tol = Math.max(bwdFreq * trackingRatio, minTrackingHz);
+ let bestIdx = -1, bDist = tol;
+ for (let j = 0; j < frames[i].peaks.length; ++j) {
+ const d = Math.abs(frames[i].peaks[j].freq - bwdFreq);
+ if (d < bDist) { bDist = d; bestIdx = j; }
+ }
+ if (bestIdx >= 0) {
+ const pk = frames[i].peaks[bestIdx];
+ bwdTimes.unshift(frames[i].time);
+ bwdFreqs.unshift(pk.freq);
+ bwdAmps.unshift(pk.amp);
+ bwdPhases.unshift(pk.phase);
+ bwdFreq = pk.freq; bwdAge = 0;
+ } else {
+ bwdAge++;
+ if (bwdAge > deathAge) break;
+ }
+ }
+
+ const allTimes = [...bwdTimes, ...times];
+ const allFreqs = [...bwdFreqs, ...freqs];
+ const allAmps = [...bwdAmps, ...amps];
+ const allPhases = [...bwdPhases, ...phases];
+
+ if (allTimes.length < 2) return null;
+
+ const freqCurve = fitBezier(allTimes, allFreqs);
+ const ac = fitBezier(allTimes, allAmps);
+ freqCurve.a0 = ac.v0; freqCurve.a1 = ac.v1;
+ freqCurve.a2 = ac.v2; freqCurve.a3 = ac.v3;
+
+ return {
+ times: allTimes, freqs: allFreqs, amps: allAmps, phases: allPhases,
+ muted: false, freqCurve,
+ replicas: { decay_alpha: 0.1, jitter: 0.05, spread_above: 0.02, spread_below: 0.02 },
+ };
+}
+
+// Track an iso-energy contour starting from (seedTime, seedFreq).
+// Instead of following spectral peaks, follows where energy ≈ seedEnergy.
+// Useful for broad/diffuse bass regions with no detectable peaks.
+// Returns a partial with large default spread, or null if seed energy is zero.
+function trackIsoContour(stftCache, seedTime, seedFreq, params) {
+ const { sampleRate, deathAge = 8 } = params;
+ const numFrames = stftCache.getNumFrames();
+ const fftSize = stftCache.fftSize;
+ const binHz = sampleRate / fftSize;
+ const halfBins = fftSize / 2;
+
+ // Find seed frame
+ let seedFrameIdx = 0, bestDt = Infinity;
+ for (let i = 0; i < numFrames; ++i) {
+ const dt = Math.abs(stftCache.getFrameAtIndex(i).time - seedTime);
+ if (dt < bestDt) { bestDt = dt; seedFrameIdx = i; }
+ }
+
+ const seedFrame = stftCache.getFrameAtIndex(seedFrameIdx);
+ const seedBin = Math.max(1, Math.min(halfBins - 2, Math.round(seedFreq / binHz)));
+ const targetSq = seedFrame.squaredAmplitude[seedBin];
+ if (targetSq <= 0) return null;
+ const targetDB = 10 * Math.log10(targetSq);
+
+ const trackingRatio = 0.15; // larger search window than peak tracker
+ const minTrackHz = 30;
+ const maxDbDev = 15; // dB: declare miss if nothing within this range
+
+ // Find bin minimizing |dB(b) - targetDB| near refBin, with mild position bias.
+ function findContourBin(sq, refBin) {
+ const tol = Math.max(refBin * binHz * trackingRatio, minTrackHz);
+ const tolBins = Math.ceil(tol / binHz);
+ const lo = Math.max(1, refBin - tolBins);
+ const hi = Math.min(halfBins - 2, refBin + tolBins);
+ let bestBin = -1, bestCost = Infinity;
+ for (let b = lo; b <= hi; ++b) {
+ const dE = Math.abs(10 * Math.log10(Math.max(sq[b], 1e-20)) - targetDB);
+ if (dE > maxDbDev) continue;
+ const dPos = Math.abs(b - refBin) / Math.max(1, tolBins);
+ const cost = dE + 3 * dPos; // energy match dominates, position breaks ties
+ if (cost < bestCost) { bestCost = cost; bestBin = b; }
+ }
+ return bestBin;
+ }
+
+ const times = [seedFrame.time];
+ const freqs = [seedBin * binHz];
+ const amps = [Math.sqrt(Math.max(0, targetSq))];
+
+ // Forward pass
+ let fwdBin = seedBin, fwdAge = 0;
+ for (let i = seedFrameIdx + 1; i < numFrames; ++i) {
+ const frame = stftCache.getFrameAtIndex(i);
+ const b = findContourBin(frame.squaredAmplitude, fwdBin);
+ if (b >= 0) {
+ times.push(frame.time);
+ freqs.push(b * binHz);
+ amps.push(Math.sqrt(Math.max(0, frame.squaredAmplitude[b])));
+ fwdBin = b; fwdAge = 0;
+ } else { if (++fwdAge > deathAge) break; }
+ }
+
+ // Backward pass
+ const bwdTimes = [], bwdFreqs = [], bwdAmps = [];
+ let bwdBin = seedBin, bwdAge = 0;
+ for (let i = seedFrameIdx - 1; i >= 0; --i) {
+ const frame = stftCache.getFrameAtIndex(i);
+ const b = findContourBin(frame.squaredAmplitude, bwdBin);
+ if (b >= 0) {
+ bwdTimes.unshift(frame.time);
+ bwdFreqs.unshift(b * binHz);
+ bwdAmps.unshift(Math.sqrt(Math.max(0, frame.squaredAmplitude[b])));
+ bwdBin = b; bwdAge = 0;
+ } else { if (++bwdAge > deathAge) break; }
+ }
+
+ const allTimes = [...bwdTimes, ...times];
+ const allFreqs = [...bwdFreqs, ...freqs];
+ const allAmps = [...bwdAmps, ...amps];
+ if (allTimes.length < 2) return null;
+
+ const freqCurve = fitBezier(allTimes, allFreqs);
+ const ac = fitBezier(allTimes, allAmps);
+ freqCurve.a0 = ac.v0; freqCurve.a1 = ac.v1;
+ freqCurve.a2 = ac.v2; freqCurve.a3 = ac.v3;
+
+ return {
+ times: allTimes, freqs: allFreqs, amps: allAmps,
+ phases: new Array(allTimes.length).fill(0),
+ muted: false, freqCurve,
+ replicas: { decay_alpha: 0.1, jitter: 0.05, spread_above: 0.15, spread_below: 0.15 },
+ };
+}
+
+// Fit interpolating curve to trajectory via least-squares for inner control point values.
+// Inner knots fixed at u=1/3 and u=2/3 (t = t0+dt/3, t0+2*dt/3).
+// The curve passes through all 4 control points (Lagrange interpolation).
+// TODO: support arbitrary number of inner control points
function fitBezier(times, values) {
const n = times.length - 1;
const t0 = times[0], v0 = values[0];
@@ -291,43 +537,35 @@ function fitBezier(times, values) {
const dt = t3 - t0;
if (dt <= 1e-9 || n < 2) {
- // Linear fallback for too few points or zero duration
return {t0, v0, t1: t0 + dt / 3, v1: v0 + (v3 - v0) / 3, t2: t0 + 2 * dt / 3, v2: v0 + 2 * (v3 - v0) / 3, t3, v3};
}
- // Least squares solve for v1, v2
- // Bezier: B(u) = (1-u)^3*v0 + 3(1-u)^2*u*v1 + 3(1-u)*u^2*v2 + u^3*v3
- // Target_i = val_i - (1-u)^3*v0 - u^3*v3
- // Model_i = A_i*v1 + B_i*v2
- // A_i = 3(1-u)^2*u
- // B_i = 3(1-u)*u^2
+ // Lagrange basis with inner knots at u1=1/3, u2=2/3
+ // l1(u) = u*(u-2/3)*(u-1) / ((1/3)*(1/3-2/3)*(1/3-1)) = 13.5*u*(u-2/3)*(u-1)
+ // l2(u) = u*(u-1/3)*(u-1) / ((2/3)*(2/3-1/3)*(2/3-1)) = -13.5*u*(u-1/3)*(u-1)
+ // l0(u) = (u-1/3)*(u-2/3)*(u-1) / ((-1/3)*(-2/3)*(-1)) = -4.5*(u-1/3)*(u-2/3)*(u-1)
+ // l3(u) = u*(u-1/3)*(u-2/3) / ((2/3)*(1/3)) = 4.5*u*(u-1/3)*(u-2/3)
+ // Least-squares: minimize Σ(l1*v1 + l2*v2 - target_i)^2
+ // target_i = values[i] - l0*v0 - l3*v3
let sA2 = 0, sB2 = 0, sAB = 0, sAT = 0, sBT = 0;
for (let i = 0; i <= n; ++i) {
- const u = (times[i] - t0) / dt;
- const u2 = u * u;
- const u3 = u2 * u;
- const invU = 1.0 - u;
- const invU2 = invU * invU;
- const invU3 = invU2 * invU;
-
- const A = 3 * invU2 * u;
- const B = 3 * invU * u2;
- const target = values[i] - (invU3 * v0 + u3 * v3);
-
- sA2 += A * A;
- sB2 += B * B;
- sAB += A * B;
- sAT += A * target;
- sBT += B * target;
+ const u = (times[i] - t0) / dt;
+ const l0 = -4.5 * (u - 1/3) * (u - 2/3) * (u - 1);
+ const l1 = 13.5 * u * (u - 2/3) * (u - 1);
+ const l2 = -13.5 * u * (u - 1/3) * (u - 1);
+ const l3 = 4.5 * u * (u - 1/3) * (u - 2/3);
+ const A = l1, B = l2;
+ const target = values[i] - l0 * v0 - l3 * v3;
+ sA2 += A * A; sB2 += B * B; sAB += A * B;
+ sAT += A * target; sBT += B * target;
}
const det = sA2 * sB2 - sAB * sAB;
let v1, v2;
if (Math.abs(det) < 1e-9) {
- // Fallback to simple 1/3, 2/3 heuristic if matrix is singular
const idx1 = Math.round(n / 3);
const idx2 = Math.round(2 * n / 3);
v1 = values[idx1];
generated by cgit v1.2.3 (git 2.51.0) at 2026-02-18 20:54:40 +0000