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.
20 */
22interface InitMessage {
23 type: "init";
24 freq: Float64Array;
25 sigma: Float64Array;
26 gain: Float32Array; // (Nm, 2) complex listener gain, re/im interleaved
27}
29/** New frequencies and decay rates for the same modes, keeping the ringing state: this is what
30 * the design panel sends while damping or material sliders move. */
31interface ParamsMessage {
32 type: "params";
33 freq: Float64Array;
34 sigma: Float64Array;
35}
37interface GainMessage {
38 type: "gain";
39 gain: Float32Array;
40}
42interface HitMessage {
43 type: "hit";
44 id: number;
45 phi: Float32Array; // (P, Nm) point-major mode values at the patch points
46 offsets: Float32Array; // (P,) surface height minus exciter face plane
47 kPoint: number; // contact stiffness per point [N/m^alpha]
48 alpha: number;
49 mass: number;
50 gap: number;
51 velocity: number; // approach speed, positive
52}
54type InMessage = InitMessage | ParamsMessage | GainMessage | HitMessage | { type: "reset" } | { type: "capture"; on: boolean };
56interface Exciter {
57 id: number;
58 phi: Float32Array;
59 offsets: Float32Array;
60 nPoints: number;
61 kPoint: number;
62 alpha: number;
63 mass: number;
64 gap: number;
65 velocity: number;
66 released: boolean;
67 zs: number;
68 vs: number;
69 startFrame: number;
70 peakForce: number;
71 contactSamples: number;
72 newtonIters: number;
73 touched: boolean;
74 silent: boolean; // warm-up exciter: no messages
75 wp: Float64Array; // free end-of-step surface height under the points (with offsets)
76 Fp: Float64Array; // current forces
77 c: Float64Array; // (P, P) one-step compliance of the surface
78 J: Float64Array; // Newton work: (P, P)
79 g: Float64Array;
80 dF: Float64Array;
81}
83const HOLD = 1e3;
84const now: () => number = typeof performance !== "undefined" ? () => performance.now() : () => Date.now();
85const RELEASE_CLEARANCE = 2e-3; // exciter is dropped once it has rebounded this far above the surface
86const NEWTON_MAX = 12;
88class CymbalModalProcessor extends AudioWorkletProcessor {
89 private nModes = 0;
90 private zr = new Float64Array(0);
91 private zi = new Float64Array(0);
92 private lamR = new Float64Array(0);
93 private lamI = new Float64Array(0);
94 private pr = new Float64Array(0); // one-sample pole
95 private pi = new Float64Array(0);
96 private br = new Float64Array(0); // one-sample forcing coefficient
97 private bi = new Float64Array(0);
98 private cr = new Float64Array(0); // output coefficient C = G lambda^2
99 private ci = new Float64Array(0);
100 private dr = new Float64Array(0); // direct term Re(G)
101 private f = new Float64Array(0);
102 private exciters: Exciter[] = [];
103 private capture = false;
104 private captureBuf: Float32Array[] = [];
105 private loadAcc = 0;
106 private loadBlocks = 0;
107 private lastLoadReport = 0;
108 private lastGain: Float32Array<ArrayBufferLike> = new Float32Array(0);
110 constructor() {
111 super();
112 this.port.onmessage = (ev: MessageEvent<InMessage>) => this.onMessage(ev.data);
113 }
115 private onMessage(msg: InMessage) {
116 switch (msg.type) {
117 case "init":
118 this.init(msg);
119 break;
120 case "params":
121 this.setParams(msg.freq, msg.sigma);
122 break;
123 case "gain":
124 this.setGain(msg.gain);
125 break;
126 case "hit":
127 this.addHit(msg);
128 break;
129 case "reset":
130 this.zr.fill(0);
131 this.zi.fill(0);
132 this.exciters = [];
133 break;
134 case "capture":
135 this.capture = msg.on;
136 if (!msg.on) {
137 this.port.postMessage({ type: "captured", chunks: this.captureBuf });
138 this.captureBuf = [];
139 }
140 break;
141 }
142 }
144 private init(msg: InitMessage) {
145 const Nm = msg.freq.length;
146 this.nModes = Nm;
147 this.zr = new Float64Array(Nm);
148 this.zi = new Float64Array(Nm);
149 this.lamR = new Float64Array(Nm);
150 this.lamI = new Float64Array(Nm);
151 this.pr = new Float64Array(Nm);
152 this.pi = new Float64Array(Nm);
153 this.br = new Float64Array(Nm);
154 this.bi = new Float64Array(Nm);
155 this.cr = new Float64Array(Nm);
156 this.ci = new Float64Array(Nm);
157 this.dr = new Float64Array(Nm);
158 this.f = new Float64Array(Nm);
159 this.lastGain = msg.gain;
160 this.setParams(msg.freq, msg.sigma);
161 this.exciters = [];
162 this.warmUp();
163 }
165 /** Poles and forcing coefficients of every mode; the ringing state is left alone. */
166 private setParams(freq: Float64Array, sigma: Float64Array) {
167 const Nm = this.nModes;
168 if (freq.length !== Nm) return;
169 const dt = 1 / sampleRate;
170 for (let m = 0; m < Nm; m++) {
171 const w0 = 2 * Math.PI * freq[m];
172 const s = Math.min(sigma[m], 0.999 * w0);
173 const wd = Math.sqrt(Math.max(w0 * w0 - s * s, 1e-6));
174 const lr = -s, li = wd;
175 this.lamR[m] = lr;
176 this.lamI[m] = li;
177 const e = Math.exp(lr * dt);
178 const Pr = e * Math.cos(li * dt), Pi = e * Math.sin(li * dt);
179 this.pr[m] = Pr;
180 this.pi[m] = Pi;
181 // b = (P - 1) / (lambda omega_d)
182 const nr = Pr - 1, ni = Pi;
183 const den = wd * (lr * lr + li * li);
184 this.br[m] = (nr * lr + ni * li) / den;
185 this.bi[m] = (ni * lr - nr * li) / den;
186 }
187 this.setGain(this.lastGain);
188 }
190 /** Run realistic fake hits through the contact code path so that the JIT has compiled and
191 * optimised it before the first real hit (the first contact blocks were otherwise several
192 * times over the deadline). Two hits (a one-point soft one and a seven-point stiff one with
193 * staggered offsets) contact, rebound and are dropped, then the bank rings freely. */
194 private warmUp() {
195 const Nm = this.nModes;
196 const scratch = new Float32Array(128);
197 const fake = (P: number, kPoint: number, alpha: number, mass: number, velocity: number) => {
198 const phi = new Float32Array(P * Nm), offsets = new Float32Array(P);
199 for (let p = 0; p < P; p++) {
200 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);
201 offsets[p] = 1e-4 * (p - 0.5 * (P - 1));
202 }
203 this.addHit({ type: "hit", id: -1, phi, offsets, kPoint, alpha, mass, gap: 1e-4, velocity });
204 this.exciters[this.exciters.length - 1].silent = true;
205 };
206 fake(1, 1.8e8, 2.5, 0.3, 8);
207 fake(7, 5e8 / 7, 1.5, 0.03, 7);
208 for (let k = 0; k < 80; k++) {
209 this.renderBlock(scratch);
210 if (k === 40) fake(7, 5e8 / 7, 1.5, 0.03, 5);
211 }
212 this.exciters = [];
213 this.zr.fill(0);
214 this.zi.fill(0);
215 }
217 private setGain(gain: Float32Array) {
218 const Nm = this.nModes;
219 this.lastGain = gain;
220 if (gain.length < 2 * Nm) return;
221 for (let m = 0; m < Nm; m++) {
222 const gr = gain[2 * m], gi = gain[2 * m + 1];
223 const lr = this.lamR[m], li = this.lamI[m];
224 const l2r = lr * lr - li * li, l2i = 2 * lr * li;
225 this.cr[m] = gr * l2r - gi * l2i;
226 this.ci[m] = gr * l2i + gi * l2r;
227 this.dr[m] = gr;
228 }
229 }
231 private addHit(msg: HitMessage) {
232 const P = msg.offsets.length, Nm = this.nModes;
233 if (msg.phi.length !== P * Nm) return;
234 // one-step compliance c_pq = sum_m phi_m(p) phi_m(q) Im(b_m)
235 const c = new Float64Array(P * P);
236 const phi = msg.phi, bi = this.bi;
237 for (let p = 0; p < P; p++) {
238 for (let q = p; q < P; q++) {
239 let s = 0;
240 const op = p * Nm, oq = q * Nm;
241 for (let m = 0; m < Nm; m++) s += bi[m] * phi[op + m] * phi[oq + m];
242 c[p * P + q] = s;
243 c[q * P + p] = s;
244 }
245 }
246 this.exciters.push({
247 id: msg.id, phi, offsets: msg.offsets, nPoints: P, kPoint: msg.kPoint, alpha: msg.alpha, mass: msg.mass,
248 gap: msg.gap, velocity: msg.velocity, released: false, zs: HOLD, vs: 0, startFrame: 0, peakForce: 0,
249 contactSamples: 0, newtonIters: 0, touched: false, silent: false, wp: new Float64Array(P), Fp: new Float64Array(P), c,
250 J: new Float64Array(P * P), g: new Float64Array(P), dF: new Float64Array(P),
251 });
252 }
254 /** Surface height (with offsets) under the patch points from the current state; returns the max. */
255 private surface(ex: Exciter): number {
256 const Nm = this.nModes, zi = this.zi, phi = ex.phi, P = ex.nPoints, wp = ex.wp;
257 let top = -Infinity;
258 for (let p = 0; p < P; p++) {
259 const o = p * Nm;
260 let s = 0;
261 for (let m = 0; m < Nm; m++) s += phi[o + m] * zi[m];
262 s += ex.offsets[p];
263 wp[p] = s;
264 if (s > top) top = s;
265 }
266 return top;
267 }
269 /** Implicit contact forces for this sample given the free end-of-step surface heights ex.wp. */
270 private solveContact(ex: Exciter, dt: number): number {
271 const P = ex.nPoints, c = ex.c, F = ex.Fp, J = ex.J, g = ex.g, dF = ex.dF;
272 const base = ex.zs + dt * ex.vs; // exciter face at the end of the step without contact
273 const em = dt * dt / ex.mass;
274 const k = ex.kPoint, a = ex.alpha;
275 // quick exit: no point can touch even with zero force
276 let touch = false;
277 for (let p = 0; p < P; p++) if (ex.wp[p] > base) touch = true;
278 if (!touch) {
279 F.fill(0);
280 return 0;
281 }
282 for (let it = 0; it < NEWTON_MAX; it++) {
283 ex.newtonIters++;
284 let Ftot = 0;
285 for (let q = 0; q < P; q++) Ftot += F[q];
286 let worst = 0;
287 for (let p = 0; p < P; p++) {
288 let d = ex.wp[p] - base - em * Ftot;
289 for (let q = 0; q < P; q++) d -= c[p * P + q] * F[q];
290 if (d > 0) {
291 const dm1 = Math.pow(d, a - 1);
292 g[p] = F[p] - k * dm1 * d;
293 const s = k * a * dm1;
294 for (let q = 0; q < P; q++) J[p * P + q] = (p === q ? 1 : 0) + s * (c[p * P + q] + em);
295 } else {
296 g[p] = F[p];
297 for (let q = 0; q < P; q++) J[p * P + q] = p === q ? 1 : 0;
298 }
299 }
300 solveSmall(J, g, dF, P);
301 for (let p = 0; p < P; p++) {
302 F[p] = Math.max(0, F[p] - dF[p]);
303 worst = Math.max(worst, Math.abs(dF[p]));
304 }
305 let scale = 1;
306 for (let p = 0; p < P; p++) scale = Math.max(scale, F[p]);
307 if (worst < 1e-7 * scale) break;
308 }
309 let Ftot = 0;
310 for (let p = 0; p < P; p++) Ftot += F[p];
311 return Ftot;
312 }
314 process(_inputs: Float32Array[][], outputs: Float32Array[][]): boolean {
315 const out = outputs[0][0];
316 if (!out || this.nModes === 0) return true;
317 const t0 = now();
318 this.renderBlock(out);
319 if (this.capture) this.captureBuf.push(Float32Array.from(out));
320 const cost = now() - t0;
321 this.loadAcc += cost / (1000 * out.length / sampleRate);
322 this.loadBlocks++;
323 if (currentFrame - this.lastLoadReport >= sampleRate / 4) {
324 this.port.postMessage({ type: "load", mean: this.loadAcc / this.loadBlocks, frame: currentFrame });
325 this.loadAcc = 0;
326 this.loadBlocks = 0;
327 this.lastLoadReport = currentFrame;
328 }
329 return true;
330 }
332 private renderBlock(out: Float32Array): void {
333 const B = out.length;
334 const Nm = this.nModes;
335 const zr = this.zr, zi = this.zi, cr = this.cr, ci = this.ci, pr = this.pr, pi = this.pi;
336 const dt = 1 / sampleRate;
338 // release any exciter that is waiting: place it just above the current surface
339 for (const ex of this.exciters) {
340 if (!ex.released) {
341 ex.zs = this.surface(ex) + ex.gap;
342 ex.vs = -ex.velocity;
343 ex.released = true;
344 ex.startFrame = currentFrame;
345 if (!ex.silent) this.port.postMessage({ type: "hitStarted", id: ex.id, frame: currentFrame, time: currentFrame / sampleRate });
346 }
347 }
349 if (this.exciters.length === 0) {
350 // free ringing: mode-outer over the block
351 out.fill(0);
352 for (let m = 0; m < Nm; m++) {
353 let r = zr[m], i = zi[m];
354 const a = pr[m], b = pi[m], c = cr[m], d = ci[m];
355 for (let n = 0; n < B; n++) {
356 const r2 = a * r - b * i;
357 i = a * i + b * r;
358 r = r2;
359 out[n] += c * i + d * r;
360 }
361 zr[m] = r;
362 zi[m] = i;
363 }
364 } else {
365 const f = this.f, dr = this.dr, br = this.br, bi = this.bi;
366 const exs = this.exciters;
367 for (let n = 0; n < B; n++) {
368 // free step of the bank
369 for (let m = 0; m < Nm; m++) {
370 const r = zr[m], i = zi[m];
371 zr[m] = pr[m] * r - pi[m] * i;
372 zi[m] = pr[m] * i + pi[m] * r;
373 }
374 // free end-of-step surface under every patch, implicit contact forces, exciter update
375 let any = false;
376 for (let e = 0; e < exs.length; e++) {
377 const ex = exs[e];
378 this.surface(ex);
379 const F = this.solveContact(ex, dt);
380 if (F > 0) {
381 any = true;
382 ex.touched = true;
383 ex.contactSamples++;
384 if (F > ex.peakForce) ex.peakForce = F;
385 }
386 ex.vs += dt * F / ex.mass;
387 ex.zs += dt * ex.vs;
388 }
389 // modal forces, forced part of the step, output
390 let y = 0;
391 if (any) {
392 f.fill(0);
393 for (let e = 0; e < exs.length; e++) {
394 const ex = exs[e], P = ex.nPoints, phi = ex.phi, F = ex.Fp;
395 for (let p = 0; p < P; p++) {
396 const Fp = F[p];
397 if (Fp === 0) continue;
398 const o = p * Nm;
399 for (let m = 0; m < Nm; m++) f[m] -= Fp * phi[o + m];
400 }
401 }
402 for (let m = 0; m < Nm; m++) {
403 const fm = f[m];
404 const r = zr[m] + br[m] * fm, i = zi[m] + bi[m] * fm;
405 zr[m] = r;
406 zi[m] = i;
407 y += cr[m] * i + ci[m] * r + dr[m] * fm;
408 }
409 } else {
410 for (let m = 0; m < Nm; m++) y += cr[m] * zi[m] + ci[m] * zr[m];
411 }
412 out[n] = y;
413 }
414 // drop exciters that have rebounded clear of the surface
415 this.exciters = exs.filter((ex) => {
416 const top = this.surface(ex);
417 const done = ex.touched && ex.vs > 0 && ex.zs - top > RELEASE_CLEARANCE;
418 const lost = !ex.touched && currentFrame - ex.startFrame > sampleRate * 0.1;
419 if (done || lost) {
420 if (!ex.silent) this.port.postMessage({
421 type: "hitDone", id: ex.id, peakForce: ex.peakForce, contactTime: ex.contactSamples * dt, rebound: ex.vs,
422 newtonIters: ex.newtonIters, contactSamples: ex.contactSamples,
423 });
424 return false;
425 }
426 return true;
427 });
428 }
429 }
430}
432/** Solve A x = b for a small dense system (row-major A, n <= 16) by Gaussian elimination with
433 * partial pivoting; A and b are destroyed. */
434function solveSmall(A: Float64Array, b: Float64Array, x: Float64Array, n: number): void {
435 for (let c = 0; c < n; c++) {
436 let piv = c, best = Math.abs(A[c * n + c]);
437 for (let r = c + 1; r < n; r++) {
438 const v = Math.abs(A[r * n + c]);
439 if (v > best) { best = v; piv = r; }
440 }
441 if (piv !== c) {
442 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; }
443 const t = b[c]; b[c] = b[piv]; b[piv] = t;
444 }
445 const d = A[c * n + c] || 1e-300;
446 for (let r = c + 1; r < n; r++) {
447 const fct = A[r * n + c] / d;
448 if (fct === 0) continue;
449 for (let k = c; k < n; k++) A[r * n + k] -= fct * A[c * n + k];
450 b[r] -= fct * b[c];
451 }
452 }
453 for (let r = n - 1; r >= 0; r--) {
454 let s = b[r];
455 for (let k = r + 1; k < n; k++) s -= A[r * n + k] * x[k];
456 x[r] = s / (A[r * n + r] || 1e-300);
457 }
458}
460registerProcessor("cymbal-modal", CymbalModalProcessor);