Harden

Author: dy

.eslintrc.json

@@ -1,20 +1,110 @@
-import rfft from '../index.js'
+import { createRequire } from 'node:module'
+import { fft } from '../index.js'
 
-const sizes = [256, 1024, 4096, 16384]
-const warmup = 1000
-const iterations = 50000
+const require = createRequire(import.meta.url)
 
-for (const N of sizes) {
-	const input = new Float32Array(N)
-	for (let i = 0; i < N; i++) input[i] = Math.random() * 2 - 1
+const N = 4096
+const WARMUP = 500
+const ITERATIONS = 20000
 
-	// warmup (also primes twiddle cache)
-	for (let i = 0; i < warmup; i++) rfft(input)
+// Shared input
+const input = new Float32Array(N)
+for (let i = 0; i < N; i++) input[i] = Math.random() * 2 - 1
+
+// --- Setup each library outside the timed loop ---
+
+// fourier-transform (this package)
+const ftRun = () => fft(input)
+
+// fft.js (indutny, radix-4) — widely considered fastest pure-JS FFT
+const FFTJS = require('fft.js')
+const fftjs = new FFTJS(N)
+const fftjsOut = fftjs.createComplexArray()
+const fftjsRun = () => fftjs.realTransform(fftjsOut, input)
+
+// ml-fft (in-place, pre-allocated arrays)
+const { FFT: MLFFT } = require('ml-fft')
+MLFFT.init(N)
+const mlRe = new Array(N), mlIm = new Array(N)
+const mlfftRun = () => {
+	for (let i = 0; i < N; i++) { mlRe[i] = input[i]; mlIm[i] = 0 }
+	MLFFT.fft(mlRe, mlIm)
+}
+
+// dsp.js (our ancestor, the original split-radix)
+const dsp = require('dsp.js')
+const dspfft = new dsp.FFT(N, 44100)
+const dspRun = () => dspfft.forward(input)
+
+// ndarray-fft
+const ndfft = require('ndarray-fft')
+const ndarray = require('ndarray')
+const ndRe = new Float64Array(N), ndIm = new Float64Array(N)
+const ndx = ndarray(ndRe), ndy = ndarray(ndIm)
+const ndfftRun = () => {
+	for (let i = 0; i < N; i++) { ndRe[i] = input[i]; ndIm[i] = 0 }
+	ndfft(1, ndx, ndy)
+}
+
+// ooura (port of Ooura's C FFT, radix-4)
+const Ooura = require('ooura')
+const oo = new Ooura(N, { type: 'real', radix: 4 })
+const ooIn = new Float64Array(N)
+const ooRe = oo.scalarArrayFactory()
+const ooIm = oo.scalarArrayFactory()
+const oouraRun = () => {
+	for (let i = 0; i < N; i++) ooIn[i] = input[i]
+	oo.fft(ooIn.buffer, ooRe.buffer, ooIm.buffer)
+}
+
+// kissfft-wasm (WASM KissFFT)
+const { rfft: kissRfft } = await import('kissfft-wasm')
+const kissInput = new Float64Array(input)
+const kissRun = () => kissRfft(kissInput)
+
+// fft-js (naive Cooley-Tukey, educational)
+const fftjs2 = require('fft-js')
+const fftjs2Input = Array.from(input)
+const fftjs2Run = () => fftjs2.fft(fftjs2Input)
+
+// --- Benchmark runner ---
+
+const libs = [
+	['fourier-transform', ftRun],
+	['fft.js (indutny)', fftjsRun],
+	['ml-fft', mlfftRun],
+	['dsp.js', dspRun],
+	['ndarray-fft', ndfftRun],
+	['ooura', oouraRun],
+	['kissfft-wasm', kissRun],
+	['fft-js', fftjs2Run],
+]
+
+console.log(`\nFFT benchmark — N=${N}, ${ITERATIONS} iterations\n`)
+
+const results = []
+
+for (const [name, fn] of libs) {
+	// warmup
+	for (let i = 0; i < WARMUP; i++) fn()
 
 	const start = performance.now()
-	for (let i = 0; i < iterations; i++) rfft(input)
+	for (let i = 0; i < ITERATIONS; i++) fn()
 	const elapsed = performance.now() - start
 
-	const us = (elapsed / iterations * 1000).toFixed(1)
-	console.log(`N=${String(N).padStart(5)}  ${us.padStart(6)}µs/call  (${iterations} iters, ${elapsed.toFixed(0)}ms)`)
+	const us = elapsed / ITERATIONS * 1000
+	results.push({ name, us, elapsed })
+}
+
+// Sort by speed
+results.sort((a, b) => a.us - b.us)
+
+const fastest = results[0].us
+for (const { name, us, elapsed } of results) {
+	const ratio = us / fastest
+	const bar = '█'.repeat(Math.min(40, Math.round(ratio * 4)))
+	console.log(
+		`${name.padEnd(22)} ${us.toFixed(1).padStart(7)}µs/call  ${ratio === 1 ? '      ' : ('×' + ratio.toFixed(1)).padStart(6)}  ${bar}`
+	)
 }
+console.log()

.travis.yml

@@ -1,11 +1,7 @@
 /**
  * Real-valued split-radix FFT.
- * Returns magnitude spectrum as Float64Array of length N/2.
  *
  * @module fourier-transform
- * @param {Float32Array|Float64Array|Array} input - Time-domain signal, length must be power of 2 (>=2).
- * @param {Float64Array} [output] - Optional pre-allocated output buffer (length N/2). If omitted, returns internal buffer (overwritten on next call with same N).
- * @returns {Float64Array} Magnitude spectrum, length N/2.
  */
 
 const { sqrt, sin, cos, abs, SQRT1_2 } = Math
@@ -16,12 +12,16 @@ const cache = new Map()
 let lastN = 0, lastEntry = null
 
 function init(N) {
+	const half = N >>> 1
 	const x = new Float64Array(N)
-	const spectrum = new Float64Array(N >>> 1)
+	const spectrum = new Float64Array(half)
+	const re = new Float64Array(half + 1)
+	const im = new Float64Array(half + 1)
+	const complex = { re, im }
 	const bSi = 2 / N
 
 	// Count twiddle factors per stage
-	let total = 0, n2 = 2, nn = N >>> 1
+	let total = 0, n2 = 2, nn = half
 	const stages = []
 	while ((nn = nn >>> 1)) {
 		n2 = n2 << 1
@@ -33,7 +33,7 @@ function init(N) {
 
 	// Interleaved twiddle table: [cc1, ss1, cc3, ss3] per entry
 	const tw = new Float64Array(total << 2)
-	n2 = 2; nn = N >>> 1
+	n2 = 2; nn = half
 	let si = 0
 	while ((nn = nn >>> 1)) {
 		n2 = n2 << 1
@@ -52,17 +52,26 @@ function init(N) {
 		si++
 	}
 
-	const entry = { x, spectrum, bSi, tw, stages }
+	const entry = { x, spectrum, complex, bSi, tw, stages }
 	cache.set(N, entry)
 	return entry
 }
 
-export default function rfft(input, output) {
+function getEntry(N) {
+	if (N === lastN) return lastEntry
+	const entry = cache.get(N) || init(N)
+	lastN = N
+	lastEntry = entry
+	return entry
+}
+
+// Shared butterfly computation
+function transform(input) {
 	const N = input.length
 	if (N < 2 || (N & (N - 1))) throw Error('Input length must be a power of 2 (>= 2).')
 
-	const entry = N === lastN ? lastEntry : (lastEntry = cache.get(N) || init(N), lastN = N, lastEntry)
-	const { x, spectrum, bSi, tw, stages } = entry
+	const entry = getEntry(N)
+	const { x, tw, stages } = entry
 
 	reverseBinPermute(N, x, input)
 
@@ -171,8 +180,21 @@ export default function rfft(input, output) {
 		si++
 	}
 
-	// Magnitude spectrum: Re(X[k]) = x[k], Im(X[k]) = x[N-k]
+	return entry
+}
+
+/**
+ * Compute magnitude spectrum of real-valued input.
+ * @param {ArrayLike<number>} input - length must be power of 2 (>= 2).
+ * @param {Float64Array} [output] - Optional buffer (length N/2). If omitted, returns internal view (overwritten on next call with same N).
+ * @returns {Float64Array} Magnitude spectrum, length N/2.
+ */
+export default function rfft(input, output) {
+	const entry = transform(input)
+	const N = input.length
+	const { x, spectrum, bSi } = entry
 	const out = output || spectrum
+
 	let i = N >>> 1
 	while (--i) {
 		const rval = x[i], ival = x[N - i]
@@ -183,6 +205,32 @@ export default function rfft(input, output) {
 	return out
 }
 
+/**
+ * Compute complex spectrum of real-valued input (unnormalized DFT).
+ * @param {ArrayLike<number>} input - length must be power of 2 (>= 2).
+ * @param {{re: Float64Array, im: Float64Array}} [output] - Optional buffers (length N/2+1 each). If omitted, returns internal view.
+ * @returns {{re: Float64Array, im: Float64Array}} Complex spectrum, N/2+1 bins (DC through Nyquist).
+ */
+export function fft(input, output) {
+	const entry = transform(input)
+	const N = input.length
+	const half = N >>> 1
+	const { x, complex } = entry
+	const re = output ? output.re : complex.re
+	const im = output ? output.im : complex.im
+
+	re[0] = x[0]
+	im[0] = 0
+	for (let k = 1; k < half; k++) {
+		re[k] = x[k]
+		im[k] = x[N - k]
+	}
+	re[half] = x[half]
+	im[half] = 0
+
+	return output || complex
+}
+
 function reverseBinPermute(N, dest, source) {
 	const halfSize = N >>> 1
 	const nm1 = N - 1