1// In-place radix-2 complex FFT (length must be a power of two).
2export function fft(re: Float64Array, im: Float64Array): void {
3 const n = re.length;
4 for (let i = 1, j = 0; i < n; i++) {
5 let bit = n >> 1;
6 for (; j & bit; bit >>= 1) j ^= bit;
7 j ^= bit;
8 if (i < j) {
9 [re[i], re[j]] = [re[j], re[i]];
10 [im[i], im[j]] = [im[j], im[i]];
11 }
12 }
13 for (let len = 2; len <= n; len <<= 1) {
14 const ang = (-2 * Math.PI) / len;
15 const wr = Math.cos(ang);
16 const wi = Math.sin(ang);
17 for (let i = 0; i < n; i += len) {
18 let cr = 1;
19 let ci = 0;
20 for (let k = 0; k < len / 2; k++) {
21 const a = i + k;
22 const b = a + len / 2;
23 const tr = re[b] * cr - im[b] * ci;
24 const ti = re[b] * ci + im[b] * cr;
25 re[b] = re[a] - tr;
26 im[b] = im[a] - ti;
27 re[a] += tr;
28 im[a] += ti;
29 const nr = cr * wr - ci * wi;
30 ci = cr * wi + ci * wr;
31 cr = nr;
32 }
33 }
34 }
35}
37export interface Stft {
38 mags: Float32Array[]; // dB per frame, nfft/2 + 1 bins
39 hop: number;
40 nfft: number;
41 fs: number;
42}
44export function stft(x: Float32Array, fs: number, nfft = 2048, hop = 256): Stft {
45 const win = new Float64Array(nfft);
46 for (let i = 0; i < nfft; i++) win[i] = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / nfft);
47 const mags: Float32Array[] = [];
48 const re = new Float64Array(nfft);
49 const im = new Float64Array(nfft);
50 for (let start = -nfft / 2; start < x.length - nfft / 2; start += hop) {
51 for (let i = 0; i < nfft; i++) {
52 const j = start + i;
53 re[i] = j >= 0 && j < x.length ? x[j] * win[i] : 0;
54 im[i] = 0;
55 }
56 fft(re, im);
57 const m = new Float32Array(nfft / 2 + 1);
58 for (let k = 0; k <= nfft / 2; k++) m[k] = 10 * Math.log10(re[k] * re[k] + im[k] * im[k] + 1e-30);
59 mags.push(m);
60 }
61 return { mags, hop, nfft, fs };
62}
64// Power spectrum of the whole signal (Hann window, zero-padded), in dB.
65export function spectrum(x: Float32Array, fs: number): { freq: Float64Array; db: Float64Array } {
66 let n = 1;
67 while (n < x.length) n <<= 1;
68 const re = new Float64Array(n);
69 const im = new Float64Array(n);
70 for (let i = 0; i < x.length; i++) re[i] = x[i] * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / x.length));
71 fft(re, im);
72 const freq = new Float64Array(n / 2 + 1);
73 const db = new Float64Array(n / 2 + 1);
74 for (let k = 0; k <= n / 2; k++) {
75 freq[k] = (k * fs) / n;
76 db[k] = 10 * Math.log10(re[k] * re[k] + im[k] * im[k] + 1e-30);
77 }
78 return { freq, db };
79}