1/**
2 * In-place iterative radix-2 complex FFT on separate real and imaginary
3 * arrays. Sizes must be powers of two. The inverse transform is scaled by 1/n.
4 */
5export class FFT {
6 readonly n: number
7 private rev: Uint32Array
8 private cos: Float64Array
9 private sin: Float64Array
11 constructor(n: number) {
12 if (n < 2 || (n & (n - 1)) !== 0) throw new Error(`FFT size must be a power of two, got ${n}`)
13 this.n = n
14 const bits = Math.round(Math.log2(n))
15 this.rev = new Uint32Array(n)
16 for (let i = 0; i < n; i++) {
17 let r = 0
18 let x = i
19 for (let b = 0; b < bits; b++) {
20 r = (r << 1) | (x & 1)
21 x >>= 1
22 }
23 this.rev[i] = r >>> 0
24 }
25 this.cos = new Float64Array(n / 2)
26 this.sin = new Float64Array(n / 2)
27 for (let i = 0; i < n / 2; i++) {
28 const a = (-2 * Math.PI * i) / n
29 this.cos[i] = Math.cos(a)
30 this.sin[i] = Math.sin(a)
31 }
32 }
34 transform(re: Float64Array, im: Float64Array, inverse = false): void {
35 const n = this.n
36 const rev = this.rev
37 for (let i = 0; i < n; i++) {
38 const j = rev[i]
39 if (j > i) {
40 const tr = re[i]
41 re[i] = re[j]
42 re[j] = tr
43 const ti = im[i]
44 im[i] = im[j]
45 im[j] = ti
46 }
47 }
48 const sgn = inverse ? -1 : 1
49 for (let size = 2; size <= n; size <<= 1) {
50 const half = size >> 1
51 const step = n / size
52 for (let start = 0; start < n; start += size) {
53 for (let k = 0, tw = 0; k < half; k++, tw += step) {
54 const wr = this.cos[tw]
55 const wi = sgn * this.sin[tw]
56 const i = start + k
57 const j = i + half
58 const tr = re[j] * wr - im[j] * wi
59 const ti = re[j] * wi + im[j] * wr
60 re[j] = re[i] - tr
61 im[j] = im[i] - ti
62 re[i] += tr
63 im[i] += ti
64 }
65 }
66 }
67 if (inverse) {
68 for (let i = 0; i < n; i++) {
69 re[i] /= n
70 im[i] /= n
71 }
72 }
73 }
74}
76export function nextPow2(n: number): number {
77 let p = 1
78 while (p < n) p <<= 1
79 return p
80}
82const cache = new Map<number, FFT>()
84/** Shared FFT plans, one per size. */
85export function fftPlan(n: number): FFT {
86 let p = cache.get(n)
87 if (!p) {
88 p = new FFT(n)
89 cache.set(n, p)
90 }
91 return p
92}