1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 1// Incremental short-time Fourier transform of the listener signal: a Hann window of nfft samples
2// advanced by hop, one power column per hop, plus the spectral centroid of each column.
4/** In-place iterative radix-2 FFT of (re, im); n must be a power of two. */
5export function fft(re: Float32Array, im: Float32Array) {
6 const n = re.length;
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) { [re[i], re[j]] = [re[j], re[i]]; [im[i], im[j]] = [im[j], im[i]]; }
12 }
13 for (let len = 2; len <= n; len <<= 1) {
14 const ang = (-2 * Math.PI) / len, wr = Math.cos(ang), wi = Math.sin(ang);
15 for (let i = 0; i < n; i += len) {
16 let cr = 1, ci = 0;
17 for (let k = 0; k < len / 2; k++) {
18 const a = i + k, b = a + len / 2;
19 const tr = re[b] * cr - im[b] * ci, ti = re[b] * ci + im[b] * cr;
20 re[b] = re[a] - tr; im[b] = im[a] - ti; re[a] += tr; im[a] += ti;
21 const nr = cr * wr - ci * wi; ci = cr * wi + ci * wr; cr = nr;
22 }
23 }
24 }
25}
27export class Stft {
28 readonly nfft: number; readonly hop: number; readonly fs: number; readonly nBins: number;
29 readonly columns: Float32Array[] = []; // power per bin
30 readonly centroid: number[] = []; // Hz per column
31 private readonly window: Float32Array;
32 private buffer = new Float32Array(0);
33 private consumed = 0; // samples already removed from the front of buffer
34 private readonly re: Float32Array; private readonly im: Float32Array;
36 constructor(fs: number, nfft = 1024, hop = 256) {
37 this.fs = fs; this.nfft = nfft; this.hop = hop; this.nBins = nfft / 2 + 1;
38 this.window = new Float32Array(nfft);
39 for (let i = 0; i < nfft; i++) this.window[i] = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / nfft);
40 this.re = new Float32Array(nfft); this.im = new Float32Array(nfft);
41 }
43 /** Time of column c (its window centre) in seconds. */
44 columnTime(c: number): number { return (c * this.hop + this.nfft / 2) / this.fs; }
45 binFrequency(b: number): number { return (b * this.fs) / this.nfft; }
47 /** Append samples; returns the number of new columns. */
48 push(samples: Float32Array): number {
49 const joined = new Float32Array(this.buffer.length + samples.length);
50 joined.set(this.buffer); joined.set(samples, this.buffer.length);
51 this.buffer = joined;
52 let n = 0, start = 0;
53 while (start + this.nfft <= this.buffer.length) {
54 for (let i = 0; i < this.nfft; i++) { this.re[i] = this.buffer[start + i] * this.window[i]; this.im[i] = 0; }
55 fft(this.re, this.im);
56 const col = new Float32Array(this.nBins);
57 let num = 0, den = 0;
58 for (let b = 0; b < this.nBins; b++) {
59 col[b] = this.re[b] * this.re[b] + this.im[b] * this.im[b];
60 num += col[b] * this.binFrequency(b); den += col[b];
61 }
62 this.columns.push(col); this.centroid.push(den > 0 ? num / den : 0);
63 start += this.hop; n++;
64 }
65 this.buffer = this.buffer.slice(start); this.consumed += start;
66 return n;
67 }
68}