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