fft

.gitignore

@@ -44,3 +44,5 @@ dist
 
 .DS_Store
 .eslintcache
+
+CPU.*

benchmark/reimFFT.ts

@@ -0,0 +1,28 @@
+import { reimFFT } from '../src/reim/reimFFT.ts';
+
+const size = 2 ** 16;
+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();
+}
+// Warmup
+reimFFT({ re, im });
+
+// Benchmark: repeat until ~5s
+const targetMs = 5000;
+let iterations = 0;
+
+console.time('reimFFT');
+const start = performance.now();
+while (performance.now() - start < targetMs) {
+  reimFFT({ re, im });
+  iterations++;
+}
+console.timeEnd('reimFFT');
+
+console.log(`${iterations} iterations, ${size} points each`);
+console.log(
+  `${((performance.now() - start) / iterations).toFixed(2)} ms per iteration`,
+);

eslint.config.js

@@ -1,4 +1,7 @@
 import { defineConfig, globalIgnores } from 'eslint/config';
 import cheminfo from 'eslint-config-cheminfo-typescript';
 
-export default defineConfig(globalIgnores(['coverage', 'lib']), cheminfo);
+export default defineConfig(
+  globalIgnores(['coverage', 'lib', 'benchmark']),
+  cheminfo,
+);

src/reim/__tests__/reimFFT.test.ts

@@ -12,4 +12,17 @@ test('reimFFT', () => {
   });
 
   expect(inverse.re).toStrictEqual(re);
+  expect(inverse.im).toStrictEqual(im);
+  // check pointer are different
+  expect(inverse.re).not.toBe(re);
+  expect(inverse.im).not.toBe(im);
+
+  const transformed2 = reimFFT({ re, im }, { inPlace: true });
+  const inverse2 = reimFFT(transformed2, { inverse: true, inPlace: true });
+
+  expect(inverse2.re).toStrictEqual(re);
+  expect(inverse2.im).toStrictEqual(im);
+  // check pointer are the same
+  expect(inverse2.re).toBe(re);
+  expect(inverse2.im).toBe(im);
 });

src/reim/reimFFT.ts

@@ -3,9 +3,13 @@ import FFT from 'fft.js';
 import type { DataReIm } from '../types/index.ts';
 import { xRotate } from '../x/index.ts';
 
+let output: Float64Array | undefined;
+let complexArray: Float64Array | undefined;
+let fft: FFT | undefined;
 export interface ReimFFTOptions {
   inverse?: boolean;
   applyZeroShift?: boolean;
+  inPlace?: boolean;
 }
 
 /**
@@ -18,20 +22,21 @@ export function reimFFT(
   data: DataReIm,
   options: ReimFFTOptions = {},
 ): DataReIm<Float64Array> {
-  const { inverse = false, applyZeroShift = false } = options;
-
   const { re, im } = data;
   const size = re.length;
   const csize = size << 1;
 
-  let complexArray = new Float64Array(csize);
+  const { inPlace = false, inverse = false, applyZeroShift = false } = options;
+
+  if (fft?.size !== size) fft = new FFT(size);
+
+  if (complexArray?.length !== csize) complexArray = new Float64Array(csize);
   for (let i = 0; i < csize; i += 2) {
     complexArray[i] = re[i >>> 1];
     complexArray[i + 1] = im[i >>> 1];
   }
 
-  const fft = new FFT(size);
-  let output = new Float64Array(csize);
+  if (output?.length !== csize) output = new Float64Array(csize);
   if (inverse) {
     if (applyZeroShift) complexArray = zeroShift(complexArray, true);
     fft.inverseTransform(output, complexArray);
@@ -40,14 +45,21 @@ export function reimFFT(
     if (applyZeroShift) output = zeroShift(output);
   }
 
-  const newRe = new Float64Array(size);
-  const newIm = new Float64Array(size);
-  for (let i = 0; i < csize; i += 2) {
-    newRe[i >>> 1] = output[i];
-    newIm[i >>> 1] = output[i + 1];
+  if (inPlace) {
+    for (let i = 0; i < csize; i += 2) {
+      re[i >>> 1] = output[i];
+      im[i >>> 1] = output[i + 1];
+    }
+    return data as DataReIm<Float64Array>;
+  } else {
+    const newRe = new Float64Array(size);
+    const newIm = new Float64Array(size);
+    for (let i = 0; i < csize; i += 2) {
+      newRe[i >>> 1] = output[i];
+      newIm[i >>> 1] = output[i + 1];
+    }
+    return { re: newRe, im: newIm };
   }
-
-  return { re: newRe, im: newIm };
 }
 
 function zeroShift(