2 * The spectrum of what the microphone heard, against where the partials
3 * *should* be — twice over.
4 *
5 * For every object two predictions are drawn. The dashed markers are the
6 * continuous equations' partials (a string's n·f0·sqrt(1+Bn²), a bar's
7 * n²·f1, a membrane's inharmonic mode set). The solid markers are the
8 * *discrete scheme's own* dispersion relation: the same operators the GPU
9 * runs have sine eigenmodes with eigenvalue lam_n = (4/h²) sin²(n π h / 2L)
10 * per axis, so mode n of a string rings at sqrt(cs² lam_n + kap² lam_n²)/2π
11 * exactly. Where the two disagree, that gap *is* the numerical dispersion —
12 * on a bar's coarse grid it is dozens of cents, and you can hear it. The
13 * measured peaks land on the solid markers, which is the point: the solver
14 * is faithful to its scheme, and the scheme is honest about its distance
15 * from the continuum.
16 */
17import type { ObjectGrid } from '../grids.ts';
18import { templateOf } from './objects.ts';
19import type { ObjectSpec } from './types.ts';
21/** Log-magnitude spectrum of a recorded trace, by radix-2 FFT with a Hann
22 * window. Returns dB relative to the peak bin, and the bin spacing. */
23export function computeSpectrum(
24 trace: Float32Array,
25 dt: number,
26 fmax: number,
27): { db: Float32Array; df: number } {
28 // Enough samples for ~5 Hz resolution, capped to keep the FFT instant.
29 let n = 1;
30 while (n < trace.length && n < 1 << 17) n <<= 1;
31 if (n > trace.length) n >>= 1;
32 n = Math.max(n, 2);
33 const re = new Float64Array(n);
34 const im = new Float64Array(n);
35 for (let i = 0; i < n; i++) {
36 const w = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (n - 1));
37 re[i] = trace[i] * w;
38 }
39 fft(re, im);
40 const df = 1 / (n * dt);
41 const bins = Math.max(1, Math.min(n / 2, Math.ceil(fmax / df)));
42 const db = new Float32Array(bins);
43 let peak = 1e-300;
44 for (let i = 0; i < bins; i++) {
45 const m = Math.hypot(re[i], im[i]);
46 db[i] = m;
47 if (m > peak) peak = m;
48 }
49 for (let i = 0; i < bins; i++) db[i] = 20 * Math.log10(Math.max(db[i] / peak, 1e-8));
50 return { db, df };
51}
53/** In-place radix-2 FFT. */
54function fft(re: Float64Array, im: Float64Array): void {
55 const n = re.length;
56 for (let i = 1, j = 0; i < n; i++) {
57 let bit = n >> 1;
58 for (; j & bit; bit >>= 1) j ^= bit;
59 j ^= bit;
60 if (i < j) {
61 [re[i], re[j]] = [re[j], re[i]];
62 [im[i], im[j]] = [im[j], im[i]];
63 }
64 }
65 for (let len = 2; len <= n; len <<= 1) {
66 const ang = (-2 * Math.PI) / len;
67 const wr = Math.cos(ang);
68 const wi = Math.sin(ang);
69 for (let i = 0; i < n; i += len) {
70 let cr = 1;
71 let ci = 0;
72 for (let k = 0; k < len / 2; k++) {
73 const ur = re[i + k];
74 const ui = im[i + k];
75 const vr = re[i + k + len / 2] * cr - im[i + k + len / 2] * ci;
76 const vi = re[i + k + len / 2] * ci + im[i + k + len / 2] * cr;
77 re[i + k] = ur + vr;
78 im[i + k] = ui + vi;
79 re[i + k + len / 2] = ur - vr;
80 im[i + k + len / 2] = ui - vi;
81 const nr = cr * wr - ci * wi;
82 ci = cr * wi + ci * wr;
83 cr = nr;
84 }
85 }
86 }
87}
89export interface PartialSet {
90 id: string;
91 label: string;
92 /** The continuous equations' partials. */
93 continuous: number[];
94 /** The discrete scheme's — where the peaks actually land. */
95 discrete: number[];
96}
98/** The Dirichlet second-difference eigenvalue for mode n. */
99const lam = (n: number, h: number, L: number): number =>
100 (4 / (h * h)) * Math.sin((n * Math.PI * h) / (2 * L)) ** 2;
102/** Both predictions for one object, up to fmax. */
103export function predictedPartials(
104 spec: ObjectSpec,
105 grid: ObjectGrid,
106 params: Record<string, number>,
107 fmax: number,
108): PartialSet {
109 const t = templateOf(spec);
110 const p = (key: string): number => params[`${key}_${spec.id}`] ?? 0;
111 const within = (fs: number[]): number[] => fs.filter((f) => f > 0 && f <= fmax);
112 let discrete: number[] = [];
113 const specNow: ObjectSpec = { ...spec, params: Object.fromEntries(
114 t.paramSpecs.map((s) => [s.key, p(s.key)]),
115 ) };
116 const continuous = within(t.partials(specNow, 16));
118 if (spec.type === 'string') {
119 const cs = 2 * grid.L1 * p('f0');
120 const kap2 = (p('B') * cs * cs * grid.L1 * grid.L1) / (Math.PI * Math.PI);
121 // Modes exist up to n = n1 - 2 (the interior dimension).
122 for (let n = 1; n <= grid.n1 - 2; n++) {
123 const l = lam(n, grid.h1, grid.L1);
124 discrete.push(Math.sqrt(cs * cs * l + kap2 * l * l) / (2 * Math.PI));
125 }
126 } else if (spec.type === 'bar') {
127 const kap = (2 * p('f1') * grid.L1 * grid.L1) / Math.PI;
128 for (let n = 1; n <= grid.n1 - 2; n++) {
129 discrete.push((kap * lam(n, grid.h1, grid.L1)) / (2 * Math.PI));
130 }
131 } else {
132 const L2 = grid.L2;
133 const cm =
134 (2 * p('f0')) / Math.sqrt(1 / (grid.L1 * grid.L1) + 1 / (L2 * L2));
135 for (let m = 1; m <= Math.min(6, grid.n1 - 2); m++) {
136 for (let n = 1; n <= Math.min(6, grid.n2 - 2); n++) {
137 const l = lam(m, grid.h1, grid.L1) + lam(n, grid.h2, grid.L2);
138 discrete.push((cm * Math.sqrt(l)) / (2 * Math.PI));
139 }
140 }
141 discrete.sort((a, b) => a - b);
142 }
143 discrete = within(discrete).slice(0, 16);
144 return { id: spec.id, label: `${t.label} ${spec.id}`, continuous, discrete };
145}
147/** One hue per object, shared by the overlay and the spectrum. */
148export const OBJECT_COLORS = [
149 '#ffd076', '#7fd1ff', '#98e698', '#ff9e9e', '#d2a8ff', '#ffc4e6', '#c9d67f',
150];
151export const objectColor = (index: number): string =>
152 OBJECT_COLORS[index % OBJECT_COLORS.length];
154export interface SpectrumDrawOptions {
155 trace: Float32Array;
156 dt: number;
157 fmax: number;
158 partials: PartialSet[];
159 dark: boolean;
160}
162/** Draw the spectrum with the partial markers into a 2D canvas. */
163export function drawSpectrum(canvas: HTMLCanvasElement, o: SpectrumDrawOptions): void {
164 const ctx = canvas.getContext('2d');
165 if (!ctx) return;
166 const dpr = devicePixelRatio;
167 const rect = canvas.getBoundingClientRect();
168 const w = Math.max(1, Math.round(rect.width * dpr));
169 const h = Math.max(1, Math.round(rect.height * dpr));
170 if (canvas.width !== w || canvas.height !== h) {
171 canvas.width = w;
172 canvas.height = h;
173 }
174 ctx.clearRect(0, 0, w, h);
175 const ink = o.dark ? 'rgba(230,233,236,0.92)' : 'rgba(31,35,40,0.92)';
176 const faint = o.dark ? 'rgba(154,164,175,0.3)' : 'rgba(87,96,106,0.3)';
177 const padL = 34 * dpr;
178 const padB = 18 * dpr;
179 const padT = 6 * dpr;
180 const plotW = w - padL - 6 * dpr;
181 const plotH = h - padT - padB;
182 const dbFloor = -90;
183 const fx = (f: number): number => padL + (plotW * f) / o.fmax;
184 const fy = (db: number): number =>
185 padT + plotH * Math.min(1, Math.max(0, -db / -dbFloor));
187 // Grid lines and labels every kilohertz (every 200 Hz under 1.5 kHz).
188 ctx.font = `${10 * dpr}px system-ui, sans-serif`;
189 ctx.fillStyle = faint;
190 ctx.strokeStyle = faint;
191 ctx.lineWidth = dpr;
192 const tick = o.fmax > 1500 ? 1000 : 200;
193 for (let f = tick; f < o.fmax; f += tick) {
194 ctx.beginPath();
195 ctx.moveTo(fx(f), padT);
196 ctx.lineTo(fx(f), padT + plotH);
197 ctx.stroke();
198 ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, fx(f) + 2 * dpr, h - 6 * dpr);
199 }
200 for (const db of [0, -30, -60]) {
201 ctx.fillText(`${db}`, 4 * dpr, fy(db) + 3 * dpr);
202 }
204 // Partial markers behind the curve: solid = the scheme's own modes,
205 // dashed = the continuous equations'.
206 o.partials.forEach((set, i) => {
207 const color = objectColor(i);
208 ctx.strokeStyle = color;
209 ctx.globalAlpha = 0.85;
210 ctx.setLineDash([]);
211 ctx.lineWidth = 1.4 * dpr;
212 for (const f of set.discrete) {
213 if (f > o.fmax) continue;
214 ctx.beginPath();
215 ctx.moveTo(fx(f), padT);
216 ctx.lineTo(fx(f), padT + plotH * 0.98);
217 ctx.stroke();
218 }
219 ctx.setLineDash([3 * dpr, 4 * dpr]);
220 ctx.globalAlpha = 0.55;
221 ctx.lineWidth = dpr;
222 for (const f of set.continuous) {
223 if (f > o.fmax) continue;
224 ctx.beginPath();
225 ctx.moveTo(fx(f), padT);
226 ctx.lineTo(fx(f), padT + plotH * 0.98);
227 ctx.stroke();
228 }
229 ctx.setLineDash([]);
230 ctx.globalAlpha = 1;
231 });
233 // The measured spectrum.
234 const { db, df } = computeSpectrum(o.trace, o.dt, o.fmax);
235 ctx.strokeStyle = ink;
236 ctx.lineWidth = 1.4 * dpr;
237 ctx.beginPath();
238 for (let i = 0; i < db.length; i++) {
239 const xx = fx(i * df);
240 const yy = fy(db[i]);
241 if (i === 0) ctx.moveTo(xx, yy);
242 else ctx.lineTo(xx, yy);
243 }
244 ctx.stroke();
245}