1/**
2 * A plain radix-2 FFT, enough for the spectrum view: Hann-windowed
3 * magnitude of the first power-of-two chunk of the trace, in dB below the
4 * peak. Nothing here is performance-critical; a quarter-million-point
5 * transform is well under a frame.
6 */
8/** In-place complex FFT, re/im interleaved as two arrays, n a power of 2. */
9export function fft(re: Float64Array, im: Float64Array): void {
10 const n = re.length;
11 for (let i = 1, j = 0; i < n; i++) {
12 let bit = n >> 1;
13 for (; j & bit; bit >>= 1) j ^= bit;
14 j ^= bit;
15 if (i < j) {
16 [re[i], re[j]] = [re[j], re[i]];
17 [im[i], im[j]] = [im[j], im[i]];
18 }
19 }
20 for (let len = 2; len <= n; len <<= 1) {
21 const ang = (-2 * Math.PI) / len;
22 const wr = Math.cos(ang);
23 const wi = Math.sin(ang);
24 for (let i = 0; i < n; i += len) {
25 let cr = 1;
26 let ci = 0;
27 for (let k = 0; k < len / 2; k++) {
28 const ur = re[i + k];
29 const ui = im[i + k];
30 const vr = re[i + k + len / 2] * cr - im[i + k + len / 2] * ci;
31 const vi = re[i + k + len / 2] * ci + im[i + k + len / 2] * cr;
32 re[i + k] = ur + vr;
33 im[i + k] = ui + vi;
34 re[i + k + len / 2] = ur - vr;
35 im[i + k + len / 2] = ui - vi;
36 const nr = cr * wr - ci * wi;
37 ci = cr * wi + ci * wr;
38 cr = nr;
39 }
40 }
41 }
42}
44export interface Spectrum {
45 /** Magnitude in dB relative to the strongest bin (all <= 0). */
46 db: Float32Array;
47 /** Hz per bin. */
48 binHz: number;
49}
51/** The spectrum of a trace sampled at 1/dt: Hann window over the first
52 * power-of-two chunk (at most 2^18 samples). */
53export function spectrumOf(trace: Float32Array, dt: number): Spectrum {
54 let n = 1;
55 while (n * 2 <= Math.min(trace.length, 1 << 18)) n *= 2;
56 const re = new Float64Array(n);
57 const im = new Float64Array(n);
58 for (let i = 0; i < n; i++) {
59 re[i] = trace[i] * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (n - 1)));
60 }
61 fft(re, im);
62 const half = n / 2;
63 const mag = new Float64Array(half);
64 let peak = 0;
65 for (let i = 0; i < half; i++) {
66 mag[i] = Math.hypot(re[i], im[i]);
67 if (mag[i] > peak) peak = mag[i];
68 }
69 const db = new Float32Array(half);
70 for (let i = 0; i < half; i++) {
71 db[i] = peak > 0 ? 20 * Math.log10(Math.max(mag[i] / peak, 1e-8)) : -160;
72 }
73 return { db, binHz: 1 / (n * dt) };
74}