1// Signal analysis for display: FFT, spectrogram, spectra, mode-energy bands.
3export function fft(re: Float64Array, im: Float64Array) {
4 const n = re.length;
5 for (let i = 1, j = 0; i < n; i++) {
6 let bit = n >> 1;
7 for (; j & bit; bit >>= 1) j ^= bit;
8 j ^= bit;
9 if (i < j) {
10 [re[i], re[j]] = [re[j], re[i]];
11 [im[i], im[j]] = [im[j], im[i]];
12 }
13 }
14 for (let len = 2; len <= n; len <<= 1) {
15 const ang = (-2 * Math.PI) / len;
16 const wr = Math.cos(ang);
17 const wi = Math.sin(ang);
18 for (let i = 0; i < n; i += len) {
19 let cr = 1;
20 let ci = 0;
21 for (let k = 0; k < len / 2; k++) {
22 const a = i + k;
23 const b = a + len / 2;
24 const tr = re[b] * cr - im[b] * ci;
25 const ti = re[b] * ci + im[b] * cr;
26 re[b] = re[a] - tr;
27 im[b] = im[a] - ti;
28 re[a] += tr;
29 im[a] += ti;
30 const nr = cr * wr - ci * wi;
31 ci = cr * wi + ci * wr;
32 cr = nr;
33 }
34 }
35 }
36}
38export interface Spectrogram {
39 nFrames: number;
40 nBins: number;
41 hop: number;
42 win: number;
43 db: Float32Array; // frame-major, nFrames x nBins
44 maxDb: number;
45}
47export function spectrogram(x: Float32Array, win = 2048, hop = 256): Spectrogram {
48 const nFrames = Math.max(1, Math.floor((x.length - win) / hop) + 1);
49 const nBins = win / 2;
50 const db = new Float32Array(nFrames * nBins);
51 const re = new Float64Array(win);
52 const im = new Float64Array(win);
53 const w = new Float64Array(win);
54 for (let i = 0; i < win; i++) w[i] = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / win);
55 let maxDb = -Infinity;
56 for (let f = 0; f < nFrames; f++) {
57 for (let i = 0; i < win; i++) {
58 re[i] = (x[f * hop + i] ?? 0) * w[i];
59 im[i] = 0;
60 }
61 fft(re, im);
62 for (let b = 0; b < nBins; b++) {
63 const v = 10 * Math.log10(re[b] * re[b] + im[b] * im[b] + 1e-300);
64 db[f * nBins + b] = v;
65 if (v > maxDb) maxDb = v;
66 }
67 }
68 return { nFrames, nBins, hop, win, db, maxDb };
69}
71// Power spectrum (dB) of x[start, start+len) with a Hann window, zero padded to a power of two.
72export function spectrum(x: Float32Array, start: number, len: number): { db: Float64Array; df: number; nfft: number } {
73 len = Math.max(2, Math.min(len, x.length - start));
74 let nfft = 1;
75 while (nfft < 2 * len) nfft <<= 1;
76 const re = new Float64Array(nfft);
77 const im = new Float64Array(nfft);
78 for (let i = 0; i < len; i++) re[i] = x[start + i] * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (len - 1)));
79 fft(re, im);
80 const db = new Float64Array(nfft / 2);
81 for (let b = 0; b < nfft / 2; b++) db[b] = 10 * Math.log10(re[b] * re[b] + im[b] * im[b] + 1e-300);
82 return { db, df: 1 / nfft, nfft };
83}
85// Sum mode energies into log-spaced frequency bands.
86export function modeBands(
87 freqs: Float64Array,
88 modeEnergy: Float32Array,
89 nFrames: number,
90 nBands = 48
91): { edges: Float64Array; e: Float32Array; nBands: number } {
92 const nModes = freqs.length;
93 // modes at (nearly) 0 Hz, from free rigid motion, go in the lowest band
94 const fLow = freqs.find(f => f > 1e-6 * freqs[nModes - 1]) ?? freqs[nModes - 1];
95 const fmin = fLow * 0.97;
96 const fmax = freqs[nModes - 1] * 1.03;
97 const edges = new Float64Array(nBands + 1);
98 for (let b = 0; b <= nBands; b++) edges[b] = fmin * Math.pow(fmax / fmin, b / nBands);
99 const band = new Int32Array(nModes);
100 for (let k = 0; k < nModes; k++) {
101 const fk = Math.max(freqs[k], fmin);
102 band[k] = Math.min(nBands - 1, Math.max(0, Math.floor((Math.log(fk / fmin) / Math.log(fmax / fmin)) * nBands)));
103 }
104 const e = new Float32Array(nFrames * nBands);
105 for (let f = 0; f < nFrames; f++) {
106 for (let k = 0; k < nModes; k++) e[f * nBands + band[k]] += modeEnergy[f * nModes + k];
107 }
108 return { edges, e, nBands };
109}
111// Relative RMS difference between two signals over [0, end).
112export function relDiff(a: Float32Array, b: Float32Array, end: number): number {
113 let num = 0;
114 let den = 0;
115 for (let i = 0; i < Math.min(end, a.length, b.length); i++) {
116 num += (a[i] - b[i]) ** 2;
117 den += a[i] ** 2;
118 }
119 return den > 0 ? Math.sqrt(num / den) : 0;
120}
122export function wavBlob(x: Float32Array, fs: number): Blob {
123 let peak = 0;
124 for (const v of x) peak = Math.max(peak, Math.abs(v));
125 const g = peak > 0 ? 0.9 / peak : 1;
126 const buf = new ArrayBuffer(44 + 2 * x.length);
127 const dv = new DataView(buf);
128 const w = (o: number, s: string) => [...s].forEach((c, i) => dv.setUint8(o + i, c.charCodeAt(0)));
129 w(0, "RIFF");
130 dv.setUint32(4, 36 + 2 * x.length, true);
131 w(8, "WAVE");
132 w(12, "fmt ");
133 dv.setUint32(16, 16, true);
134 dv.setUint16(20, 1, true);
135 dv.setUint16(22, 1, true);
136 dv.setUint32(24, fs, true);
137 dv.setUint32(28, 2 * fs, true);
138 dv.setUint16(32, 2, true);
139 dv.setUint16(34, 16, true);
140 w(36, "data");
141 dv.setUint32(40, 2 * x.length, true);
142 for (let i = 0; i < x.length; i++) dv.setInt16(44 + 2 * i, Math.round(32767 * Math.max(-1, Math.min(1, g * x[i]))), true);
143 return new Blob([buf], { type: "audio/wav" });
144}
146// Exponential growth rate of a separation curve: least-squares slope of
147// ln(sep) over the stretch from 100*eps up to 1e-2 (relative to the motion's
148// size). Returns null when the separation never reaches 1e-2 or the stretch is
149// too short to fit.
150export function growthRate(times: Float64Array, sep: Float64Array, eps: number): number | null {
151 let i0 = -1;
152 let i1 = -1;
153 for (let i = 0; i < sep.length; i++) {
154 if (i0 < 0 && sep[i] > 100 * eps) i0 = i;
155 if (i0 >= 0 && sep[i] > 1e-2) {
156 i1 = i;
157 break;
158 }
159 }
160 if (i0 < 0 || i1 < 0 || i1 - i0 < 4) return null;
161 let st = 0, sy = 0, stt = 0, sty = 0;
162 const n = i1 - i0 + 1;
163 for (let i = i0; i <= i1; i++) {
164 const y = Math.log(sep[i]);
165 st += times[i];
166 sy += y;
167 stt += times[i] * times[i];
168 sty += times[i] * y;
169 }
170 return (n * sty - st * sy) / (n * stt - st * st);
171}