1 files changed, 20 insertions, 25 deletions

diff --git a/tools/mq_editor/mq_extract.js b/tools/mq_editor/mq_extract.js
index ff38d63..18897fb 100644
--- a/tools/mq_editor/mq_extract.js
+++ b/tools/mq_editor/mq_extract.js
@@ -526,7 +526,10 @@ function trackIsoContour(stftCache, seedTime, seedFreq, params) {
   };
 }
 
-// Fit cubic bezier to trajectory using least-squares for inner control points
+// 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];
@@ -534,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];