4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1// Minimal iterative radix-2 complex FFT, in-place on separate re/im arrays.
3export function fft(re: Float64Array, im: Float64Array): void {
4 const n = re.length;
5 if ((n & (n - 1)) !== 0) throw new Error('FFT length must be a power of 2');
6 // bit-reversal permutation
7 for (let i = 1, j = 0; i < n; i++) {
8 let bit = n >> 1;
9 for (; j & bit; bit >>= 1) j ^= bit;
10 j ^= bit;
11 if (i < j) {
12 const tr = re[i]; re[i] = re[j]; re[j] = tr;
13 const ti = im[i]; im[i] = im[j]; im[j] = ti;
14 }
15 }
16 for (let len = 2; len <= n; len <<= 1) {
17 const ang = (-2 * Math.PI) / len;
18 const wr = Math.cos(ang);
19 const wi = Math.sin(ang);
20 for (let i = 0; i < n; i += len) {
21 let cwr = 1, cwi = 0;
22 for (let k = 0; k < len / 2; k++) {
23 const ur = re[i + k], ui = im[i + k];
24 const vr = re[i + k + len / 2] * cwr - im[i + k + len / 2] * cwi;
25 const vi = re[i + k + len / 2] * cwi + im[i + k + len / 2] * cwr;
26 re[i + k] = ur + vr;
27 im[i + k] = ui + vi;
28 re[i + k + len / 2] = ur - vr;
29 im[i + k + len / 2] = ui - vi;
30 const nwr = cwr * wr - cwi * wi;
31 cwi = cwr * wi + cwi * wr;
32 cwr = nwr;
33 }
34 }
35 }
36}
38// Magnitude spectrum of a real signal with a Hann window, zero-padded to the
39// next power of two (at least minSize). Returns bins 0..N/2.
40export function magnitudeSpectrum(x: Float32Array | Float64Array, minSize = 16384): {
41 mag: Float64Array;
42 size: number;
43} {
44 let size = minSize;
45 while (size < x.length) size <<= 1;
46 const re = new Float64Array(size);
47 const im = new Float64Array(size);
48 const n = x.length;
49 for (let i = 0; i < n; i++) {
50 const w = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (n - 1)));
51 re[i] = x[i] * w;
52 }
53 fft(re, im);
54 const half = size / 2;
55 const mag = new Float64Array(half + 1);
56 for (let i = 0; i <= half; i++) mag[i] = Math.hypot(re[i], im[i]);
57 return { mag, size };
58}
60// Simple STFT magnitude in dB for the spectrogram, frames x bins.
61export function stft(
62 x: Float32Array,
63 win = 512,
64 hop = 128,
65): { frames: number; bins: number; db: Float32Array; floorDb: number } {
66 const nFrames = Math.max(1, Math.floor((x.length - win) / hop) + 1);
67 const bins = win / 2 + 1;
68 const db = new Float32Array(nFrames * bins);
69 const re = new Float64Array(win);
70 const im = new Float64Array(win);
71 const hann = new Float64Array(win);
72 for (let i = 0; i < win; i++) hann[i] = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (win - 1)));
73 let maxDb = -Infinity;
74 for (let f = 0; f < nFrames; f++) {
75 const off = f * hop;
76 for (let i = 0; i < win; i++) {
77 re[i] = (off + i < x.length ? x[off + i] : 0) * hann[i];
78 im[i] = 0;
79 }
80 fft(re, im);
81 for (let b = 0; b < bins; b++) {
82 const m = Math.hypot(re[b], im[b]);
83 const v = 20 * Math.log10(m + 1e-12);
84 db[f * bins + b] = v;
85 if (v > maxDb) maxDb = v;
86 }
87 }
88 // normalize so 0 dB is the frame-bin maximum
89 for (let i = 0; i < db.length; i++) db[i] -= maxDb;
90 return { frames: nFrames, bins, db, floorDb: -90 };
91}