1// Modal synthesis of a strike. Each mode m of the linearized plate is a
2// damped oscillator q'' + 2 sigma q' + omega^2 q = f, with sigma the sum of
3// the metal's internal friction, the thermoelastic loss, the radiation of
4// the mode (damping.ts). The stick is
5// the Hertz contact of the full simulation, stepped against the modal
6// system with an exact update of every oscillator under a piecewise
7// constant force. Once the stick has left, every mode is a damped sinusoid
8// with a known complex amplitude, and the delayed Rayleigh sum over the
9// plate collapses to one complex gain per mode. During the attack, while the
10// modes are being driven, the delayed sum is evaluated directly from a
11// table of the nodal accelerations. The parts that are sums over the modes
12// (gains, the acceleration table, the audio tail, the movie) accept a model
13// restricted to a subset of the modes, so that a pool of workers can split
14// the modes between them (synthPool.ts).
16import { type LinearSystem } from "../model/assemble.ts";
17import { type ModeSet } from "./modes.ts";
18import { type Plate, type SimParams, STICK_RELEASE, listenerTerms, strikePatch, nodeMass } from "../params.ts";
20/** What the synthesis needs of a plate and its modes; plain arrays so that it can be posted to a worker. */
21export interface ModalModel {
22 nGrid: number;
23 /** Number of plate nodes (rows of wShapes). */
24 nRows: number;
25 count: number;
26 freqs: Float64Array;
27 /** Transverse component of the mass-normalized mode shapes, node-major: wShapes[row * count + m]. */
28 wShapes: Float32Array;
29 /** Row of each grid node in wShapes, or -1 outside the plate. */
30 wRow: Int32Array;
31 /** Grid node of each row. */
32 rowNode: Int32Array;
33 /** Decay rate of each mode from the metal's internal friction [1/s]. */
34 sigmaMat: Float64Array;
35 /** Decay rate of each mode from thermoelastic damping [1/s]; filled by damping.ts, zero until then. */
36 sigmaTe: Float64Array;
37 /** Decay rate of each mode from its radiation into the air [1/s]; filled by damping.ts, zero until then. */
38 sigmaRad: Float64Array;
39 /** Node masses on the grid [kg]. */
40 masses: Float64Array;
41}
43export function modalModel(pl: Plate, sys: LinearSystem, modes: ModeSet): ModalModel {
44 const nGrid = pl.grid.n;
45 const wRow = new Int32Array(nGrid).fill(-1);
46 const rows: number[] = [];
47 for (let k = 0; k < nGrid; k++) if (sys.nodeDof[k] >= 0) { wRow[k] = rows.length; rows.push(k); }
48 const nRows = rows.length, M = modes.count;
49 const wShapes = new Float32Array(nRows * M);
50 for (let r = 0; r < nRows; r++) {
51 const d = sys.nodeDof[rows[r]] + (sys.fields - 1);
52 for (let m = 0; m < M; m++) wShapes[r * M + m] = modes.shapes[m * sys.n + d];
53 }
54 const masses = new Float64Array(nGrid);
55 for (let j = 0; j < pl.grid.ny; j++) for (let i = 0; i < pl.grid.nx; i++) masses[j * pl.grid.nx + i] = nodeMass(pl, i, j);
56 // Hysteretic damping: the metal's loss factor eta on the mode's strain
57 // energy, omega^2 / 2 for a mass-normalized mode, gives the rate
58 // sigma = eta omega / 2. The thermoelastic and the radiation rates need
59 // the shapes on the grid and are filled in by damping.ts.
60 const sigmaMat = new Float64Array(M);
61 for (let m = 0; m < M; m++) {
62 const omega = 2 * Math.PI * modes.freqs[m];
63 if (!(omega > 0)) continue;
64 sigmaMat[m] = (pl.mat.eta * omega) / 2;
65 }
66 return { nGrid, nRows, count: M, freqs: modes.freqs, wShapes, wRow, rowNode: Int32Array.from(rows), sigmaMat, sigmaTe: new Float64Array(M), sigmaRad: new Float64Array(M), masses };
67}
69/** The model restricted to the modes [from, to). */
70export function sliceModel(model: ModalModel, from: number, to: number): ModalModel {
71 const M = to - from;
72 const wShapes = new Float32Array(model.nRows * M);
73 for (let r = 0; r < model.nRows; r++) wShapes.set(model.wShapes.subarray(r * model.count + from, r * model.count + to), r * M);
74 return { ...model, count: M, freqs: model.freqs.slice(from, to), wShapes, sigmaMat: model.sigmaMat.slice(from, to), sigmaTe: model.sigmaTe.slice(from, to), sigmaRad: model.sigmaRad.slice(from, to) };
75}
77/** The transverse shape of mode m on the grid (zero outside the plate). */
78export function modeShapeOnGrid(model: ModalModel, m: number): Float32Array {
79 const w = new Float32Array(model.nGrid);
80 for (let r = 0; r < model.nRows; r++) w[model.rowNode[r]] = model.wShapes[r * model.count + m];
81 return w;
82}
84/** Per-mode decay rates: internal friction, thermoelastic and radiation, kept below the mode's frequency (underdamped). */
85export function modalRates(model: ModalModel): Float64Array {
86 const sigma = new Float64Array(model.count);
87 for (let m = 0; m < model.count; m++) {
88 const w = 2 * Math.PI * model.freqs[m];
89 sigma[m] = Math.min(model.sigmaMat[m] + model.sigmaTe[m] + model.sigmaRad[m], 0.99 * w);
90 }
91 return sigma;
92}
94/** Damped angular frequency of a mode. */
95function dampedOmega(f: number, sigma: number): number {
96 const w = 2 * Math.PI * f;
97 return Math.sqrt(w * w - sigma * sigma);
98}
100/** The contact: modal gains g_m = sum over the patch of (m_k / M_patch) phi_m(k), so that the modal force is -F g_m and the contact coordinate is sum_m g_m q_m. */
101export function strikeGains(model: ModalModel, pl: Plate, p: SimParams): Float64Array {
102 const patch = strikePatch(pl, p);
103 const g = new Float64Array(model.count);
104 for (const k of patch.nodes) {
105 const r = model.wRow[k];
106 if (r < 0) continue;
107 const f = model.masses[k] / patch.mass;
108 for (let m = 0; m < model.count; m++) g[m] += f * model.wShapes[r * model.count + m];
109 }
110 return g;
111}
113/** The state of every mode at the end of the contact phase and the tables recorded during it. */
114export interface ContactResult {
115 /** Time at which the stick is withdrawn, STICK_RELEASE after the contact last opened [s]. */
116 tEnd: number;
117 /** Time at which the contact last opened [s]; every mode decays freely from here on. */
118 tContact: number;
119 /** Complex modal amplitudes A at tEnd: q(t) = Re[A exp(s (t - tEnd))], s = -sigma + i omega_d. */
120 aRe: Float64Array;
121 aIm: Float64Array;
122 /** Modal accelerations sampled at `tableRate` from t = 0 to tEnd, mode-major: qdd[m * tableLen + i]. */
123 qdd: Float32Array;
124 tableLen: number;
125 tableRate: number;
126 /** Modal displacements at the movie frame times before tEnd, mode-major: qFrames[m * nFrames + f]. */
127 qFrames: Float32Array;
128 nFrames: number;
129 /** Stick force history at the fine step [N] and the step [s], for display. */
130 force: Float32Array;
131 dt: number;
132 peakForce: number;
133 contactTime: number;
134}
136/** The contact result restricted to the modes [from, to). */
137export function sliceContact(c: ContactResult, from: number, to: number): ContactResult {
138 return {
139 ...c,
140 aRe: c.aRe.slice(from, to),
141 aIm: c.aIm.slice(from, to),
142 qdd: c.qdd.slice(from * c.tableLen, to * c.tableLen),
143 qFrames: c.qFrames.slice(from * c.nFrames, to * c.nFrames),
144 };
145}
147/**
148 * Steps the stick against the modal system until the contact has been open
149 * for STICK_RELEASE, recording the tables the synthesis needs. Every mode
150 * is advanced by its exact state transition under a constant force over the
151 * step; the stick by symplectic Euler; the force is the Hertz law of the
152 * penetration at the start of the step.
153 */
154export function simulateContact(model: ModalModel, sigma: Float64Array, g: Float64Array, p: SimParams, fps: number, dt = 2e-6, tableRate = 2 * p.fs): ContactResult {
155 const M = model.count;
156 const omega = new Float64Array(M), omegaD = new Float64Array(M);
157 const p11 = new Float64Array(M), p12 = new Float64Array(M), p21 = new Float64Array(M), p22 = new Float64Array(M);
158 for (let m = 0; m < M; m++) {
159 const w = 2 * Math.PI * model.freqs[m];
160 const s = sigma[m];
161 const wd = dampedOmega(model.freqs[m], s);
162 omega[m] = w; omegaD[m] = wd;
163 const e = Math.exp(-s * dt), c = Math.cos(wd * dt), sn = Math.sin(wd * dt);
164 p11[m] = e * (c + (s * sn) / wd);
165 p12[m] = (e * sn) / wd;
166 p21[m] = (-e * w * w * sn) / wd;
167 p22[m] = e * (c - (s * sn) / wd);
168 }
169 const q = new Float64Array(M), qd = new Float64Array(M);
170 const tMax = 0.2;
171 const maxSteps = Math.ceil(tMax / dt);
172 const tableCap = Math.ceil(tMax * tableRate) + 2;
173 const qddTable = new Float32Array(M * tableCap);
174 const frameCap = Math.ceil(tMax * fps) + 2;
175 const qFrames = new Float32Array(M * frameCap);
176 const force = new Float32Array(maxSteps + 1);
177 let tableLen = 0, nFrames = 0, peakForce = 0, contactTime = 0;
178 let z = 0, zd = -p.strikeVelocity; // stick tip, positive away from the plate; it arrives moving toward it
179 let lastOpen = 0, everContact = false, F = 0;
180 let step = 0;
181 const record = (t: number): void => {
182 while (tableLen < tableCap && tableLen / tableRate <= t + 1e-12) {
183 for (let m = 0; m < M; m++) qddTable[m * tableCap + tableLen] = -2 * sigma[m] * qd[m] - omega[m] * omega[m] * q[m] - F * g[m];
184 tableLen++;
185 }
186 while (nFrames < frameCap && nFrames / fps <= t + 1e-12) {
187 for (let m = 0; m < M; m++) qFrames[m * frameCap + nFrames] = q[m];
188 nFrames++;
189 }
190 };
191 for (step = 0; step < maxSteps; step++) {
192 const t = step * dt;
193 let wc = 0;
194 for (let m = 0; m < M; m++) wc += g[m] * q[m];
195 const delta = wc - z; // penetration: the plate has come past the tip
196 F = delta > 0 ? p.contactStiffness * Math.pow(delta, p.contactExponent) : 0;
197 if (F > 0) { everContact = true; lastOpen = t + dt; contactTime += dt; }
198 if (F > peakForce) peakForce = F;
199 force[step] = F;
200 record(t);
201 if (everContact && t - lastOpen >= STICK_RELEASE) break;
202 for (let m = 0; m < M; m++) {
203 const qp = (-F * g[m]) / (omega[m] * omega[m]);
204 const x = q[m] - qp, xd = qd[m];
205 q[m] = p11[m] * x + p12[m] * xd + qp;
206 qd[m] = p21[m] * x + p22[m] * xd;
207 }
208 zd += (F / p.stickMass) * dt;
209 z += zd * dt;
210 }
211 const tEnd = step * dt;
212 const aRe = new Float64Array(M), aIm = new Float64Array(M);
213 for (let m = 0; m < M; m++) {
214 aRe[m] = q[m];
215 aIm[m] = -(qd[m] + sigma[m] * q[m]) / omegaD[m];
216 }
217 const qdd = new Float32Array(M * tableLen);
218 for (let m = 0; m < M; m++) qdd.set(qddTable.subarray(m * tableCap, m * tableCap + tableLen), m * tableLen);
219 const frames = new Float32Array(M * nFrames);
220 for (let m = 0; m < M; m++) frames.set(qFrames.subarray(m * frameCap, m * frameCap + nFrames), m * nFrames);
221 return { tEnd, tContact: Math.min(lastOpen, tEnd), aRe, aIm, qdd, tableLen, tableRate, qFrames: frames, nFrames, force: force.slice(0, step + 1), dt, peakForce, contactTime };
222}
224// ---------- Listener ----------
226/** Listener terms with the common delay removed, so that only the spread of the delays across the plate remains. */
227export function listenerWeights(pl: Plate, p: SimParams): { weight: Float64Array; delay: Float64Array; maxDelay: number } {
228 const { weight, delay } = listenerTerms(pl, p);
229 let dmin = Infinity;
230 for (let k = 0; k < weight.length; k++) if (weight[k] > 0) dmin = Math.min(dmin, delay[k]);
231 let maxDelay = 0;
232 for (let k = 0; k < weight.length; k++) if (weight[k] > 0) { delay[k] -= dmin; maxDelay = Math.max(maxDelay, delay[k]); } else delay[k] = 0;
233 return { weight, delay, maxDelay };
234}
236/**
237 * The time from which the tail formula holds at the listener: every node's
238 * retarded time is past the last contact. The table path covers the output
239 * samples before it.
240 */
241export function tailStart(c: ContactResult, maxDelay: number): number {
242 return Math.min(c.tEnd, c.tContact + maxDelay + 2 / c.tableRate);
243}
245// ---------- Sums over the modes (accept a sliced model) ----------
247/**
248 * The complex radiation gain of each mode for the listener:
249 * G_m = sum_k weight_k phi_m(k) exp(-s_m delay_k), with s_m = -sigma + i omega_d.
250 */
251export function radiationGains(model: ModalModel, sigma: Float64Array, weight: Float64Array, delay: Float64Array): { re: Float64Array; im: Float64Array } {
252 const M = model.count;
253 const re = new Float64Array(M), im = new Float64Array(M);
254 const wd = new Float64Array(M);
255 for (let m = 0; m < M; m++) wd[m] = dampedOmega(model.freqs[m], sigma[m]);
256 for (let r = 0; r < model.nRows; r++) {
257 const k = model.rowNode[r];
258 const wk = weight[k];
259 if (wk === 0) continue;
260 const tau = delay[k];
261 const base = r * M;
262 for (let m = 0; m < M; m++) {
263 const a = wk * model.wShapes[base + m] * Math.exp(sigma[m] * tau);
264 re[m] += a * Math.cos(wd[m] * tau);
265 im[m] -= a * Math.sin(wd[m] * tau);
266 }
267 }
268 return { re, im };
269}
271/** The nodal accelerations of the model's modes over the first `len` samples of the contact table, row-major over the plate nodes: out[r * len + i] += sum_m phi_m(r) qdd_m(i). */
272export function nodalAccelerationTable(model: ModalModel, c: ContactResult, len: number, out?: Float32Array): Float32Array {
273 const L = c.tableLen, M = model.count;
274 const table = out ?? new Float32Array(model.nRows * len);
275 for (let r = 0; r < model.nRows; r++) {
276 const row = r * len, base = r * M;
277 for (let m = 0; m < M; m++) {
278 const phi = model.wShapes[base + m];
279 if (phi === 0) continue;
280 const src = m * L;
281 for (let i = 0; i < len; i++) table[row + i] += phi * c.qdd[src + i];
282 }
283 }
284 return table;
285}
287/**
288 * The listener pressure from `tStart` on: p(t) = sum_m Re[G_m s_m^2 A_m
289 * exp(s_m (t - tEnd))], generated by the two-term recurrence of a damped
290 * sinusoid, in blocks that stay in cache. Output samples before tStart are
291 * left alone.
292 */
293export function modalTail(model: ModalModel, sigma: Float64Array, c: ContactResult, gains: { re: Float64Array; im: Float64Array }, fs: number, tStart: number, nSamples: number, out?: Float32Array): Float32Array {
294 const y = out ?? new Float32Array(nSamples);
295 const dt = 1 / fs;
296 const block = 4096;
297 const M = model.count;
298 const n0 = Math.ceil(tStart * fs - 1e-9);
299 const a = new Float64Array(M), b = new Float64Array(M), p1 = new Float64Array(M), p2 = new Float64Array(M), limit = new Float64Array(M);
300 for (let m = 0; m < M; m++) {
301 const s = sigma[m], wd = dampedOmega(model.freqs[m], s);
302 // c = G s^2 A, then rotated to the first output sample: c exp(s (n0/fs - tEnd)).
303 const s2Re = s * s - wd * wd, s2Im = -2 * s * wd;
304 const saRe = s2Re * c.aRe[m] - s2Im * c.aIm[m], saIm = s2Re * c.aIm[m] + s2Im * c.aRe[m];
305 let cRe = gains.re[m] * saRe - gains.im[m] * saIm, cIm = gains.re[m] * saIm + gains.im[m] * saRe;
306 const t0 = n0 / fs - c.tEnd;
307 const e0 = Math.exp(-s * t0), c0 = Math.cos(wd * t0), s0 = Math.sin(wd * t0);
308 const r = e0 * (cRe * c0 - cIm * s0), i = e0 * (cRe * s0 + cIm * c0);
309 cRe = r; cIm = i;
310 const e = Math.exp(-s * dt), zRe = e * Math.cos(wd * dt), zIm = e * Math.sin(wd * dt);
311 a[m] = 2 * zRe; b[m] = -(e * e);
312 p2[m] = cRe; // p(n0)
313 p1[m] = cRe * zRe - cIm * zIm; // p(n0 + 1)
314 limit[m] = s > 0 ? Math.ceil((Math.log(1e7) / s) * fs) : Infinity; // 140 dB down: dropped
315 }
316 for (let start = n0; start < nSamples; start += block) {
317 const end = Math.min(nSamples, start + block);
318 for (let m = 0; m < M; m++) {
319 if (start - n0 > limit[m]) continue;
320 const am = a[m], bm = b[m];
321 let x1 = p1[m], x2 = p2[m];
322 let n = start;
323 if (start === n0) {
324 y[n] += x2; n++;
325 if (n < end) { y[n] += x1; n++; }
326 }
327 for (; n < end; n++) {
328 const x = am * x1 + bm * x2;
329 y[n] += x;
330 x2 = x1; x1 = x;
331 }
332 p1[m] = x1; p2[m] = x2;
333 }
334 }
335 return y;
336}
338/**
339 * Movie frames of the transverse displacement: frames[f * nGrid + k] +=
340 * sum_m phi_m(k) q_m(t_f), with q from the contact tables before tEnd and
341 * from the free decay after it.
342 */
343export function movieFrames(model: ModalModel, sigma: Float64Array, c: ContactResult, fps: number, nFrames: number, out?: Float32Array): Float32Array {
344 const frames = out ?? new Float32Array(nFrames * model.nGrid);
345 const M = model.count;
346 const q = new Float64Array(M);
347 const wd = new Float64Array(M);
348 for (let m = 0; m < M; m++) wd[m] = dampedOmega(model.freqs[m], sigma[m]);
349 for (let f = 0; f < nFrames; f++) {
350 const tt = f / fps - c.tEnd;
351 for (let m = 0; m < M; m++) {
352 if (f < c.nFrames) { q[m] = c.qFrames[m * c.nFrames + f]; continue; }
353 const e = Math.exp(-sigma[m] * tt);
354 q[m] = e * (c.aRe[m] * Math.cos(wd[m] * tt) - c.aIm[m] * Math.sin(wd[m] * tt));
355 }
356 const base = f * model.nGrid;
357 for (let r = 0; r < model.nRows; r++) {
358 const rb = r * M;
359 let s = 0;
360 for (let m = 0; m < M; m++) s += model.wShapes[rb + m] * q[m];
361 frames[base + model.rowNode[r]] += s;
362 }
363 }
364 return frames;
365}
367// ---------- Assembly of the parts ----------
369/**
370 * The listener pressure over the attack from the nodal acceleration table:
371 * p(t) = sum_k weight_k wdd_k(t - delay_k), linearly interpolated, for the
372 * output samples before tStart.
373 */
374export function attackPressure(model: ModalModel, table: Float32Array, len: number, tableRate: number, weight: Float64Array, delay: Float64Array, fs: number, tStart: number, out: Float32Array): void {
375 const nOut = Math.min(out.length, Math.ceil(tStart * fs - 1e-9));
376 for (let r = 0; r < model.nRows; r++) {
377 const k = model.rowNode[r];
378 const wk = weight[k];
379 if (wk === 0) continue;
380 const row = r * len;
381 const dk = delay[k];
382 for (let n = 0; n < nOut; n++) {
383 const tr = n / fs - dk;
384 if (tr < 0) continue;
385 const x = tr * tableRate;
386 const i = Math.floor(x);
387 if (i + 1 >= len) continue;
388 const fr = x - i;
389 out[n] += wk * ((1 - fr) * table[row + i] + fr * table[row + i + 1]);
390 }
391 }
392}
394export interface Synthesis {
395 audio: Float32Array;
396 frames: Float32Array;
397 nFrames: number;
398 contact: ContactResult;
399}
401/** Runs the whole synthesis on the calling thread (the tests, and browsers without workers). */
402export function synthesize(model: ModalModel, pl: Plate, p: SimParams, withMovie = true): Synthesis {
403 const sigma = modalRates(model);
404 const g = strikeGains(model, pl, p);
405 const contact = simulateContact(model, sigma, g, p, p.fps);
406 const { weight, delay, maxDelay } = listenerWeights(pl, p);
407 const tStart = tailStart(contact, maxDelay);
408 const len = Math.min(contact.tableLen, Math.ceil(tStart * contact.tableRate) + 2);
409 const nSamples = Math.round(p.duration * p.fs);
410 const audio = new Float32Array(nSamples);
411 const table = nodalAccelerationTable(model, contact, len);
412 attackPressure(model, table, len, contact.tableRate, weight, delay, p.fs, tStart, audio);
413 const gains = radiationGains(model, sigma, weight, delay);
414 modalTail(model, sigma, contact, gains, p.fs, tStart, nSamples, audio);
415 const nFrames = withMovie ? Math.round(p.duration * p.fps) : 0;
416 const frames = withMovie ? movieFrames(model, sigma, contact, p.fps, nFrames) : new Float32Array(0);
417 return { audio, frames, nFrames, contact };
418}