5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 1// The modal synthesis: a single mode struck by an impulse gives the closed
2// form damped sinusoid; the radiation gain of a uniform piston equals the
3// analytic value; the attack and the tail join; the sliced parts sum to the
4// whole.
5import { test } from "node:test";
6import assert from "node:assert/strict";
7import { makeGrid, defaultParams, type SimParams, listenerTerms, RHO_AIR, SPEED_OF_SOUND, nodeWeight } from "../src/params.ts";
8import { uniformPlate } from "../src/model/reference.ts";
9import { assemble } from "../src/model/assemble.ts";
10import { computeModes, localRunner } from "../src/modal/modes.ts";
11import { constantLoss, LN1000 } from "../src/loss.ts";
12import {
13 modalModel, modalRates, strikeGains, simulateContact, listenerWeights, radiationGains, tailStart, nodalAccelerationTable, attackPressure, modalTail, movieFrames, synthesize, sliceModel, sliceContact, type ModalModel,
14} from "../src/modal/synth.ts";
16const steel = { E: 200e9, nu: 0.3, rho: 7850 };
18async function smallModel(): Promise<{ model: ModalModel; pl: ReturnType<typeof uniformPlate> }> {
19 const g = makeGrid(0.2, 0.15, 0.0125);
20 const pl = uniformPlate(g, steel, 1e-3);
21 pl.loss = constantLoss(1.0);
22 const sys = assemble(pl);
23 const modes = await computeModes(sys, 4000, localRunner(sys), { sliceSize: 20 });
24 return { model: modalModel(pl, sys, modes), pl };
25}
27test("a single mode under an impulse follows the closed-form damped sinusoid", async () => {
28 const { model, pl } = await smallModel();
29 const p: SimParams = { ...defaultParams, fs: 48000, duration: 0.1, fps: 100 };
30 // Keep only mode 0 by slicing the model.
31 const one = sliceModel(model, 0, 1);
32 const sigma = modalRates(one, pl.loss);
33 const g = strikeGains(one, pl, p);
34 const c = simulateContact(one, sigma, g, p, p.fps);
35 assert.ok(c.peakForce > 0 && c.contactTime > 0);
36 // After tEnd the mode is q(t) = Re[A exp(s (t - tEnd))]; check the movie frames against it.
37 const nFrames = 10;
38 const frames = movieFrames(one, sigma, c, p.fps, nFrames);
39 const w = 2 * Math.PI * one.freqs[0], s = sigma[0], wd = Math.sqrt(w * w - s * s);
40 for (let f = 1; f < nFrames; f++) {
41 const t = f / p.fps - c.tEnd;
42 if (t < 0) continue;
43 const q = Math.exp(-s * t) * (c.aRe[0] * Math.cos(wd * t) - c.aIm[0] * Math.sin(wd * t));
44 // Compare at a node away from the edge.
45 const k = Math.floor(pl.grid.ny / 2) * pl.grid.nx + Math.floor(pl.grid.nx / 3);
46 const phi = one.wShapes[one.wRow[k] * one.count];
47 assert.ok(Math.abs(frames[f * one.nGrid + k] - phi * q) < 1e-9 + 1e-6 * Math.abs(phi * q));
48 }
49 assert.ok(Math.abs(sigma[0] - LN1000 / 1.0) < 1e-9, "decay rate from T60 = 1 s");
50});
52test("the radiation gain of a uniform piston equals the analytic Rayleigh sum", async () => {
53 const { model, pl } = await smallModel();
54 const p: SimParams = { ...defaultParams, listenerDist: 1.0, listenerAngle: 0 };
55 // A fake model with one mode of uniform shape (mass-normalized to 1/sqrt(M)).
56 const uniform: ModalModel = { ...sliceModel(model, 0, 1), freqs: new Float64Array([1000]), padSigma: new Float64Array([0]) };
57 let M = 0;
58 for (let k = 0; k < uniform.nGrid; k++) M += uniform.masses[k];
59 const phi = 1 / Math.sqrt(M);
60 for (let r = 0; r < uniform.nRows; r++) uniform.wShapes[r] = phi;
61 const { weight, delay } = listenerWeights(pl, p);
62 const sigma = new Float64Array([0]);
63 const gains = radiationGains(uniform, sigma, weight, delay);
64 // Analytic: G = sum_k rho dA / (2 pi r_k) phi exp(-i omega tau_k).
65 const lt = listenerTerms(pl, p);
66 let dmin = Infinity;
67 for (let k = 0; k < lt.weight.length; k++) if (lt.weight[k] > 0) dmin = Math.min(dmin, lt.delay[k]);
68 let re = 0, im = 0;
69 const w = 2 * Math.PI * 1000;
70 for (let k = 0; k < lt.weight.length; k++) {
71 if (lt.weight[k] === 0) continue;
72 const tau = lt.delay[k] - dmin;
73 re += lt.weight[k] * phi * Math.cos(w * tau);
74 im -= lt.weight[k] * phi * Math.sin(w * tau);
75 }
76 // The shapes are float32, so agreement is to about 1e-7.
77 assert.ok(Math.abs(gains.re[0] - re) < 1e-6 * Math.abs(re) && Math.abs(gains.im[0] - im) < 1e-6 * Math.max(Math.abs(im), 1e-6));
78 // On axis at 1 m the weights sum to rho A / (2 pi r) within the delay spread.
79 let area = 0;
80 for (let j = 0; j < pl.grid.ny; j++) for (let i = 0; i < pl.grid.nx; i++) area += nodeWeight(pl.grid, i, j) * pl.grid.hx * pl.grid.hy;
81 const approx = (RHO_AIR * area) / (2 * Math.PI * 1.0) * phi;
82 assert.ok(Math.abs(Math.hypot(re, im) / approx - 1) < 0.05);
83 void SPEED_OF_SOUND;
84});
86test("the attack and the tail join, and the sliced parts sum to the whole", async () => {
87 const { model, pl } = await smallModel();
88 const p: SimParams = { ...defaultParams, fs: 48000, duration: 0.2, fps: 120, strikeVelocity: 2 };
89 const whole = synthesize(model, pl, pl.loss, p, true);
90 const nSamples = Math.round(p.duration * p.fs);
91 assert.equal(whole.audio.length, nSamples);
92 // No discontinuity at the tail start beyond the local sample-to-sample variation.
93 const { maxDelay, weight, delay } = listenerWeights(pl, p);
94 const n0 = Math.ceil(tailStart(whole.contact, maxDelay) * p.fs - 1e-9);
95 let maxStep = 0, jump = Math.abs(whole.audio[n0] - whole.audio[n0 - 1]);
96 for (let n = n0 + 2; n < n0 + 200; n++) maxStep = Math.max(maxStep, Math.abs(whole.audio[n] - whole.audio[n - 1]));
97 assert.ok(jump < 3 * maxStep, `jump ${jump} vs typical step ${maxStep}`);
98 // The tail decays: the last 20 ms is quieter than the first 20 ms.
99 const rms = (a: Float32Array) => Math.sqrt(a.reduce((s, v) => s + v * v, 0) / a.length);
100 assert.ok(rms(whole.audio.subarray(nSamples - 960)) < rms(whole.audio.subarray(0, 960)));
101 // Split the modes into three slices and sum the parts.
102 const sigma = modalRates(model, pl.loss);
103 const g = strikeGains(model, pl, p);
104 const c = simulateContact(model, sigma, g, p, p.fps);
105 const tStart = tailStart(c, maxDelay);
106 const len = Math.min(c.tableLen, Math.ceil(tStart * c.tableRate) + 2);
107 const audio = new Float32Array(nSamples);
108 const table = new Float32Array(model.nRows * len);
109 const nFrames = Math.round(p.duration * p.fps);
110 const frames = new Float32Array(nFrames * model.nGrid);
111 const cuts = [0, Math.floor(model.count / 3), Math.floor((2 * model.count) / 3), model.count];
112 for (let s = 0; s < 3; s++) {
113 const sub = sliceModel(model, cuts[s], cuts[s + 1]);
114 const cc = sliceContact(c, cuts[s], cuts[s + 1]);
115 const sg = sigma.slice(cuts[s], cuts[s + 1]);
116 nodalAccelerationTable(sub, cc, len, table);
117 modalTail(sub, sg, cc, radiationGains(sub, sg, weight, delay), p.fs, tStart, nSamples, audio);
118 movieFrames(sub, sg, cc, p.fps, nFrames, frames);
119 }
120 attackPressure(model, table, len, c.tableRate, weight, delay, p.fs, tStart, audio);
121 let err = 0, nrm = 0;
122 for (let n = 0; n < nSamples; n++) { err += (audio[n] - whole.audio[n]) ** 2; nrm += whole.audio[n] ** 2; }
123 assert.ok(Math.sqrt(err / nrm) < 1e-5, `audio split error ${Math.sqrt(err / nrm)}`);
124 err = 0; nrm = 0;
125 for (let i = 0; i < frames.length; i++) { err += (frames[i] - whole.frames[i]) ** 2; nrm += whole.frames[i] ** 2; }
126 assert.ok(Math.sqrt(err / nrm) < 1e-5, `frame split error ${Math.sqrt(err / nrm)}`);
127});
129test("a lighter hit gives a proportionally smaller signal in the linear model", async () => {
130 const { model, pl } = await smallModel();
131 const p1: SimParams = { ...defaultParams, fs: 48000, duration: 0.05, strikeVelocity: 0.2 };
132 const p2: SimParams = { ...p1, strikeVelocity: 0.1 };
133 const a1 = synthesize(model, pl, pl.loss, p1, false).audio, a2 = synthesize(model, pl, pl.loss, p2, false).audio;
134 const peak = (a: Float32Array) => a.reduce((m, v) => Math.max(m, Math.abs(v)), 0);
135 // The Hertz contact is nonlinear (a faster hit is shorter and brighter), so the ratio exceeds 2 somewhat.
136 const ratio = peak(a1) / peak(a2);
137 assert.ok(ratio > 1.8 && ratio < 3.5, `ratio ${ratio}`);
138});