1/**
2 * Every object in closed form: exact modes of the discrete schemes, and the
3 * mic trace synthesized from them against precomputed room kernels.
4 *
5 * This is the mathematics behind the fast render path (src/studio/room.ts
6 * records the kernels it consumes; src/studio/fast.ts strings the two
7 * together). The facts it rests on, verified numerically in
8 * scripts/modal-poc.ts and test/checks.ts:
9 *
10 * - Every object's discrete operator has exact sine eigenmodes. The pinned
11 * 5-point Laplacian's are product sines; the second difference's are
12 * sines with lam_m = (4/h²) sin²(mπ/2(n−1)); and the fourth difference
13 * with antisymmetric ghosts is exactly the square of the pinned second
14 * difference (the README's own point about it), so it shares those sines
15 * with eigenvalue lam². The damping terms are built from the same
16 * operators. Each modal amplitude therefore obeys a scalar two-term
17 * recursion — the leapfrog restricted to one mode — whose solution is a
18 * closed-form geometric series. No object is ever timestepped.
19 * - The coupling into the air is linear and one-way, so every source the
20 * air sees is a sum of *fixed spatial patterns* times scalar series in
21 * those modal amplitudes: an object's rad route is one pattern per mode
22 * (prof .* rad(phi_m)), and a string's bridge route is a single pattern
23 * (bprof) times the slope its modes arrive at the far end with.
24 * - The whole system is LTI, so the mic trace is the sum over (pattern,
25 * series) pairs of the series convolved with that pattern's impulse
26 * response of the room — the kernel, which depends only on geometry.
27 *
28 * Each series being exactly Re(c·z^k) per characteristic root, its
29 * convolution reduces to a complex one-pole recursion over the kernel: the
30 * whole synthesis is O(poles × samples) scalar work on the CPU, with the
31 * closed-form exponential tail falling out for free once the kernel ends.
32 *
33 * Everything here matches the compiled f32 schemes to their own rounding
34 * noise; the approximations are the membrane mode set kept (the caller's
35 * eps/fmax — line objects have few modes and always keep all) and the
36 * kernel length.
37 */
38import type { AirGrid, ObjectGrid } from '../grids.ts';
39import type { ObjectSpec, Params } from './types.ts';
40import { suffixed } from './types.ts';
42/** One root of a mode's recursion, carrying the displacement series:
43 * alpha[k] += Re(a · z^k). A real root is the zIm = 0 case. */
44export interface BasePole {
45 aRe: number;
46 aIm: number;
47 zRe: number;
48 zIm: number;
49}
51/** One term of a source series, ready to filter a kernel: Re(c · z^k). */
52export interface SeriesPole {
53 cRe: number;
54 cIm: number;
55 zRe: number;
56 zIm: number;
57}
59/** One mode of one object: which sine it is (n = 0 on a line grid), its
60 * discrete frequency and decay under the current parameters, and its
61 * displacement series as poles. */
62export interface ModalMode {
63 /** Sine indices, 1-based; a line mode has n = 0. */
64 m: number;
65 n: number;
66 /** Eigenvalue of minus the discrete second difference / Laplacian. */
67 lam: number;
68 /** Discrete frequency (Hz) and 1/e decay time (s) under these params. */
69 freq: number;
70 tau: number;
71 /** Peak source strength |acc coefficient|, for truncation decisions. */
72 amp: number;
73 poles: BasePole[];
74}
76/** sin(mπi/(n−1)) for m = 1..n−2, i = 0..n−1, row-major by m. */
77export function sineTable(n: number): Float64Array {
78 const t = new Float64Array((n - 2) * n);
79 for (let m = 1; m <= n - 2; m++)
80 for (let i = 0; i < n; i++) t[(m - 1) * n + i] = Math.sin((m * Math.PI * i) / (n - 1));
81 return t;
82}
84/* ----------------------------------------------------- the recursions --- */
86/**
87 * Solve one mode's recursion in closed form.
88 *
89 * Restricted to an eigenvalue, every object's step is
90 * a[k+1] = cA·a[k] + cB·a[k−1],
91 * and the excitation is the pair (a[0], a[−1]): a strike is a velocity
92 * (a[0] = 0, a[−1] = b), a pluck a held displacement (a[0] = a[−1] = b).
93 * The characteristic roots give a[k] as a geometric series; underdamped
94 * modes (all the audible ones) as one conjugate pair written Re(a·z^k),
95 * modes pushed past critical damping near the grid's Nyquist as two real
96 * roots.
97 */
98export function solveBasePoles(
99 cA: number,
100 cB: number,
101 a0: number,
102 am1: number,
103): BasePole[] {
104 const disc = cA * cA + 4 * cB;
105 if (disc < 0) {
106 const zRe = cA / 2;
107 const zIm = Math.sqrt(-disc) / 2;
108 const zAbs2 = zRe * zRe + zIm * zIm;
109 // alpha[0] = Re(a), alpha[-1] = Re(a·conj(z))/|z|².
110 const aRe = a0;
111 const aIm = zIm !== 0 ? (am1 * zAbs2 - aRe * zRe) / zIm : 0;
112 return [{ aRe, aIm, zRe, zIm }];
113 }
114 const sq = Math.sqrt(Math.max(disc, 1e-300));
115 const r1 = (cA + sq) / 2;
116 const r2 = (cA - sq) / 2;
117 if (r1 === r2 || r1 === 0 || r2 === 0) return [];
118 // a1 + a2 = a0; a1/r1 + a2/r2 = am1.
119 const det = 1 / r1 - 1 / r2;
120 const a1 = (am1 - a0 / r2) / det;
121 const a2 = a0 - a1;
122 return [
123 { aRe: a1, aIm: 0, zRe: r1, zIm: 0 },
124 { aRe: a2, aIm: 0, zRe: r2, zIm: 0 },
125 ];
126}
128const freqOf = (poles: BasePole[], dt: number): { freq: number; tau: number } => {
129 if (poles.length === 0) return { freq: 0, tau: 0 };
130 const p = poles[0];
131 if (p.zIm !== 0) {
132 const zAbs = Math.hypot(p.zRe, p.zIm);
133 return {
134 freq: Math.abs(Math.atan2(p.zIm, p.zRe)) / (2 * Math.PI * dt),
135 tau: zAbs < 1 ? -dt / Math.log(zAbs) : Infinity,
136 };
137 }
138 const rMax = Math.max(...poles.map((q) => Math.abs(q.zRe)));
139 return {
140 freq: poles.some((q) => q.zRe < 0) ? 1 / (2 * dt) : 0,
141 tau: rMax > 0 && rMax < 1 ? -dt / Math.log(rMax) : Infinity,
142 };
143};
145/** Complex-multiply a base pole's coefficient by f(z). */
146const scalePole = (p: BasePole, fRe: number, fIm: number): SeriesPole => ({
147 cRe: p.aRe * fRe - p.aIm * fIm,
148 cIm: p.aRe * fIm + p.aIm * fRe,
149 zRe: p.zRe,
150 zIm: p.zIm,
151});
153/** The acceleration series (un − 2u + um)/dt² of a mode, as poles:
154 * coefficient scaled by (z − 2 + 1/z)/dt² per root. */
155export function accPoles(mode: ModalMode, dt: number, gain = 1): SeriesPole[] {
156 return mode.poles.map((p) => {
157 const zAbs2 = p.zRe * p.zRe + p.zIm * p.zIm;
158 const fRe = (gain * (p.zRe - 2 + p.zRe / zAbs2)) / (dt * dt);
159 const fIm = (gain * (p.zIm - p.zIm / zAbs2)) / (dt * dt);
160 return scalePole(p, fRe, fIm);
161 });
162}
164/** The series alpha[k+1] (the object's *new* displacement at step k, which
165 * is what a string's bridge slope is read from): coefficient scaled by z. */
166export function nextPoles(mode: ModalMode, gain = 1): SeriesPole[] {
167 return mode.poles.map((p) => scalePole(p, gain * p.zRe, gain * p.zIm));
168}
170/* ------------------------------------------------------ the mode sets --- */
172export interface ModeSetOptions {
173 /** Membranes only: keep modes whose source strength is at least eps × the
174 * strongest. 0 keeps everything the strike touches. Line objects have
175 * few modes and always keep all. */
176 eps?: number;
177 /** Membranes only: keep modes below this discrete frequency, Hz. */
178 fmax?: number;
179}
181const ampOf = (poles: BasePole[], dt: number): number => {
182 let amp = 0;
183 for (const p of accPoles({ poles } as ModalMode, dt)) amp += Math.hypot(p.cRe, p.cIm);
184 return amp;
185};
187/** Eigenvalue of minus the pinned second difference for sine mode m. */
188export const lineLam = (g: ObjectGrid, m: number): number => {
189 const s = Math.sin((m * Math.PI) / (2 * (g.n1 - 1)));
190 return (4 / (g.h1 * g.h1)) * s * s;
191};
193/** Eigenvalue of minus the pinned 5-point Laplacian for sine mode (m, n). */
194export const membraneLam = (g: ObjectGrid, m: number, n: number): number => {
195 const s1 = Math.sin((m * Math.PI) / (2 * (g.n1 - 1)));
196 const s2 = Math.sin((n * Math.PI) / (2 * (g.n2 - 1)));
197 return (4 / (g.h1 * g.h1)) * s1 * s1 + (4 / (g.h2 * g.h2)) * s2 * s2;
198};
200/** Project a line field onto every interior sine mode, normalized. */
201export function projectLine(g: ObjectGrid, v: Float64Array): Float64Array {
202 const n = g.n1;
203 const t = sineTable(n);
204 const out = new Float64Array(n - 2);
205 for (let m = 0; m < n - 2; m++) {
206 let s = 0;
207 for (let i = 0; i < n; i++) s += t[m * n + i] * v[i];
208 out[m] = s / ((n - 1) / 2);
209 }
210 return out;
211}
213/** Project a membrane field onto every interior sine mode: two separable
214 * contractions, out[(m−1) + (n1−2)(n−1)], normalized. */
215export function projectModes(g: ObjectGrid, v: Float64Array): Float64Array {
216 const { n1, n2 } = g;
217 const M1 = n1 - 2;
218 const M2 = n2 - 2;
219 const s1 = sineTable(n1);
220 const s2 = sineTable(n2);
221 const tmp = new Float64Array(M1 * n2);
222 for (let m = 0; m < M1; m++) {
223 for (let j = 0; j < n2; j++) {
224 let s = 0;
225 for (let i = 0; i < n1; i++) s += s1[m * n1 + i] * v[i + n1 * j];
226 tmp[m + M1 * j] = s;
227 }
228 }
229 const out = new Float64Array(M1 * M2);
230 const norm = ((n1 - 1) / 2) * ((n2 - 1) / 2);
231 for (let m = 0; m < M1; m++) {
232 for (let n = 0; n < M2; n++) {
233 let s = 0;
234 for (let j = 0; j < n2; j++) s += s2[n * n2 + j] * tmp[m + M1 * j];
235 out[m + M1 * n] = s / norm;
236 }
237 }
238 return out;
239}
241/** One object's parameter, unsuffixed, defaulting hard on absence. */
242const par = (params: Params, id: string, key: string): number => {
243 const v = params[suffixed(key, id)];
244 if (v === undefined) throw new Error(`fast path: parameter '${suffixed(key, id)}' missing`);
245 return v;
246};
248/**
249 * The mode set for one object under its current parameters, sorted
250 * strongest first: everything the synthesis needs. For strings and bars
251 * every interior mode is kept (there are a few dozen); for membranes
252 * `opts` bounds the set, which with the kernel length is the fast path's
253 * only approximation.
254 *
255 * The excitations are exactly the generated init's: a string's triangular
256 * pluck held then released (a[0] = a[−1]), a bar's and membrane's
257 * mallet-shaped initial velocity (a[0] = 0).
258 */
259export function objectModeSet(
260 obj: ObjectSpec,
261 g: ObjectGrid,
262 dt: number,
263 params: Params,
264 opts: ModeSetOptions = {},
265): ModalMode[] {
266 const id = obj.id;
267 const t60 = par(params, id, 't60');
268 const sig1 = par(params, id, 'sig1');
269 const d = (6.91 / t60) * dt;
270 /** cA/cB from the stiffness's squared angular frequency at eigenvalue lam. */
271 const coeffs = (w2: number, lam: number): { cA: number; cB: number } => ({
272 cA: (2 - dt * dt * w2 - 2 * sig1 * dt * lam) / (1 + d),
273 cB: (-(1 - d) + 2 * sig1 * dt * lam) / (1 + d),
274 });
276 if (obj.type === 'membrane') {
277 const eps = opts.eps ?? 1e-6;
278 const fmax = opts.fmax ?? Infinity;
279 const f0 = par(params, id, 'f0');
280 const cm = (2 * f0) / Math.sqrt(1 / (g.L1 * g.L1) + 1 / (g.L2 * g.L2));
281 const um0 = new Float64Array(g.npts);
282 const sx0 = par(params, id, 'strikex') * g.L1;
283 const sy0 = par(params, id, 'strikey') * g.L2;
284 const sr = par(params, id, 'striker');
285 const vel = par(params, id, 'vel');
286 for (let k = 0; k < g.npts; k++) {
287 const dx = g.xs64[k] - sx0;
288 const dy = g.ys64[k] - sy0;
289 um0[k] = -(dt * vel) * g.pin[k] * Math.exp(-(dx * dx + dy * dy) / (sr * sr));
290 }
291 const b = projectModes(g, um0);
292 const M1 = g.n1 - 2;
293 const all: ModalMode[] = [];
294 for (let n = 1; n <= g.n2 - 2; n++) {
295 for (let m = 1; m <= M1; m++) {
296 const bmn = b[m - 1 + M1 * (n - 1)];
297 if (bmn === 0) continue;
298 const lam = membraneLam(g, m, n);
299 const { cA, cB } = coeffs(cm * cm * lam, lam);
300 const poles = solveBasePoles(cA, cB, 0, bmn);
301 const { freq, tau } = freqOf(poles, dt);
302 if (poles.length === 0 || freq > fmax) continue;
303 all.push({ m, n, lam, freq, tau, amp: ampOf(poles, dt), poles });
304 }
305 }
306 all.sort((a, bb) => bb.amp - a.amp);
307 const ampMax = all.length ? all[0].amp : 0;
308 return all.filter((md) => md.amp >= eps * ampMax);
309 }
311 // A line object: pluck (string) or strike (bar), then the same recursion
312 // with the type's own stiffness.
313 const n = g.n1;
314 const u0 = new Float64Array(n);
315 let w2: (lam: number) => number;
316 if (obj.type === 'string') {
317 const f0 = par(params, id, 'f0');
318 const B = par(params, id, 'B');
319 const cs = 2 * g.L1 * f0;
320 const kap2 = (B * (cs * g.L1) ** 2) / (Math.PI * Math.PI);
321 w2 = (lam) => cs * cs * lam + kap2 * lam * lam;
322 const pp = par(params, id, 'pluckpos') * g.L1;
323 const amp = par(params, id, 'amp');
324 for (let i = 0; i < n; i++) {
325 const a = g.xs64[i] / pp;
326 const bb = (g.L1 - g.xs64[i]) / (g.L1 - pp);
327 u0[i] = amp * 0.5 * (a + bb - Math.abs(a - bb));
328 }
329 } else {
330 const f1 = par(params, id, 'f1');
331 const kap = (2 * f1 * g.L1 * g.L1) / Math.PI;
332 w2 = (lam) => kap * kap * lam * lam;
333 const sp = par(params, id, 'strikepos') * g.L1;
334 const sr = par(params, id, 'striker');
335 const vel = par(params, id, 'vel');
336 for (let i = 0; i < n; i++) {
337 const r = (g.xs64[i] - sp) / sr;
338 u0[i] = -(dt * vel) * g.pin[i] * Math.exp(-r * r);
339 }
340 }
341 const b = projectLine(g, u0);
342 const out: ModalMode[] = [];
343 for (let m = 1; m <= n - 2; m++) {
344 if (b[m - 1] === 0) continue;
345 const lam = lineLam(g, m);
346 const { cA, cB } = coeffs(w2(lam), lam);
347 // A pluck starts held (a[0] = a[−1]); a strike starts moving (a[0] = 0).
348 const poles =
349 obj.type === 'string'
350 ? solveBasePoles(cA, cB, b[m - 1], b[m - 1])
351 : solveBasePoles(cA, cB, 0, b[m - 1]);
352 if (poles.length === 0) continue;
353 const { freq, tau } = freqOf(poles, dt);
354 out.push({ m, n: 0, lam, freq, tau, amp: ampOf(poles, dt), poles });
355 }
356 return out;
357}
359/* -------------------------------------------------------- air tables ---- */
361/** The bilinearly interpolated sine tables the coupling ops effectively
362 * apply, per air column: rad of a product of sines factorizes per axis,
363 * and a line object's rad is the x factor alone. */
364export function radSineTables(
365 g: ObjectGrid,
366 air: AirGrid,
367): { sx: Float32Array; sy: Float32Array; M1: number; M2: number } {
368 const { n1, n2, h1, h2, L1, L2 } = g;
369 const M1 = n1 - 2;
370 const s1 = sineTable(n1);
371 const sx = new Float32Array(M1 * air.nx);
372 for (let ix = 0; ix < air.nx; ix++) {
373 const xw = -air.Lx / 2 + (ix + 0.5) * air.h;
374 const u = (xw + L1 / 2) / h1;
375 const j = Math.floor(u);
376 if (j < 0 || j >= n1 - 1) continue;
377 const f = u - j;
378 for (let m = 0; m < M1; m++)
379 sx[m * air.nx + ix] = (1 - f) * s1[m * n1 + j] + f * s1[m * n1 + j + 1];
380 }
381 if (g.kind === 'line') {
382 // No second axis: the y factor is 1 everywhere (the tube profile in the
383 // scene's prof carries the y and z structure).
384 return { sx, sy: new Float32Array(air.ny).fill(1), M1, M2: 1 };
385 }
386 const M2 = n2 - 2;
387 const s2 = sineTable(n2);
388 const sy = new Float32Array(M2 * air.ny);
389 for (let iy = 0; iy < air.ny; iy++) {
390 const yw = -air.Ly / 2 + (iy + 0.5) * air.h;
391 const v = (yw + L2 / 2) / h2;
392 const j = Math.floor(v);
393 if (j < 0 || j >= n2 - 1) continue;
394 const f = v - j;
395 for (let n = 0; n < M2; n++)
396 sy[n * air.ny + iy] = (1 - f) * s2[n * n2 + j] + f * s2[n * n2 + j + 1];
397 }
398 return { sx, sy, M1, M2 };
399}
401/* --------------------------------------------------------- synthesis ---- */
403/** One (series, kernel) pair: the poles' series convolved with the kernel. */
404export interface Contribution {
405 poles: SeriesPole[];
406 /** One kernel row, `kernelSteps` samples. */
407 kernel: Float32Array;
408}
410/**
411 * Synthesize the mic trace: every contribution's kernel filtered by its
412 * poles, summed. The trace may be longer than the kernels: a pole keeps
413 * ringing the room's response down after the kernel ends, which is exactly
414 * the closed-form tail.
415 */
416export function synthesize(
417 contribs: Contribution[],
418 kernelSteps: number,
419 totalSteps: number,
420): Float32Array {
421 const y = new Float64Array(totalSteps);
422 for (const c of contribs) {
423 const h = c.kernel;
424 for (const pl of c.poles) {
425 const { cRe, cIm, zRe, zIm } = pl;
426 let wRe = 0;
427 let wIm = 0;
428 let wPk = 0;
429 const lim = Math.min(kernelSteps, totalSteps);
430 for (let t = 0; t < lim; t++) {
431 const nRe = zRe * wRe - zIm * wIm + h[t];
432 const nIm = zRe * wIm + zIm * wRe;
433 wRe = nRe;
434 wIm = nIm;
435 const m2 = wRe * wRe + wIm * wIm;
436 if (m2 > wPk) wPk = m2;
437 y[t] += cRe * wRe - cIm * wIm;
438 }
439 // The closed-form tail: the pole rings the kernel's response down with
440 // no further input. Cut it 240 dB under its own peak — far below both
441 // f32 and the .wav — which also keeps the recursion out of subnormal
442 // floats, where every multiply is slow.
443 const cut = Math.max(wPk * 1e-24, 1e-280);
444 for (let t = lim; t < totalSteps; t++) {
445 const nRe = zRe * wRe - zIm * wIm;
446 const nIm = zRe * wIm + zIm * wRe;
447 wRe = nRe;
448 wIm = nIm;
449 if (wRe * wRe + wIm * wIm < cut) break;
450 y[t] += cRe * wRe - cIm * wIm;
451 }
452 }
453 }
454 const out = new Float32Array(totalSteps);
455 for (let t = 0; t < totalSteps; t++) out[t] = y[t];
456 return out;
457}