improve algorithm

src/VowelChart.tsx

@@ -21,9 +21,6 @@ export function VowelChart(props: {
         mousePath: [],
     }
 
-    const [width, setWidth] = Solid.createSignal(700)
-    const [height, setHeight] = Solid.createSignal(400)
-
 
     Solid.onMount(() => {
         state.ctx = state.canvas.getContext("2d")!
@@ -150,7 +147,9 @@ function mouseMove(
     state: State,
     ev: MouseEvent | TouchEvent)
 {
-    ev.preventDefault()
+    if (state.mouseDown)
+        ev.preventDefault()
+
     updateMousePos(state, ev)
 }
 
@@ -159,7 +158,9 @@ function mouseUp(
     state: State,
     ev: MouseEvent | TouchEvent)
 {
-    ev.preventDefault()
+    if (state.mouseDown)
+        ev.preventDefault()
+
     state.mouseDown = false
     state.synth.setGain(0)
     updateMousePos(state, ev)

src/formantExtractor.ts

@@ -1,27 +1,31 @@
-import { forwardLinearPrediction } from "./lpc.ts"
+import { forwardLinearPrediction, praatBurgMethod } from "./lpc.ts"
 import { findRoots, Complex } from "./roots.ts"
 
 
 // From: https://www.mathworks.com/help/signal/ug/formant-estimation-with-lpc-coefficients.html
 // From: https://github.com/praat/praat/blob/master/fon/Sound_to_Formant.cpp
 export function extractFormants(
-    /// An array of samples in the range [-1, 1]
+    /// An array of sound samples in the range [-1, 1]
     sample: Float32Array,
     /// The sampling frequency in hertz
     samplingFrequency: number)
     /// Returns the frequencies of the formants in hertz
     : number[]
 {
-    const sampleWindowed =
-        //sample.map((s, i) => s * hammingWindow(i, sample.length))
-        sample.map((s, i) => s * praatGaussianWindow(i, sample.length))
+    if (sample.every(s => s === 0))
+        return []
 
-    const sampleFiltered =
+    const samplePreemphasized =
         //preemphasisFilter(sampleWindowed)
-        praatPreemphasis(sampleWindowed, samplingFrequency)
+        praatPreemphasis(sample, samplingFrequency)
+
+    const sampleWindowed =
+        //sample.map((s, i) => s * hammingWindow(i, sample.length))
+        samplePreemphasized.map((s, i) => s * praatGaussianWindow(i, sample.length))
 
     const lpc =
-        forwardLinearPrediction(sampleFiltered, 10)
+        //forwardLinearPrediction(sampleWindowed, 10)
+        praatBurgMethod(sampleWindowed, 10)
 
     const roots = findRoots(lpc)
         .filter(c => c.imag >= 0)
@@ -30,11 +34,13 @@ export function extractFormants(
     const angles = roots
         .map(c => Math.atan2(c.imag, c.real))
 
+    const nyquistFrequency = samplingFrequency / 2
+
     const frequencies = angles
-        .map(a => a * samplingFrequency / 2 / Math.PI)
+        .map(a => a * nyquistFrequency / Math.PI)
 
     const bandwidths = roots
-        .map(r => -Math.log(complexMagnitude(r)) * samplingFrequency / 2 / Math.PI)
+        .map(r => -Math.log(complexMagnitude(r)) * nyquistFrequency / Math.PI)
 
     const formants = []
     for (let i = 0; i < angles.length; i++)
@@ -118,11 +124,11 @@ function praatPreemphasis(
 {
     const result = new Float32Array(array)
 
-    const frequency = 50
+    const preemphasisFrequency = 50
     const dx = 1 / samplingFrequency
-    const preEmphasis = Math.exp(-2.0 * Math.PI * frequency * dx)
+    const preEmphasis = Math.exp(-2.0 * Math.PI * preemphasisFrequency * dx)
 
-    for (let i = array.length - 1; i >= 2; i--)
+    for (let i = array.length - 1; i >= 1; i--)
         result[i] -= preEmphasis * result[i - 1]
 
     return result

src/lpc.ts

@@ -1,3 +1,4 @@
+// Verified equivalent to Praat's Burg method.
 export function forwardLinearPrediction(
     sample: Float32Array,
     numCoeffs: number)
@@ -48,6 +49,101 @@ export function forwardLinearPrediction(
 }
 
 
+// From: https://github.com/praat/praat/blob/master/dwsys/NUM2.cpp#L1370
+export function praatBurgMethod(
+    samples: Float32Array,
+    numCoeffs: number)
+    : number[]
+{
+    const a = new Array<number>(numCoeffs).fill(0)
+
+	if (samples.length <= 2)
+    {
+		a[1] = -1.0
+        return a
+		//return ( n == 2 ? 0.5 * (x [1] * x [1] + x [2] * x [2]) : x [1] * x [1] );
+	}
+
+    const b1 = new Array<number>(samples.length).fill(0)
+    const b2 = new Array<number>(samples.length).fill(0)
+    const aa = new Array<number>(numCoeffs).fill(0)
+
+    const n = samples.length
+    const m = numCoeffs
+
+	// (3)
+
+	let p = 0
+	for (let j = 1; j <= n; j++)
+		p += samples[j - 1] * samples[j - 1]
+
+	let xms = p / n
+	if (xms <= 0)
+		return a // warning empty
+
+	// (9)
+
+	b1[1 - 1] = samples[1 - 1]
+	if (n < 2)
+		return a
+
+	b2[n - 1 - 1] = samples[n - 1]
+	for (let j = 2; j <= n - 1; j++)
+    {
+		b1[j - 1] = samples[j - 1]
+        b2[j - 1 - 1] = samples[j - 1]
+    }
+
+	for (let i = 1; i <= m; i++)
+    {
+		// (7)
+
+		let num = 0
+        let denum = 0
+		for (let j = 1; j <= n - i; j++)
+        {
+			num += b1[j - 1] * b2[j - 1]
+			denum += b1[j - 1] * b1[j - 1] + b2[j - 1] * b2[j - 1]
+		}
+
+		if (denum <= 0)
+			return a // warning ill-conditioned
+
+		a[i - 1] = 2 * num / denum
+
+		// (10)
+
+		xms *= 1 - a[i - 1] * a[i - 1]
+
+		// (5)
+
+		for (let j = 1; j <= i - 1; j++)
+			a[j - 1] = aa[j - 1] - a[i - 1] * aa[i - j - 1]
+
+		if (i < m)
+        {
+			// (8) Watch out: i -> i+1
+
+			for (let j = 1; j <= i; j++)
+				aa[j - 1] = a[j - 1]
+
+			for (let j = 1; j <= n - i - 1; j++)
+            {
+				b1[j - 1] -= aa[i - 1] * b2[j - 1]
+				b2[j - 1] = b2[j + 1 - 1] - aa[i - 1] * b1[j + 1 - 1]
+			}
+		}
+	}
+
+    for (let i = 0; i < a.length; i++)
+        a[i] *= -1
+
+    a.unshift(1)
+
+	return a
+}
+
+
 /*void ForwardLinearPrediction( vector<double> &coeffs, const vector<double> &x )
 {
     // GET SIZE FROM INPUT VECTORS

src/main.tsx

@@ -7,6 +7,8 @@ import { VowelChart } from "./VowelChart.tsx"
 import { AnalysisChart } from "./AnalysisChart.tsx"
 import { RecordingPanel } from "./RecordingPanel.tsx"
 
+//import "./test.ts"
+
 
 function Page()
 {