perf(fft): cache FFT instances

benchmark/reimFFT.ts

@@ -1,99 +1,95 @@
 /* eslint-disable no-console */
+import FFT from 'fft.js';
+import { XSadd } from 'ml-xsadd';
+
 import { reimFFT } from '../src/reim/reimFFT.ts';
 import { reimArrayFFT } from '../src/reimArray/reimArrayFFT.ts';
+import type { DataReIm } from '../src/types/index.ts';
 
-const size = 2 ** 16;
-const count = 10; // number of spectra in the array benchmark
+const size = 2 ** 16; // 64k-point transform: FFT setup dominates the cost
+const count = 10; // number of spectra processed per round
+const targetMs = 5000;
 
-// Build input data
+// Deterministic, reproducible input so every section runs on identical data.
+const { random } = new XSadd(42);
 const spectra = Array.from({ length: count }, () => {
   const re = new Float64Array(size);
   const im = new Float64Array(size);
   for (let i = 0; i < size; i++) {
-    re[i] = Math.random();
-    im[i] = Math.random();
+    re[i] = random();
+    im[i] = random();
   }
   return { re, im };
 });
 
-// Warmup
-for (const s of spectra) reimFFT(s);
-for (const s of spectra) reimFFT(s, { inPlace: true });
-reimArrayFFT(spectra);
-reimArrayFFT(spectra, { inPlace: true });
-
-const targetMs = 5000;
+/**
+ * `reimFFT` as it was *before* the cache fix: a fresh `FFT` instance is built on
+ * every call. Kept here as the baseline to confirm the cached version is faster.
+ * @param data - complex spectrum.
+ * @returns FFT of the complex spectrum.
+ */
+function reimFFTNoCache(data: DataReIm): DataReIm<Float64Array> {
+  const { re, im } = data;
+  const length = re.length;
+  const csize = length << 1;
 
-// --- reimFFT (loop over each spectrum individually) ---
-{
-  let iterations = 0;
-  const start = performance.now();
-  console.time('reimFFT (loop)');
-  while (performance.now() - start < targetMs) {
-    for (const s of spectra) reimFFT(s);
-    iterations++;
+  const complexArray = new Float64Array(csize);
+  for (let i = 0; i < csize; i += 2) {
+    complexArray[i] = re[i >>> 1];
+    complexArray[i + 1] = im[i >>> 1];
   }
-  const elapsed = performance.now() - start;
-  console.timeEnd('reimFFT (loop)');
-  console.log(
-    `  ${iterations * count} total FFTs, ${count} spectra × ${iterations} rounds`,
-  );
-  console.log(`  ${(elapsed / (iterations * count)).toFixed(3)} ms per FFT`);
-}
 
-console.log('');
+  const fft = new FFT(length);
+  const output = new Float64Array(csize);
+  fft.transform(output, complexArray);
 
-// --- reimFFT inPlace (loop over each spectrum individually) ---
-{
-  let iterations = 0;
-  const start = performance.now();
-  console.time('reimFFT inPlace (loop)');
-  while (performance.now() - start < targetMs) {
-    for (const s of spectra) reimFFT(s, { inPlace: true });
-    iterations++;
+  const newRe = new Float64Array(length);
+  const newIm = new Float64Array(length);
+  for (let i = 0; i < csize; i += 2) {
+    newRe[i >>> 1] = output[i];
+    newIm[i >>> 1] = output[i + 1];
   }
-  const elapsed = performance.now() - start;
-  console.timeEnd('reimFFT inPlace (loop)');
-  console.log(
-    `  ${iterations * count} total FFTs, ${count} spectra × ${iterations} rounds`,
-  );
-  console.log(`  ${(elapsed / (iterations * count)).toFixed(3)} ms per FFT`);
+  return { re: newRe, im: newIm };
 }
 
-console.log('');
-
-// --- reimArrayFFT (single call for the whole array) ---
-{
-  let iterations = 0;
+/**
+ * Run `task` repeatedly for `targetMs` and report the time per FFT. Each round
+ * performs `count` transforms.
+ * @param label - section name.
+ * @param task - one round of work (transforms all `count` spectra).
+ */
+function bench(label: string, task: () => void): void {
+  task(); // warmup
+  let rounds = 0;
   const start = performance.now();
-  console.time('reimArrayFFT');
   while (performance.now() - start < targetMs) {
-    reimArrayFFT(spectra);
-    iterations++;
+    task();
+    rounds++;
   }
   const elapsed = performance.now() - start;
-  console.timeEnd('reimArrayFFT');
+  const totalFFTs = rounds * count;
+  console.log(label);
   console.log(
-    `  ${iterations * count} total FFTs, ${count} spectra × ${iterations} rounds`,
+    `  ${(elapsed / totalFFTs).toFixed(3)} ms per FFT  (${totalFFTs} FFTs over ${rounds} rounds)`,
   );
-  console.log(`  ${(elapsed / (iterations * count)).toFixed(3)} ms per FFT`);
+  console.log('');
 }
 
+console.log(`FFT size: ${size} (2^16), ${count} spectra per round`);
 console.log('');
 
-// --- reimArrayFFT inPlace (single call for the whole array) ---
-{
-  let iterations = 0;
-  const start = performance.now();
-  console.time('reimArrayFFT inPlace');
-  while (performance.now() - start < targetMs) {
-    reimArrayFFT(spectra, { inPlace: true });
-    iterations++;
-  }
-  const elapsed = performance.now() - start;
-  console.timeEnd('reimArrayFFT inPlace');
-  console.log(
-    `  ${iterations * count} total FFTs, ${count} spectra × ${iterations} rounds`,
-  );
-  console.log(`  ${(elapsed / (iterations * count)).toFixed(3)} ms per FFT`);
-}
+// Before the fix: new FFT instance per call.
+bench('reimFFT — before fix (new FFT per call)', () => {
+  for (const spectrum of spectra) reimFFTNoCache(spectrum);
+});
+
+// After the fix: FFT instance cached per size and reused across calls.
+bench('reimFFT — after fix (cached FFT instance)', () => {
+  for (const spectrum of spectra) reimFFT(spectrum);
+});
+
+// reimArrayFFT reuses a single FFT instance (and working buffers) for the whole
+// array in one call.
+bench('reimArrayFFT (single shared FFT instance)', () => {
+  reimArrayFFT(spectra);
+});

benchmark/xHilbertTransform.ts

@@ -0,0 +1,94 @@
+/* eslint-disable no-console */
+import FFT from 'fft.js';
+import { XSadd } from 'ml-xsadd';
+
+import { xHilbertTransform } from '../src/x/xHilbertTransform.ts';
+
+const size = 2 ** 16; // power of two => the FFT path (hilbertTransformWithFFT)
+const count = 10; // signals processed per round
+const targetMs = 5000;
+
+// Deterministic, reproducible input.
+const { random } = new XSadd(42);
+const signals = Array.from({ length: count }, () => {
+  const array = new Float64Array(size);
+  for (let i = 0; i < size; i++) array[i] = random() * 2 - 1;
+  return array;
+});
+
+/**
+ * Hilbert transform via FFT as it was *before* the shared cache: a fresh `FFT`
+ * instance is built on every call. Mirrors `hilbertTransformWithFFT`, used as
+ * the baseline to confirm the cached version is faster.
+ * @param array - real input signal whose length is a power of two.
+ * @returns the Hilbert transform (90° phase-shifted signal).
+ */
+function hilbertNoCache(array: Float64Array): Float64Array {
+  const length = array.length;
+  const fft = new FFT(length);
+
+  const spectrum = new Float64Array(length * 2);
+  fft.realTransform(spectrum, array);
+  fft.completeSpectrum(spectrum);
+
+  const half = length >> 1;
+  const nyquist = half << 1;
+  spectrum[nyquist] = 0;
+  spectrum[nyquist + 1] = 0;
+  for (let j = (half + 1) << 1; j < spectrum.length; j += 2) {
+    spectrum[j] = -spectrum[j];
+    spectrum[j + 1] = -spectrum[j + 1];
+  }
+
+  const hilbertSignal = new Float64Array(length * 2);
+  fft.inverseTransform(hilbertSignal, spectrum);
+
+  const result = new Float64Array(length);
+  for (let i = 0; i < length; i++) result[i] = hilbertSignal[i * 2 + 1];
+  return result;
+}
+
+/**
+ * Run `task` for `targetMs` and report the time per transform. Each round
+ * processes `count` signals.
+ * @param label - section name.
+ * @param task - one round of work (transforms all `count` signals).
+ */
+function bench(label: string, task: () => void): void {
+  task(); // warmup
+  let rounds = 0;
+  const start = performance.now();
+  while (performance.now() - start < targetMs) {
+    task();
+    rounds++;
+  }
+  const elapsed = performance.now() - start;
+  const total = rounds * count;
+  console.log(label);
+  console.log(
+    `  ${(elapsed / total).toFixed(3)} ms per transform  (${total} transforms over ${rounds} rounds)`,
+  );
+  console.log('');
+}
+
+// Sanity check: the cached path and the baseline must compute the same thing.
+const reference = hilbertNoCache(signals[0]);
+const cached = xHilbertTransform(signals[0]);
+let maxDiff = 0;
+for (let i = 0; i < reference.length; i++) {
+  const diff = Math.abs(reference[i] - cached[i]);
+  if (diff > maxDiff) maxDiff = diff;
+}
+
+console.log(
+  `xHilbertTransform: size ${size} (2^16), ${count} signals per round`,
+);
+console.log(`equivalence check: max abs diff ${maxDiff.toExponential(2)}\n`);
+
+bench('before shared cache (new FFT per call)', () => {
+  for (const signal of signals) hilbertNoCache(signal);
+});
+
+bench('after shared cache (reused FFT instance)', () => {
+  for (const signal of signals) xHilbertTransform(signal);
+});

src/__tests__/__snapshots__/index.test.ts.snap

@@ -181,6 +181,8 @@ exports[`existence of exported functions 1`] = `
   "matrixZRescale",
   "matrixZRescalePerColumn",
   "matrixTranspose",
+  "clearFFTCache",
+  "setFFTCacheMaxSize",
   "createNumberArray",
   "createDoubleArray",
   "createFromToArray",

src/matrix/matrixHilbertTransform.ts

@@ -1,5 +1,4 @@
-import FFT from 'fft.js';
-
+import { getFFT } from '../utils/fftCache.ts';
 import { isPowerOfTwo } from '../utils/index.ts';
 
 import { matrixCreateEmpty } from './matrixCreateEmpty.ts';
@@ -41,7 +40,7 @@ export function matrixHilbertTransform(
   }
 
   // Single FFT instance reused across all rows
-  const fft = new FFT(size);
+  const fft = getFFT(size);
 
   // Multiplier computed once — identical for every row of the same length
   const half = size >> 1;

src/reim/reimFFT.ts

@@ -1,6 +1,5 @@
-import FFT from 'fft.js';
-
 import type { DataReIm } from '../types/index.ts';
+import { getFFT } from '../utils/fftCache.ts';
 
 import { zeroShift } from './zeroShift.ts';
 
@@ -36,7 +35,7 @@ export function reimFFT(
     complexArray[i + 1] = im[i >>> 1];
   }
 
-  const fft = new FFT(size);
+  const fft = getFFT(size);
   let output = new Float64Array(csize);
   if (inverse) {
     if (applyZeroShift) complexArray = zeroShift(complexArray, true);

src/reimArray/reimArrayFFT.ts

@@ -1,7 +1,6 @@
-import FFT from 'fft.js';
-
 import { zeroShift } from '../reim/zeroShift.ts';
 import type { DataReIm } from '../types/index.ts';
+import { getFFT } from '../utils/fftCache.ts';
 
 export interface ReimArrayFFTOptions {
   inverse?: boolean;
@@ -41,7 +40,7 @@ export function reimArrayFFT(
   }
 
   // Single FFT instance and working buffers reused across all spectra
-  const fft = new FFT(size);
+  const fft = getFFT(size);
   const complexArray = new Float64Array(csize);
   const output = new Float64Array(csize);
 

src/reimMatrix/reimMatrixFFT.ts

@@ -1,7 +1,6 @@
-import FFT from 'fft.js';
-
 import { zeroShift } from '../reim/zeroShift.ts';
 import type { DataReImMatrix } from '../types/index.ts';
+import { getFFT } from '../utils/fftCache.ts';
 
 export interface ReimMatrixFFTOptions {
   inverse?: boolean;
@@ -44,7 +43,7 @@ export function reimMatrixFFT(
   }
 
   // Single FFT instance and working buffers reused across all rows
-  const fft = new FFT(size);
+  const fft = getFFT(size);
   const complexArray = new Float64Array(csize);
   const output = new Float64Array(csize);
 

src/reimMatrix/reimMatrixFFTByColumns.ts

@@ -1,7 +1,6 @@
-import FFT from 'fft.js';
-
 import { zeroShift } from '../reim/zeroShift.ts';
 import type { DataReImMatrix } from '../types/index.ts';
+import { getFFT } from '../utils/fftCache.ts';
 
 export interface ReimMatrixFFTByColumnsOptions {
   inverse?: boolean;
@@ -49,7 +48,7 @@ export function reimMatrixFFTByColumns(
   }
 
   // Single FFT instance and working buffers reused across all columns
-  const fft = new FFT(numRows);
+  const fft = getFFT(numRows);
   const complexArray = new Float64Array(csize);
   const output = new Float64Array(csize);