1// In-place radix-2 complex FFT and a Hann window. The twiddle factors and
2// the bit-reversal permutation are cached per transform length, so calling
3// fft repeatedly at one length (as the spectrogram does) costs only the
4// butterflies.
6interface Tables {
7 cos: Float32Array;
8 sin: Float32Array;
9 rev: Uint32Array;
10}
12const cache = new Map<number, Tables>();
14function tables(n: number): Tables {
15 let t = cache.get(n);
16 if (t) return t;
17 const half = n >> 1;
18 const cos = new Float32Array(half);
19 const sin = new Float32Array(half);
20 for (let k = 0; k < half; k++) {
21 const a = (-2 * Math.PI * k) / n;
22 cos[k] = Math.cos(a);
23 sin[k] = Math.sin(a);
24 }
25 const rev = new Uint32Array(n);
26 const bits = 31 - Math.clz32(n);
27 for (let i = 0; i < n; i++) {
28 let r = 0;
29 for (let b = 0; b < bits; b++) r |= ((i >> b) & 1) << (bits - 1 - b);
30 rev[i] = r;
31 }
32 t = { cos, sin, rev };
33 cache.set(n, t);
34 return t;
35}
37/**
38 * Forward complex FFT of (re, im) in place, X[k] = sum_n x[n] exp(-2 pi i k n / N).
39 * The length must be a power of two and the two arrays must have equal length.
40 */
41export function fft(re: Float32Array, im: Float32Array): void {
42 const n = re.length;
43 if (n < 1 || (n & (n - 1)) !== 0) throw new Error(`fft: length ${n} is not a power of two`);
44 if (im.length !== n) throw new Error(`fft: re (${n}) and im (${im.length}) differ in length`);
45 if (n === 1) return;
46 const { cos, sin, rev } = tables(n);
47 for (let i = 0; i < n; i++) {
48 const r = rev[i];
49 if (r > i) {
50 const tr = re[i];
51 re[i] = re[r];
52 re[r] = tr;
53 const ti = im[i];
54 im[i] = im[r];
55 im[r] = ti;
56 }
57 }
58 for (let size = 2; size <= n; size <<= 1) {
59 const half = size >> 1;
60 const step = n / size;
61 for (let start = 0; start < n; start += size) {
62 for (let k = 0, w = 0; k < half; k++, w += step) {
63 const c = cos[w];
64 const s = sin[w];
65 const a = start + k;
66 const b = a + half;
67 const xr = re[b] * c - im[b] * s;
68 const xi = re[b] * s + im[b] * c;
69 re[b] = re[a] - xr;
70 im[b] = im[a] - xi;
71 re[a] += xr;
72 im[a] += xi;
73 }
74 }
75 }
76}
78/** Periodic Hann window of length n, w[k] = 0.5 (1 - cos(2 pi k / n)). */
79export function hann(n: number): Float32Array {
80 const w = new Float32Array(n);
81 for (let k = 0; k < n; k++) w[k] = 0.5 - 0.5 * Math.cos((2 * Math.PI * k) / n);
82 return w;
83}