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