Initial Commit

packages/media/audio/audio-analysis-beat/src/algorithms/spectral-flux.ts

@@ -0,0 +1,141 @@
+/**
+ * Spectral-flux onset detection function via a hand-rolled radix-2 real FFT.
+ *
+ * Per-frame magnitude spectra are differenced bin-by-bin; only positive
+ * differences count toward the ODF. No browser analyser node, no Node bindings — the
+ * implementation is pure numeric code over `Float32Array`/`Float64Array` so
+ * the module runs identically in Bun, Node.js, Deno, and browsers.
+ */
+
+/**
+ * In-place Cooley-Tukey radix-2 complex FFT.
+ *
+ * `re` and `im` carry the real and imaginary parts; both arrays must have the
+ * same length and that length must be a power of two. The transform is in
+ * place — on return `re[k]` and `im[k]` carry the `k`th complex bin.
+ */
+export function fftRadix2(re: Float64Array, im: Float64Array): void {
+    const n = re.length;
+    if (n <= 1) return;
+    // Bit-reversal permutation
+    let j = 0;
+    for (let i = 1; i < n; i++) {
+        let bit = n >> 1;
+        for (; (j & bit) !== 0; bit >>= 1) {
+            j ^= bit;
+        }
+        j ^= bit;
+        if (i < j) {
+            const tr = re[i]!;
+            const ti = im[i]!;
+            re[i] = re[j]!;
+            im[i] = im[j]!;
+            re[j] = tr;
+            im[j] = ti;
+        }
+    }
+    // Butterflies
+    for (let size = 2; size <= n; size <<= 1) {
+        const half = size >> 1;
+        const theta = (-2 * Math.PI) / size;
+        const wRe = Math.cos(theta);
+        const wIm = Math.sin(theta);
+        for (let start = 0; start < n; start += size) {
+            let curRe = 1;
+            let curIm = 0;
+            for (let k = 0; k < half; k++) {
+                const a = start + k;
+                const b = a + half;
+                const tre = re[b]! * curRe - im[b]! * curIm;
+                const tim = re[b]! * curIm + im[b]! * curRe;
+                re[b] = re[a]! - tre;
+                im[b] = im[a]! - tim;
+                re[a] = re[a]! + tre;
+                im[a] = im[a]! + tim;
+                const newRe = curRe * wRe - curIm * wIm;
+                const newIm = curRe * wIm + curIm * wRe;
+                curRe = newRe;
+                curIm = newIm;
+            }
+        }
+    }
+}
+
+/**
+ * Compute the magnitude spectrum of a real-valued frame of length
+ * `fftSize` (which must be a power of two). Returned array has length
+ * `fftSize / 2 + 1` (the non-redundant half plus Nyquist). A Hann window is
+ * applied in-place to a working copy of the frame before transform.
+ */
+export function magnitudeSpectrum(
+    frame: Float32Array,
+    fftSize: number
+): Float64Array {
+    const re = new Float64Array(fftSize);
+    const im = new Float64Array(fftSize);
+    const n = Math.min(frame.length, fftSize);
+    const denom = Math.max(1, fftSize - 1);
+    for (let i = 0; i < n; i++) {
+        const w = 0.5 * (1 - Math.cos((2 * Math.PI * i) / denom));
+        re[i] = frame[i]! * w;
+    }
+    fftRadix2(re, im);
+    const halfPlusOne = (fftSize >> 1) + 1;
+    const mag = new Float64Array(halfPlusOne);
+    for (let k = 0; k < halfPlusOne; k++) {
+        const r = re[k]!;
+        const i = im[k]!;
+        mag[k] = Math.sqrt(r * r + i * i);
+    }
+    return mag;
+}
+
+/**
+ * Single-frame spectral flux: sum of positive bin-to-bin differences between
+ * the current magnitude spectrum and the previous one. `prevMag` is mutated
+ * in place to become the current spectrum so the caller can chain.
+ */
+export function spectralFlux(
+    frame: Float32Array,
+    prevMag: Float64Array,
+    fftSize: number
+): number {
+    const mag = magnitudeSpectrum(frame, fftSize);
+    let flux = 0;
+    for (let k = 0; k < mag.length; k++) {
+        const m = mag[k]!;
+        const d = m - prevMag[k]!;
+        if (d > 0) flux += d;
+        prevMag[k] = m;
+    }
+    return flux;
+}
+
+/**
+ * Build a spectral-flux ODF over a sample stream by splitting it into
+ * `hop`-sample frames of `fftSize` samples each (zero-padded if necessary).
+ * Returns a `Float64Array` of length `ceil(length / hop)`.
+ */
+export function spectralFluxOdf(
+    mono: Float32Array,
+    fftSize: number,
+    hop: number
+): Float64Array {
+    if (mono.length === 0 || fftSize <= 0 || hop <= 0) {
+        return new Float64Array(0);
+    }
+    const halfPlusOne = (fftSize >> 1) + 1;
+    const prev = new Float64Array(halfPlusOne);
+    const frames = Math.ceil(mono.length / hop);
+    const odf = new Float64Array(frames);
+    const buf = new Float32Array(fftSize);
+    for (let f = 0; f < frames; f++) {
+        const start = f * hop;
+        buf.fill(0);
+        const take = Math.min(fftSize, mono.length - start);
+        for (let i = 0; i < take; i++) buf[i] = mono[start + i]!;
+        odf[f] = spectralFlux(buf, prev, fftSize);
+    }
+    odf[0] = 0; // first frame has no predecessor signal
+    return odf;
+}