ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 1/** AudioWorklet processor: a bank of damped modal resonators driven by rigid exciters through a
2 * power-law contact law, with the listener pressure as output.
3 *
4 * Mode m obeys a'' + 2 sigma a' + omega0^2 a = f_m with f_m = sum_p phi_m(p) (-F_p), the modal
5 * force of the contact springs at the patch points p. Its state is a complex number z with
6 * a = Im(z) and dz/dt = lambda z + f_m / omega_d, lambda = -sigma + i omega_d. One sample with
7 * the force held constant is z <- P z + b f_m with P = exp(lambda dt) and
8 * b = (P - 1) / (lambda omega_d), exact for piecewise-constant forcing. The real acceleration is
9 * Im(lambda^2 z) + f_m, so the listener pressure with the complex per-mode gain G_m (Rayleigh
10 * sum, delays folded into a phase at the mode frequency) is
11 * p = sum_m Im(G_m lambda_m^2 z_m) + Re(G_m) f_m .
12 *
13 * Contact is implicit: the force of a sample is the one that satisfies the contact law at the
14 * end of that sample, F_p = k max(w_p + o_p - z_s, 0)^alpha with w_p and z_s advanced by the
15 * force itself. The surface under the patch responds to the step's forces through the P x P
16 * compliance c_pq = sum_m phi_m(p) phi_m(q) Im(b_m), and the exciter through dt^2 / m, so the
17 * end-of-step penetration is affine in F and the small nonlinear system is solved by Newton
18 * iteration. This keeps the bank at one step per sample for any contact stiffness, where an
19 * explicit coupling would need sub-steps. Between hits the bank runs mode-outer over the block.
21 * A hit carries the frame at which its exciter is released, so one processor serves live play
22 * (release at once), timed playback of a recording, and an offline render, where the whole
23 * performance is handed over in processorOptions before rendering begins. A hit whose frame has
24 * not arrived costs nothing until the block it falls in.
27interface InitMessage {
28 type: "init";
29 freq: Float64Array;
30 sigma: Float64Array;
31 gain: Float32Array; // (Nm, 2) complex listener gain, re/im interleaved
32}
34/** New frequencies and decay rates for the same modes, keeping the ringing state: this is what
35 * the design panel sends while damping or material sliders move. */
36interface ParamsMessage {
37 type: "params";
38 freq: Float64Array;
39 sigma: Float64Array;
40}
42interface GainMessage {
43 type: "gain";
44 gain: Float32Array;
45}
47interface HitMessage {
48 type: "hit";
49 id: number;
50 phi: Float32Array; // (P, Nm) point-major mode values at the patch points
51 offsets: Float32Array; // (P,) surface height minus exciter face plane
52 kPoint: number; // contact stiffness per point [N/m^alpha]
53 alpha: number;
54 mass: number;
55 gap: number;
56 velocity: number; // approach speed, positive
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 57 atFrame?: number; // release frame; absent means as soon as possible
58 silent?: boolean; // do not report this hit back (playback and offline renders)
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 61type InMessage = InitMessage | ParamsMessage | GainMessage | HitMessage | { type: "reset" };
63interface Exciter {
64 id: number;
65 phi: Float32Array;
66 offsets: Float32Array;
67 nPoints: number;
68 kPoint: number;
69 alpha: number;
70 mass: number;
71 gap: number;
72 velocity: number;
74 silent: boolean;
76 zs: number;
77 vs: number;
78 startFrame: number;
79 peakForce: number;
80 contactSamples: number;
81 newtonIters: number;
82 touched: boolean;
83 wp: Float64Array; // free end-of-step surface height under the points (with offsets)
84 Fp: Float64Array; // current forces
85 c: Float64Array; // (P, P) one-step compliance of the surface
86 J: Float64Array; // Newton work: (P, P)
87 g: Float64Array;
88 dF: Float64Array;
89}
91const HOLD = 1e3;
92const now: () => number = typeof performance !== "undefined" ? () => performance.now() : () => Date.now();
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 93const RELEASE_CLEARANCE = 2e-3; // an exciter is dropped once it has rebounded this far clear
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 94const NEWTON_MAX = 12;
96class CymbalModalProcessor extends AudioWorkletProcessor {
97 private nModes = 0;
98 private zr = new Float64Array(0);
99 private zi = new Float64Array(0);
100 private lamR = new Float64Array(0);
101 private lamI = new Float64Array(0);
102 private pr = new Float64Array(0); // one-sample pole
103 private pi = new Float64Array(0);
104 private br = new Float64Array(0); // one-sample forcing coefficient
105 private bi = new Float64Array(0);
106 private cr = new Float64Array(0); // output coefficient C = G lambda^2
107 private ci = new Float64Array(0);
108 private dr = new Float64Array(0); // direct term Re(G)
109 private f = new Float64Array(0);
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 110 private lastGain: Float32Array<ArrayBufferLike> = new Float32Array(0);
111 private exciters: Exciter[] = []; // due in this block, or already sounding
112 private pending: Exciter[] = []; // waiting for a later frame
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 113 private loadAcc = 0;
114 private loadBlocks = 0;
115 private lastLoadReport = 0;
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 117 constructor(options?: { processorOptions?: { init?: InitMessage; hits?: HitMessage[] } }) {
119 this.port.onmessage = (ev: MessageEvent<InMessage>) => this.onMessage(ev.data);
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 120 // An offline render hands the whole performance over here, because a port message posted
121 // before startRendering is not guaranteed to arrive before the first block.
122 const opts = options?.processorOptions;
123 if (opts?.init) this.init(opts.init);
124 if (opts?.hits) for (const h of opts.hits) this.addHit(h);
127 private onMessage(msg: InMessage) {
128 switch (msg.type) {
129 case "init":
130 this.init(msg);
131 break;
132 case "params":
133 this.setParams(msg.freq, msg.sigma);
134 break;
135 case "gain":
136 this.setGain(msg.gain);
137 break;
138 case "hit":
139 this.addHit(msg);
140 break;
141 case "reset":
142 this.zr.fill(0);
143 this.zi.fill(0);
144 this.exciters = [];
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 145 this.pending = [];
147 }
148 }
150 private init(msg: InitMessage) {
151 const Nm = msg.freq.length;
152 this.nModes = Nm;
153 this.zr = new Float64Array(Nm);
154 this.zi = new Float64Array(Nm);
155 this.lamR = new Float64Array(Nm);
156 this.lamI = new Float64Array(Nm);
157 this.pr = new Float64Array(Nm);
158 this.pi = new Float64Array(Nm);
159 this.br = new Float64Array(Nm);
160 this.bi = new Float64Array(Nm);
161 this.cr = new Float64Array(Nm);
162 this.ci = new Float64Array(Nm);
163 this.dr = new Float64Array(Nm);
164 this.f = new Float64Array(Nm);
165 this.lastGain = msg.gain;
166 this.setParams(msg.freq, msg.sigma);
167 this.exciters = [];
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 168 this.pending = [];
170 }
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 172 /** Run realistic fake hits through the contact code path so that the JIT has compiled and
173 * optimised it before the first real one; the first contact blocks were otherwise several
174 * times over the deadline. */
175 private warmUp() {
176 const Nm = this.nModes;
177 const scratch = new Float32Array(128);
178 const fake = (P: number, kPoint: number, alpha: number, mass: number, velocity: number) => {
179 const phi = new Float32Array(P * Nm), offsets = new Float32Array(P);
180 for (let p = 0; p < P; p++) {
181 for (let m = 0; m < Nm; m++) phi[p * Nm + m] = 0.7 * Math.sin(0.37 * m + 1.3 * p) * Math.cos(0.11 * m);
182 offsets[p] = 1e-4 * (p - 0.5 * (P - 1));
183 }
184 this.addHit({ type: "hit", id: -1, phi, offsets, kPoint, alpha, mass, gap: 1e-4, velocity, silent: true });
185 };
186 fake(1, 1.8e8, 2.5, 0.3, 8);
187 fake(7, 5e8 / 7, 1.5, 0.03, 7);
188 for (let k = 0; k < 80; k++) {
189 this.renderBlock(scratch, 0);
190 if (k === 40) fake(7, 5e8 / 7, 1.5, 0.03, 5);
191 }
192 this.exciters = [];
193 this.pending = [];
194 this.zr.fill(0);
195 this.zi.fill(0);
196 }
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 198 /** Poles and forcing coefficients of every mode; the ringing state is left alone. */
199 private setParams(freq: Float64Array, sigma: Float64Array) {
200 const Nm = this.nModes;
201 if (freq.length !== Nm) return;
202 const dt = 1 / sampleRate;
203 for (let m = 0; m < Nm; m++) {
204 const w0 = 2 * Math.PI * freq[m];
205 const s = Math.min(sigma[m], 0.999 * w0);
206 const wd = Math.sqrt(Math.max(w0 * w0 - s * s, 1e-6));
207 const lr = -s, li = wd;
208 this.lamR[m] = lr;
209 this.lamI[m] = li;
210 const e = Math.exp(lr * dt);
211 const Pr = e * Math.cos(li * dt), Pi = e * Math.sin(li * dt);
212 this.pr[m] = Pr;
213 this.pi[m] = Pi;
214 // b = (P - 1) / (lambda omega_d)
215 const nr = Pr - 1, ni = Pi;
216 const den = wd * (lr * lr + li * li);
217 this.br[m] = (nr * lr + ni * li) / den;
218 this.bi[m] = (ni * lr - nr * li) / den;
219 }
220 this.setGain(this.lastGain);
221 }
223 private setGain(gain: Float32Array) {
224 const Nm = this.nModes;
225 this.lastGain = gain;
226 if (gain.length < 2 * Nm) return;
227 for (let m = 0; m < Nm; m++) {
228 const gr = gain[2 * m], gi = gain[2 * m + 1];
229 const lr = this.lamR[m], li = this.lamI[m];
230 const l2r = lr * lr - li * li, l2i = 2 * lr * li;
231 this.cr[m] = gr * l2r - gi * l2i;
232 this.ci[m] = gr * l2i + gi * l2r;
233 this.dr[m] = gr;
234 }
235 }
237 private addHit(msg: HitMessage) {
238 const P = msg.offsets.length, Nm = this.nModes;
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 239 if (Nm === 0 || msg.phi.length !== P * Nm) return;
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 240 // one-step compliance c_pq = sum_m phi_m(p) phi_m(q) Im(b_m)
241 const c = new Float64Array(P * P);
242 const phi = msg.phi, bi = this.bi;
243 for (let p = 0; p < P; p++) {
244 for (let q = p; q < P; q++) {
245 let s = 0;
246 const op = p * Nm, oq = q * Nm;
247 for (let m = 0; m < Nm; m++) s += bi[m] * phi[op + m] * phi[oq + m];
248 c[p * P + q] = s;
249 c[q * P + p] = s;
250 }
251 }
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 252 const ex: Exciter = {
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 253 id: msg.id, phi, offsets: msg.offsets, nPoints: P, kPoint: msg.kPoint, alpha: msg.alpha, mass: msg.mass,
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 254 gap: msg.gap, velocity: msg.velocity, atFrame: msg.atFrame ?? 0, silent: msg.silent ?? false,
255 released: false, zs: HOLD, vs: 0, startFrame: 0, peakForce: 0, contactSamples: 0, newtonIters: 0,
256 touched: false, wp: new Float64Array(P), Fp: new Float64Array(P), c,
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 257 J: new Float64Array(P * P), g: new Float64Array(P), dF: new Float64Array(P),
259 if (ex.atFrame > currentFrame + 128) this.pending.push(ex);
260 else this.exciters.push(ex);
263 /** Surface height (with offsets) under the patch points from the current state; returns the max. */
264 private surface(ex: Exciter): number {
265 const Nm = this.nModes, zi = this.zi, phi = ex.phi, P = ex.nPoints, wp = ex.wp;
266 let top = -Infinity;
267 for (let p = 0; p < P; p++) {
268 const o = p * Nm;
269 let s = 0;
270 for (let m = 0; m < Nm; m++) s += phi[o + m] * zi[m];
271 s += ex.offsets[p];
272 wp[p] = s;
273 if (s > top) top = s;
274 }
275 return top;
276 }
278 /** Implicit contact forces for this sample given the free end-of-step surface heights ex.wp. */
279 private solveContact(ex: Exciter, dt: number): number {
280 const P = ex.nPoints, c = ex.c, F = ex.Fp, J = ex.J, g = ex.g, dF = ex.dF;
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 281 const base = ex.zs + dt * ex.vs; // the exciter face at the end of the step, without contact
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 282 const em = dt * dt / ex.mass;
283 const k = ex.kPoint, a = ex.alpha;
284 let touch = false;
285 for (let p = 0; p < P; p++) if (ex.wp[p] > base) touch = true;
286 if (!touch) {
287 F.fill(0);
288 return 0;
289 }
290 for (let it = 0; it < NEWTON_MAX; it++) {
291 ex.newtonIters++;
292 let Ftot = 0;
293 for (let q = 0; q < P; q++) Ftot += F[q];
294 let worst = 0;
295 for (let p = 0; p < P; p++) {
296 let d = ex.wp[p] - base - em * Ftot;
297 for (let q = 0; q < P; q++) d -= c[p * P + q] * F[q];
298 if (d > 0) {
299 const dm1 = Math.pow(d, a - 1);
300 g[p] = F[p] - k * dm1 * d;
301 const s = k * a * dm1;
302 for (let q = 0; q < P; q++) J[p * P + q] = (p === q ? 1 : 0) + s * (c[p * P + q] + em);
303 } else {
304 g[p] = F[p];
305 for (let q = 0; q < P; q++) J[p * P + q] = p === q ? 1 : 0;
306 }
307 }
308 solveSmall(J, g, dF, P);
309 for (let p = 0; p < P; p++) {
310 F[p] = Math.max(0, F[p] - dF[p]);
311 worst = Math.max(worst, Math.abs(dF[p]));
312 }
313 let scale = 1;
314 for (let p = 0; p < P; p++) scale = Math.max(scale, F[p]);
315 if (worst < 1e-7 * scale) break;
316 }
317 let Ftot = 0;
318 for (let p = 0; p < P; p++) Ftot += F[p];
319 return Ftot;
320 }
322 process(_inputs: Float32Array[][], outputs: Float32Array[][]): boolean {
323 const out = outputs[0][0];
324 if (!out || this.nModes === 0) return true;
325 const t0 = now();
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 326 this.renderBlock(out, currentFrame);
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 327 const cost = now() - t0;
328 this.loadAcc += cost / (1000 * out.length / sampleRate);
329 this.loadBlocks++;
330 if (currentFrame - this.lastLoadReport >= sampleRate / 4) {
331 this.port.postMessage({ type: "load", mean: this.loadAcc / this.loadBlocks, frame: currentFrame });
332 this.loadAcc = 0;
333 this.loadBlocks = 0;
334 this.lastLoadReport = currentFrame;
335 }
336 return true;
337 }
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 339 private renderBlock(out: Float32Array, frame0: number): void {
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 340 const B = out.length;
341 const Nm = this.nModes;
342 const zr = this.zr, zi = this.zi, cr = this.cr, ci = this.ci, pr = this.pr, pi = this.pi;
343 const dt = 1 / sampleRate;
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 345 // anything due in this block joins the working set; the rest costs nothing yet
346 if (this.pending.length) {
347 const end = frame0 + B;
348 const still: Exciter[] = [];
349 for (const ex of this.pending) {
350 if (ex.atFrame < end) this.exciters.push(ex);
351 else still.push(ex);
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 353 this.pending = still;
356 if (this.exciters.length === 0) {
357 // free ringing: mode-outer over the block
358 out.fill(0);
359 for (let m = 0; m < Nm; m++) {
360 let r = zr[m], i = zi[m];
361 const a = pr[m], b = pi[m], c = cr[m], d = ci[m];
362 for (let n = 0; n < B; n++) {
363 const r2 = a * r - b * i;
364 i = a * i + b * r;
365 r = r2;
366 out[n] += c * i + d * r;
367 }
368 zr[m] = r;
369 zi[m] = i;
370 }
371 } else {
372 const f = this.f, dr = this.dr, br = this.br, bi = this.bi;
373 const exs = this.exciters;
374 for (let n = 0; n < B; n++) {
375 // free step of the bank
376 for (let m = 0; m < Nm; m++) {
377 const r = zr[m], i = zi[m];
378 zr[m] = pr[m] * r - pi[m] * i;
379 zi[m] = pr[m] * i + pi[m] * r;
380 }
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 381 // release whatever is due at this sample, just above the surface as it is now
382 const frame = frame0 + n;
383 for (let e = 0; e < exs.length; e++) {
384 const ex = exs[e];
385 if (ex.released || frame < ex.atFrame) continue;
386 ex.zs = this.surface(ex) + ex.gap;
387 ex.vs = -ex.velocity;
388 ex.released = true;
389 ex.startFrame = frame;
390 if (!ex.silent) this.port.postMessage({ type: "hitStarted", id: ex.id, frame, time: frame / sampleRate });
391 }
392 // implicit contact forces, then the exciter update
394 for (let e = 0; e < exs.length; e++) {
395 const ex = exs[e];
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 396 if (!ex.released) continue;
398 const F = this.solveContact(ex, dt);
399 if (F > 0) {
400 any = true;
401 ex.touched = true;
402 ex.contactSamples++;
403 if (F > ex.peakForce) ex.peakForce = F;
404 }
405 ex.vs += dt * F / ex.mass;
406 ex.zs += dt * ex.vs;
407 }
408 // modal forces, forced part of the step, output
409 let y = 0;
410 if (any) {
411 f.fill(0);
412 for (let e = 0; e < exs.length; e++) {
413 const ex = exs[e], P = ex.nPoints, phi = ex.phi, F = ex.Fp;
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 414 if (!ex.released) continue;
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 415 for (let p = 0; p < P; p++) {
416 const Fp = F[p];
417 if (Fp === 0) continue;
418 const o = p * Nm;
419 for (let m = 0; m < Nm; m++) f[m] -= Fp * phi[o + m];
420 }
421 }
422 for (let m = 0; m < Nm; m++) {
423 const fm = f[m];
424 const r = zr[m] + br[m] * fm, i = zi[m] + bi[m] * fm;
425 zr[m] = r;
426 zi[m] = i;
427 y += cr[m] * i + ci[m] * r + dr[m] * fm;
428 }
429 } else {
430 for (let m = 0; m < Nm; m++) y += cr[m] * zi[m] + ci[m] * zr[m];
431 }
432 out[n] = y;
433 }
434 // drop exciters that have rebounded clear of the surface
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 435 const frameEnd = frame0 + B;
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 436 this.exciters = exs.filter((ex) => {
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 437 if (!ex.released) return true;
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 438 const top = this.surface(ex);
439 const done = ex.touched && ex.vs > 0 && ex.zs - top > RELEASE_CLEARANCE;
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 440 const lost = !ex.touched && frameEnd - ex.startFrame > sampleRate * 0.1;
ae4a407cymplay: design and play a cymsim cymbal in the browserJeremy Magland 441 if (done || lost) {
442 if (!ex.silent) this.port.postMessage({
9c9fa54Guide the first visit, and separate recording from playingJeremy Magland 443 type: "hitDone", id: ex.id, peakForce: ex.peakForce, contactTime: ex.contactSamples * dt,
444 rebound: ex.vs, newtonIters: ex.newtonIters, contactSamples: ex.contactSamples,
446 return false;
447 }
448 return true;
449 });
450 }
451 }
452}
454/** Solve A x = b for a small dense system (row-major A, n <= 16) by Gaussian elimination with
455 * partial pivoting; A and b are destroyed. */
456function solveSmall(A: Float64Array, b: Float64Array, x: Float64Array, n: number): void {
457 for (let c = 0; c < n; c++) {
458 let piv = c, best = Math.abs(A[c * n + c]);
459 for (let r = c + 1; r < n; r++) {
460 const v = Math.abs(A[r * n + c]);
461 if (v > best) { best = v; piv = r; }
462 }
463 if (piv !== c) {
464 for (let k = 0; k < n; k++) { const t = A[c * n + k]; A[c * n + k] = A[piv * n + k]; A[piv * n + k] = t; }
465 const t = b[c]; b[c] = b[piv]; b[piv] = t;
466 }
467 const d = A[c * n + c] || 1e-300;
468 for (let r = c + 1; r < n; r++) {
469 const fct = A[r * n + c] / d;
470 if (fct === 0) continue;
471 for (let k = c; k < n; k++) A[r * n + k] -= fct * A[c * n + k];
472 b[r] -= fct * b[c];
473 }
474 }
475 for (let r = n - 1; r >= 0; r--) {
476 let s = b[r];
477 for (let k = r + 1; k < n; k++) s -= A[r * n + k] * x[k];
478 x[r] = s / (A[r * n + r] || 1e-300);
479 }
480}
482registerProcessor("cymbal-modal", CymbalModalProcessor);