Add instrument/: the cymbal behind the track-render interface
The instrument/ module implements the whole-track instrument interface
used by midiboot and now by the web app (one render path; the webapp
worker is gone, and gpuSim.ts moves in here so both share the WGSL
time loop). The loop runs on WebGPU when available, with per-chunk
progress and cooperative cancellation, and falls back to the numbl
interpreter; the preview tail is chosen per backend (4 s GPU, 1.5 s
CPU) and truncated rings get a short end fade. Finished strokes
persist in IndexedDB with LRU eviction and a 30-day age-out. Replaces
the exported bundle; export-src/ and exports/ are removed.
24 changed files+3882−912
.github/workflows/deploy.ymlmodified+4−0View file
@@ -39,6 +39,10 @@ jobs:
3939 with:
4040 node-version: 22
4141
42+ - name: Install (instrument)
43+ working-directory: instrument
44+ run: npm ci
45+
4246 - name: Install
4347 working-directory: webapp
4448 run: npm ci
README.mdmodified+1−0View file
@@ -9,6 +9,7 @@ The sibling [cymbal](../cymbal) repository models the same instrument by an enti
99 | Path | What it is |
1010 |---|---|
1111 | `matlab/` | The algorithm in MATLAB syntax, run by [numbl](https://numbl.org) |
12+| `instrument/` | The model behind the track-render instrument interface used by [midiboot](../../proof-of-concept/midiboot) and by the web app (worker, note mapping, caching, the WGSL time loop; `npm test` there runs the CPU path in node) |
1213 | `webapp/` | Interactive browser app that executes the MATLAB code client-side; its "How it works" button opens an interactive tutorial on the theory and the solver, with figures and in-page computed audio demos |
1314
1415 The MATLAB code runs from the command line with `cd matlab && numbl run cymbal2_demo.m`, or `cymbal2_validate.m` to run the checks. The web app is a Vite project; see `webapp/README.md`.
export-src/cymbal2/export.jsondeleted+0−23View file
@@ -1,23 +0,0 @@
1-{
2- "id": "cymbal2",
3- "name": "Cymbal2",
4- "description": "Time-domain nonlinear shell cymbal (Samejima 2021); pitch scales the shell around MIDI 60 (48..72), velocity onto stroke speed. Renders are slow: the crash is computed, not modeled.",
5- "source": "simulated_instruments/cymbal2",
6- "knobs": [
7- {
8- "key": "strike_r",
9- "label": "strike radius",
10- "min": 0.05,
11- "max": 0.99,
12- "default": 0.6
13- },
14- {
15- "key": "grip_spring",
16- "label": "grip stiffness (N m/rad)",
17- "min": 20,
18- "max": 2000,
19- "default": 140,
20- "scale": "log"
21- }
22- ]
23-}
export-src/cymbal2/render_note.mdeleted+0−37View file
@@ -1,37 +0,0 @@
1-function y = render_note(note)
2-% MIDI pitch scales the shell geometry: every length by
3-% s = 2^-((pitch - 60)/12), which shifts the whole spectrum by the
4-% semitone ratio exactly (shell modes go as thickness / radius^2, and the
5-% curvature and washer frequencies scale the same way when the washer
6-% circle scales too). The mallet stays fixed, so this is one player
7-% striking a family of cymbal sizes; pitch 60 is the default 18-inch
8-% shell, playable 48..72. The sound is inharmonic, so the note number is
9-% a handle on register, not a tuned fundamental. Velocity drives the hand
10-% speed of the stroke, and the crash character follows from the physics:
11-% soft notes stay nearly linear, hard notes wash. Fixed: a single stroke,
12-% a fixed seed, normalize off so velocity survives. The preview tier
13-% coarsens the azimuthal grid and shortens the tail; the full tier
14-% carries the ring to the end and takes minutes per distinct note.
15-P = cymbal2_defaults();
16-P.sample_rate = note.sample_rate;
17-P.strike_count = 1;
18-P.seed = 1;
19-P.pre_delay_ms = 0;
20-P.normalize = 0;
21-% The preview tier drops the oversampling and shortens the tail but keeps
22-% the azimuthal grid: hard strokes need the full wavenumber band to carry
23-% the cascade, and coarsening it is what destabilizes the integration.
24-if note.preview ~= 0
25- P.oversample = 1;
26- P.tail_s = 1.5;
27-end
28-s = 2 ^ (-(max(48, min(72, note.pitch)) - 60) / 12);
29-P.radius_m = P.radius_m * s;
30-P.rise_m = P.rise_m * s;
31-P.thickness_m = P.thickness_m * s;
32-P.washer_radius_m = P.washer_radius_m * s;
33-P.strike_velocity = 0.3 + (6.5 - 0.3) * note.velocity;
34-P.strike_r = note.params.strike_r;
35-P.grip_spring = note.params.grip_spring;
36-y = cymbal2_render(P);
37-end
exports/cymbal2.jsondeleted+0−74View file
@@ -1,74 +0,0 @@
1-{
2- "format": "numbl-instrument-v1",
3- "id": "cymbal2",
4- "name": "Cymbal2",
5- "description": "Time-domain nonlinear shell cymbal (Samejima 2021); pitch scales the shell around MIDI 60 (48..72), velocity onto stroke speed. Renders are slow: the crash is computed, not modeled.",
6- "source": "simulated_instruments/cymbal2",
7- "knobs": [
8- {
9- "key": "strike_r",
10- "label": "strike radius",
11- "min": 0.05,
12- "max": 0.99,
13- "default": 0.6
14- },
15- {
16- "key": "grip_spring",
17- "label": "grip stiffness (N m/rad)",
18- "min": 20,
19- "max": 2000,
20- "default": 140,
21- "scale": "log"
22- }
23- ],
24- "files": [
25- {
26- "name": "render_note.m",
27- "content": "function y = render_note(note)\n% MIDI pitch scales the shell geometry: every length by\n% s = 2^-((pitch - 60)/12), which shifts the whole spectrum by the\n% semitone ratio exactly (shell modes go as thickness / radius^2, and the\n% curvature and washer frequencies scale the same way when the washer\n% circle scales too). The mallet stays fixed, so this is one player\n% striking a family of cymbal sizes; pitch 60 is the default 18-inch\n% shell, playable 48..72. The sound is inharmonic, so the note number is\n% a handle on register, not a tuned fundamental. Velocity drives the hand\n% speed of the stroke, and the crash character follows from the physics:\n% soft notes stay nearly linear, hard notes wash. Fixed: a single stroke,\n% a fixed seed, normalize off so velocity survives. The preview tier\n% coarsens the azimuthal grid and shortens the tail; the full tier\n% carries the ring to the end and takes minutes per distinct note.\nP = cymbal2_defaults();\nP.sample_rate = note.sample_rate;\nP.strike_count = 1;\nP.seed = 1;\nP.pre_delay_ms = 0;\nP.normalize = 0;\n% The preview tier drops the oversampling and shortens the tail but keeps\n% the azimuthal grid: hard strokes need the full wavenumber band to carry\n% the cascade, and coarsening it is what destabilizes the integration.\nif note.preview ~= 0\n P.oversample = 1;\n P.tail_s = 1.5;\nend\ns = 2 ^ (-(max(48, min(72, note.pitch)) - 60) / 12);\nP.radius_m = P.radius_m * s;\nP.rise_m = P.rise_m * s;\nP.thickness_m = P.thickness_m * s;\nP.washer_radius_m = P.washer_radius_m * s;\nP.strike_velocity = 0.3 + (6.5 - 0.3) * note.velocity;\nP.strike_r = note.params.strike_r;\nP.grip_spring = note.params.grip_spring;\ny = cymbal2_render(P);\nend\n"
28- },
29- {
30- "name": "cymbal2_defaults.m",
31- "content": "function P = cymbal2_defaults()\n% CYMBAL2_DEFAULTS Default parameters, following Samejima (2021) Table 1:\n% an 18-inch medium crash (45 cm diameter) modeled as a shallow spherical\n% shell, supported at the center by a washer of finite impedance and struck\n% by a mallet whose stick is an elastic beam held in a compliant grip.\n% The output duration is not a parameter: it is\n% (strike_count - 1) / strike_rate_hz + tail_s.\n\n% Shell geometry (a shallow spherical cap; rise sets the radius of curvature\n% Rs = radius^2 / (2 * rise))\nP.radius_m = 0.230;\nP.rise_m = 0.040;\nP.thickness_m = 0.001;\n\n% Material (bronze, per the paper's measured cymbal)\nP.youngs_modulus = 1.1e11;\nP.poisson_ratio = 0.33;\nP.density = 8860.0;\n\n% Damping: decay rate sigma0 + sigma1 * a^2 * k^2 at flexural wavenumber k\nP.sigma0 = 1.34; % s^-1, frequency-independent\nP.sigma1 = 1.2e-3; % s^-1 per unit nondimensional k^2\n\n% Support at the center: an SDOF washer per unit length of the inner\n% circumference, entering as the shear boundary condition at r = washer\n% radius. support_clamped = 1 replaces it with the rigid clamp of the\n% earlier literature (Bilbao 2010), which is worth hearing once.\nP.washer_radius_m = 0.029;\nP.washer_stiffness = 1.65e5; % N/m per m\nP.washer_mass = 0.0; % kg per m (paper sets 0)\nP.washer_damping = 0.0; % N s/m per m (paper sets 0)\nP.support_clamped = 0;\n\n% Mallet head (SDOF, Hertz-type nonlinear stiffness f = K * delta^alpha)\nP.head_radius_m = 0.025;\nP.head_mass_kg = 0.028;\nP.head_stiffness = 1.6e8; % N m^-alpha\nP.head_alpha = 2.54;\n\n% Mallet stick (Euler-Bernoulli beam, FEM)\nP.stick_length_m = 0.390;\nP.stick_area_m2 = 1.43139e-4;\nP.stick_youngs = 8.8e9;\nP.stick_density = 400.0;\nP.stick_inertia_m4 = 1.63044e-9;\nP.stick_damping_rate = 10.0; % mass-proportional, R = rate * M\nP.n_beam = 20; % beam elements (21 nodes)\n\n% Grip: 0 = rotational spring hand (loose grip), 1 = clamped hand (tight)\nP.grip_clamped = 0;\nP.grip_spring = 140.0; % N m/rad\n\n% Strike gesture: the hand accelerates downward to strike_velocity, coasts\n% until the head touches, holds briefly, and retracts to its rest height.\n% The paper prescribes constant acceleration to the same terminal speed\n% but never reverses, so it cannot play a second stroke. Because the stick\n% is compliant, the tip arrives somewhat slower than the hand.\nP.strike_velocity = 5.75; % peak hand speed, m/s; reproduces the\n % character of the authors' published demo\nP.strike_accel = 250.0; % downstroke hand acceleration, m/s^2\nP.retract_accel = 2500.0; % wrist-snap retract acceleration, m/s^2\nP.strike_hold_ms = 0.5; % dwell between contact and retract\nP.strike_gap_m = -1; % rest clearance; <= 0 sizes it to the stroke\nP.strike_r = 0.6; % contact radius, 0 at the washer, 1 at the rim\nP.strike_theta_deg = 180.0;\nP.strike_count = 1;\nP.strike_rate_hz = 4.0;\nP.velocity_spread = 0.0; % relative velocity scatter of repeated strokes\nP.seed = 1;\n\n% Pickup: normal acceleration of the shell at one point, as in the paper\nP.pickup_r = 0.6;\nP.pickup_theta_deg = 0.0;\n\n% Discretization\nP.n_radial = 24; % radial grid intervals (the paper's N_r)\nP.n_theta = 60; % azimuthal grid points (the paper's N_theta)\nP.g_param = 0.2498; % free parameter of the implicit shell scheme\nP.phi_param = 0.7; % free parameter of the washer and beam schemes\n\n% Rendering / output\nP.sample_rate = 48000;\nP.oversample = 2; % simulation runs at oversample * sample_rate\nP.tail_s = 6.0; % the low modes carry T60 of about 5 s\nP.pre_delay_ms = 3.0;\nP.normalize = 1;\nP.peak_dbfs = -3.0;\nP.full_scale_accel = 4.0e4; % m/s^2 mapped to full scale when normalize = 0\nend\n"
32- },
33- {
34- "name": "cymbal2_demo.m",
35- "content": "% Render a cymbal crash with default parameters and print diagnostics.\n% Set P.strike_velocity = 0.5 for a tap on the same instrument, or\n% P.support_clamped = 1 to hear the rigid center clamp of the earlier\n% literature instead of the washer. P.strike_count = 4 plays a repeated\n% figure; the hand retracts between strokes, which the paper's\n% constant-acceleration hand cannot do.\n% Run with: numbl run cymbal2_demo.m\n\nP = cymbal2_defaults();\ntic;\n[y, dinfo] = cymbal2_render(P);\nelapsed = toc;\n\nfprintf('grid : %d radial x %d azimuthal, %.1f kHz simulation rate\\n', ...\n P.n_radial, P.n_theta, dinfo.simRate / 1000);\nfprintf('shell : a = %.0f cm, rise %.0f mm -> Rs = %.2f m, washer at %.1f cm\\n', ...\n P.radius_m * 100, P.rise_m * 1000, dinfo.RsM, P.washer_radius_m * 100);\nf0 = cymbal2_modes(P, 0); f0 = f0(f0 > 1);\nf1 = cymbal2_modes(P, 1); f1 = f1(f1 > 1);\nf2 = cymbal2_modes(P, 2); f2 = f2(f2 > 1);\nfprintf('low modes : %.1f Hz (m=1 tilt), %.1f Hz (m=2), %.1f Hz (m=0 bounce)\\n', ...\n f1(1), f2(1), f0(1));\nfprintf('stroke : hand %.1f m/s -> tip %.2f m/s at contact\\n', ...\n dinfo.meanV, abs(dinfo.contactVel(1)));\nfprintf('contact : %.0f N peak, %.2f mm compression, %.1f ms over %d event(s)\\n', ...\n dinfo.peakForceN, dinfo.peakXiMm, dinfo.contactMs, dinfo.contactEvents);\nfprintf('shell peak : %.2f mm displacement = %.1f thicknesses\\n', ...\n dinfo.peakDispM * 1000, dinfo.peakDispThickness);\nfprintf('pickup : %.2e m/s^2 peak acceleration at r = %.1f cm\\n', ...\n dinfo.peakAccel, dinfo.pickupRadiusM * 100);\nfprintf('length : %.2f s @ %d Hz\\n', dinfo.lengthS, P.sample_rate);\nfprintf('render time: %.1f s\\n', elapsed);\n"
36- },
37- {
38- "name": "cymbal2_grid.m",
39- "content": "function G = cymbal2_grid(P)\n% CYMBAL2_GRID Polar grid, radial difference matrices, and the azimuthal\n% transform. The radial grid carries l = 1..Nr+2 lines: real nodes at\n% l = 3..Nr spanning [washer radius, rim radius], and two virtual lines\n% beyond each boundary, eliminated by the boundary-condition rows of the\n% assembled systems (Samejima 2021, Fig. 4). Azimuthal operators are\n% circulant, so every linear solve block-diagonalizes over the real\n% cos/sin basis; F maps physical theta samples to coefficients (U * F)\n% and Gi maps back (Uhat * Gi), with mcol giving each column's wavenumber.\n\na = P.radius_m;\nc = P.washer_radius_m;\nNr = P.n_radial;\nnth = P.n_theta;\nif mod(nth, 2) ~= 0\n error('cymbal2: n_theta must be even');\nend\nif Nr < 8\n error('cymbal2: n_radial must be at least 8');\nend\nnrt = Nr + 2;\nhr = (a - c) / (Nr - 3);\nif c <= hr\n error('cymbal2: washer_radius too small for n_radial (need c > (a-c)/(n_radial-3))');\nend\n\nr = c + ((1:nrt)' - 3) * hr; % r(1) = c - 2hr may be tiny on coarse grids;\nr(1) = max(r(1), hr / 2); % row 1 of every operator is never used\n\n% Radial stencils (dense; nrt is small). Rows outside the valid range are\n% left zero and are overwritten by boundary-condition rows downstream.\nD1 = zeros(nrt, nrt);\nD2 = zeros(nrt, nrt);\nD3 = zeros(nrt, nrt);\nfor l = 2:nrt - 1\n D1(l, l - 1) = -1 / (2 * hr); D1(l, l + 1) = 1 / (2 * hr);\n D2(l, l - 1) = 1 / hr^2; D2(l, l) = -2 / hr^2; D2(l, l + 1) = 1 / hr^2;\nend\nfor l = 3:nrt - 2\n D3(l, l - 2) = -1 / (2 * hr^3);\n D3(l, l - 1) = 2 / (2 * hr^3);\n D3(l, l + 1) = -2 / (2 * hr^3);\n D3(l, l + 2) = 1 / (2 * hr^3);\nend\n\n% Azimuthal stencils as dense circulant matrices, applied on the right\n% (U * T2 etc., so they carry a transpose; T2 is symmetric, T1 is\n% antisymmetric and T1t below is arranged so Uth = U * T1t).\nhth = 2 * pi / nth;\nT2 = zeros(nth, nth);\nT1t = zeros(nth, nth);\nfor j = 1:nth\n jm = 1 + mod(j - 2, nth);\n jp = 1 + mod(j, nth);\n T2(j, j) = -2 / hth^2; T2(j, jm) = 1 / hth^2; T2(j, jp) = 1 / hth^2;\n % row j of T1 has -1 at jm, +1 at jp (scaled); transpose: column j\n T1t(jm, j) = -1 / (2 * hth);\n T1t(jp, j) = 1 / (2 * hth);\nend\n\n% Real transform: columns ordered m = 0, (1c,1s), (2c,2s), ..., nth/2\nth = ((1:nth)' - 1) * hth;\nF = zeros(nth, nth);\nGi = zeros(nth, nth);\nmcol = zeros(1, nth);\nF(:, 1) = 1 / nth;\nGi(1, :) = 1;\ncol = 2;\nfor m = 1:nth / 2 - 1\n F(:, col) = 2 * cos(m * th) / nth; Gi(col, :) = cos(m * th)'; mcol(col) = m; col = col + 1;\n F(:, col) = 2 * sin(m * th) / nth; Gi(col, :) = sin(m * th)'; mcol(col) = m; col = col + 1;\nend\nF(:, col) = cos((nth / 2) * th) / nth;\nGi(col, :) = cos((nth / 2) * th)';\nmcol(col) = nth / 2;\n\n% Eigenvalue of the circulant second-difference for wavenumber m\nlam2 = (2 * cos(mcol * hth) - 2) / hth^2;\n\nG.a = a; G.c = c; G.Nr = Nr; G.nrt = nrt; G.nth = nth;\nG.hr = hr; G.hth = hth;\nG.r = r; G.th = th;\nG.D1 = D1; G.D2 = D2; G.D3 = D3;\nG.T1t = T1t; G.T2 = T2;\nG.F = F; G.Gi = Gi; G.mcol = mcol; G.lam2 = lam2;\nG.Rd = 1.0 ./ r; % 1/r column, broadcast over theta\nG.R2d = 1.0 ./ r.^2;\nend\n"
40- },
41- {
42- "name": "cymbal2_mallet.m",
43- "content": "function MB = cymbal2_mallet(P)\n% CYMBAL2_MALLET Precompute the striking mechanism: the mallet stick as an\n% Euler-Bernoulli beam (Hermite cubic elements, consistent mass,\n% mass-proportional damping), held at one end by the hand, with the mallet\n% head attached at the other end as a separate SDOF body (Samejima 2021,\n% Sec. 3). The hand grips either rigidly (grip_clamped = 1: transverse\n% displacement and rotation both prescribed) or through a rotational\n% spring (grip_clamped = 0: displacement prescribed, rotation free against\n% grip_spring). The head couples to the beam tip by the rigid connection\n% u_M = u_S. The paper integrates the head as its own SDOF (Eq. 34) and\n% solves the connection force from the constraint; because that couples an\n% explicit head update to the implicit beam, the discrete pair is\n% unstable, so the head mass is lumped onto the tip node instead, which is\n% the same physics under the rigid connection (substituting u_M = u_S and\n% f_S = -f_M2 into Eq. (31) puts M_M on the tip and f_M1 on the beam).\n\nne = P.n_beam;\nnn = ne + 1;\nndof = 2 * nn;\nle = P.stick_length_m / ne;\nEI = P.stick_youngs * P.stick_inertia_m4;\nrA = P.stick_density * P.stick_area_m2;\n\nKe = EI / le^3 * [12, 6 * le, -12, 6 * le; ...\n 6 * le, 4 * le^2, -6 * le, 2 * le^2; ...\n -12, -6 * le, 12, -6 * le; ...\n 6 * le, 2 * le^2, -6 * le, 4 * le^2];\nMe = rA * le / 420 * [156, 22 * le, 54, -13 * le; ...\n 22 * le, 4 * le^2, 13 * le, -3 * le^2; ...\n 54, 13 * le, 156, -22 * le; ...\n -13 * le, -3 * le^2, -22 * le, 4 * le^2];\n\nK = zeros(ndof, ndof);\nM = zeros(ndof, ndof);\nfor e = 1:ne\n ix = (2 * e - 1):(2 * e + 2);\n K(ix, ix) = K(ix, ix) + Ke;\n M(ix, ix) = M(ix, ix) + Me;\nend\nR = P.stick_damping_rate * M;\nM(ndof - 1, ndof - 1) = M(ndof - 1, ndof - 1) + P.head_mass_kg;\n\n% Node 1 is the grip, node nn the tip. Dof 1 (grip displacement) is always\n% prescribed to the hand motion; dof 2 (grip rotation) is prescribed to 0\n% for the clamped hand, or free with the rotational spring added.\nif P.grip_clamped ~= 0\n fi = 3:ndof;\nelse\n K(2, 2) = K(2, 2) + P.grip_spring;\n fi = 2:ndof;\nend\nKf = K(fi, fi); Mf = M(fi, fi); Rf = R(fi, fi);\nKc = K(fi, 1); Mc = M(fi, 1); Rc = R(fi, 1);\n\nht = 1 / (P.sample_rate * P.oversample);\nps = P.phi_param;\n% Stiffness split with the implicit weight ps/2 on the (k+1, k-1) average,\n% as in the shell's washer row; the complementary split printed in the\n% paper is unstable for the beam's highest elements at audio rates.\nAS = Mf + (ht / 2) * Rf + (ps / 2) * ht^2 * Kf;\nBS = 2 * Mf - (1 - ps) * ht^2 * Kf;\nCS = -Mf + (ht / 2) * Rf - (ps / 2) * ht^2 * Kf;\nASi = inv(AS);\nMB.ASiB = ASi * BS;\nMB.ASiC = ASi * CS;\n\n% Hand coupling, folded so the per-step load is\n% vM * hA(k) + vR * hB(k) + vK * hC(k)\n% with hA = -(uH+ - 2 uH + uH-), hB = -(ht/2)(uH+ - uH-),\n% hC = -ht^2 (ps uH + ((1-ps)/2)(uH+ + uH-)).\nMB.vM = ASi * Mc;\nMB.vR = ASi * Rc;\nMB.vK = ASi * Kc;\n\ntipi = find(fi == ndof - 1);\nMB.tipi = tipi;\nMB.tipcol = ASi(:, tipi);\nMB.gtip = ASi(tipi, tipi);\nMB.nfree = length(fi);\nMB.ht = ht;\nMB.ps = ps;\nend\n"
44- },
45- {
46- "name": "cymbal2_modes.m",
47- "content": "function [f, W] = cymbal2_modes(P, m)\n% CYMBAL2_MODES Eigenfrequencies (Hz) of the linearized shell for one\n% azimuthal wavenumber m, with the membrane (curvature) stiffening and the\n% support boundary condition included; the washer is taken as its static\n% stiffness (its mass and damping are ignored here, as in the paper's\n% linear FEM reference). Used by the demo and the validation, not by the\n% renderer. W returns the mode shapes on the real radial nodes l = 3..Nr.\n\nG = cymbal2_grid(P);\nnrt = G.nrt;\nNr = G.Nr;\nhth = G.hth;\nlam2 = (2 * cos(m * hth) - 2) / hth^2;\nop = cymbal2_perm(P, G, lam2);\n\nE = P.youngs_modulus;\nh = P.thickness_m;\nrhoh = P.density * h;\nRs = P.radius_m^2 / (2 * P.rise_m);\n\n% Airy solve operator: Fm * phi = Z * (E h / Rs) * Lap * u, with\n% phi = phi_r = 0 at both circles and identity rows on the outermost\n% virtual lines.\nFm = zeros(nrt, nrt);\nFm(1, 1) = 1;\nFm(2, :) = G.D1(3, :);\nFm(3, 3) = 1;\nfor l = 4:Nr - 1\n Fm(l, :) = op.Bih(l, :);\nend\nFm(Nr, Nr) = 1;\nFm(Nr + 1, :) = G.D1(Nr, :);\nFm(Nr + 2, Nr + 2) = 1;\nZl = zeros(nrt, nrt);\nfor l = 4:Nr - 1\n Zl(l, :) = op.Lap(l, :);\nend\nPhimap = Fm \\ ((E * h / Rs) * Zl);\n\n% Full stiffness on the equation-of-motion rows l = 3..Nr\nKfull = zeros(nrt, nrt);\nKterm = op.Dflex * op.Bih + (1 / Rs) * op.Lap * Phimap;\nfor l = 3:Nr\n Kfull(l, :) = Kterm(l, :);\nend\n\n% Boundary rows, static. With the washer the shell's equation of motion\n% holds on l = 3..Nr and the four virtual lines are eliminated; with the\n% rigid clamp u and u_r vanish at l = 3, so u(3) itself joins the\n% eliminated set and the equations of motion start at l = 4 (the innermost\n% virtual line is then never engaged).\ne3 = zeros(1, nrt); e3(3) = 1;\nif P.support_clamped ~= 0\n Cmat = [e3; G.D1(3, :); op.MomE; op.ShE];\n keep = 4:Nr;\n elim = [2, 3, Nr + 1, Nr + 2];\nelse\n Cmat = [op.MomW; op.ShW - P.washer_stiffness * e3; op.MomE; op.ShE];\n keep = 3:Nr;\n elim = [1, 2, Nr + 1, Nr + 2];\nend\nX = -Cmat(:, elim) \\ Cmat(:, keep);\nKred = Kfull(keep, keep) + Kfull(keep, elim) * X;\n\n[V, Dm] = eig(Kred / rhoh);\nom2 = real(diag(Dm));\n[om2, idx] = sort(om2);\nom2 = max(om2, 0);\nf = sqrt(om2) / (2 * pi);\nW = real(V(:, idx));\nend\n"
48- },
49- {
50- "name": "cymbal2_output.m",
51- "content": "function [y, peakAccel] = cymbal2_output(P, acc)\n% CYMBAL2_OUTPUT Output conditioning: DC blocker on the raw pickup\n% acceleration, Kaiser-windowed-sinc anti-alias decimation from the\n% simulation rate to the output rate, tail fade, and normalization. Same\n% structure as the sibling instruments' output stages.\n\npeakAccel = max(abs(acc));\nR = P.oversample;\nT = length(acc);\n\n% One-pole DC blocker (the pickup is an acceleration, but the contact\n% leaves a small quasi-static offset while the washer relaxes)\nfc = 15.0;\na1 = exp(-2 * pi * fc / (P.sample_rate * R));\np = zeros(1, T);\nxp = 0; yp = 0;\nfor k = 1:T\n yv = acc(k) - xp + a1 * yp;\n xp = acc(k);\n yp = yv;\n p(k) = yv;\nend\n\nif R > 1\n taps = 12 * R + 1;\n M2 = (taps - 1) / 2;\n fcNorm = 0.45 / R;\n m = (0:taps - 1) - M2;\n hh = sin(2 * pi * fcNorm * m) ./ (pi * m);\n hh(M2 + 1) = 2 * fcNorm;\n r = m / M2;\n beta = 8;\n hh = hh .* (c2_i0(beta * sqrt(max(0, 1 - r.^2))) / c2_i0(beta));\n hh = hh / sum(hh);\n pf = conv(p, hh);\n outLen = floor(T / R);\n y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen);\nelse\n y = p;\nend\n\n% Squared fade on the tail\nfade = min(2048, floor(length(y) / 8));\nif fade > 1\n y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;\nend\n\nif P.normalize ~= 0\n pk = max(abs(y));\n if pk > 0\n y = y / pk * 10^(P.peak_dbfs / 20);\n end\nelse\n y = y / P.full_scale_accel;\nend\nend\n\nfunction v = c2_i0(x)\n% Modified Bessel function of the first kind, order zero (series)\nv = ones(size(x));\nterm = ones(size(x));\nhx = x / 2;\nfor k = 1:40\n term = term .* (hx / k).^2;\n v = v + term;\nend\nend\n"
52- },
53- {
54- "name": "cymbal2_perm.m",
55- "content": "function op = cymbal2_perm(P, G, lam2)\n% CYMBAL2_PERM Radial operators for one azimuthal wavenumber, dimensional.\n% lam2 is the (negative) eigenvalue of the circulant second-difference in\n% theta for this wavenumber, so every d^2/dtheta^2 becomes multiplication\n% by lam2. Returns the Laplacian, the biharmonic (as the exact discrete\n% composition Lap * Lap), and the radial-moment and Kelvin-Kirchhoff shear\n% rows used as boundary conditions at the washer (l = 3) and at the free\n% edge (l = Nr). The moment and shear expressions are those of Samejima\n% (2021), Eqs. (6)-(7) and (14)-(15).\n\nnrt = G.nrt;\nnu = P.poisson_ratio;\nDfl = P.youngs_modulus * P.thickness_m^3 / (12 * (1 - nu^2));\n\nRm = diag(G.Rd);\nR2 = diag(G.R2d);\nR3 = diag(G.Rd.^3);\n\nLap = G.D2 + Rm * G.D1 + lam2 * R2;\nBih = Lap * Lap;\n\n% Radial bending moment (per unit length, divided by -D):\n% u_rr + nu * (u_r / r + u_thth / r^2)\nMom = G.D2 + nu * (Rm * G.D1 + lam2 * R2);\n\n% Kelvin-Kirchhoff shear: V = -D * [u_rrr + u_rr / r - u_r / r^2\n% + (nu - 3) / r^3 * u_thth + (2 - nu) / r^2 * u_r_thth]\nSh = -Dfl * (G.D3 + Rm * G.D2 - R2 * G.D1 ...\n + (nu - 3) * lam2 * R3 + (2 - nu) * lam2 * R2 * G.D1);\n\nop.Lap = Lap;\nop.Bih = Bih;\nop.MomW = Mom(3, :); % moment row at the washer circle\nop.MomE = Mom(nrt - 2, :); % moment row at the free edge (l = Nr)\nop.ShW = Sh(3, :);\nop.ShE = Sh(nrt - 2, :);\nop.Dflex = Dfl;\nend\n"
56- },
57- {
58- "name": "cymbal2_render.m",
59- "content": "function [y, dinfo] = cymbal2_render(P)\n% CYMBAL2_RENDER Full render: precompute -> time integration -> output.\n\nS = cymbal2_system(P);\nMB = cymbal2_mallet(P);\n[acc, dinfo] = cymbal2_sim(P, S, MB);\n[y, peakAccel] = cymbal2_output(P, acc);\n\ndinfo.peakAccel = peakAccel;\ndinfo.strikeRadiusM = S.strikeRadiusM;\ndinfo.pickupRadiusM = S.pickupRadiusM;\ndinfo.RsM = S.Rs;\ndinfo.lengthS = length(y) / P.sample_rate;\nend\n"
60- },
61- {
62- "name": "cymbal2_sim.m",
63- "content": "function [acc, dinfo, fhist] = cymbal2_sim(P, S, MB)\n% CYMBAL2_SIM Time-domain integration of the full system: von Karman\n% shallow spherical shell (implicit in the stiff linear terms, explicit in\n% the nonlinearity), washer support built into the boundary rows, mallet\n% head + stick coupled through the Hertz contact. Returns the normal\n% acceleration of the shell at the pickup node at the simulation rate,\n% diagnostics, and the contact force history.\n\nG = S.G;\nnrt = G.nrt;\nNr = G.Nr;\nnth = G.nth;\nht = S.ht;\nrhoh = S.rhoh;\nRs = S.Rs;\nEhRs = S.Eh / Rs;\nEh2 = S.Eh / 2;\nls = S.ls; js = S.js; lp = S.lp; jp = S.jp;\nfsc = S.forceScale;\nnlr0 = S.eomRow0;\n\nK = max(1, round(P.strike_count));\n\n% Stroke schedule; per-stroke peak hand speeds with a seeded spread. Each\n% stroke is launched early by its predicted hand travel time (acceleration\n% plus coast over the rest clearance), so the contacts land near the\n% nominal beats instead of drifting with velocity. A nonpositive\n% strike_gap_m asks for the clearance to be sized to the stroke: contact\n% arrives just after the hand reaches full speed, which keeps the tip's\n% whip phase, and with it the contact speed, consistent across strokes.\ntstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);\nvstr = zeros(1, K);\nstate = P.seed * 2654435761;\nfor i = 1:K\n state = mod(state * 1664525 + 1013904223, 4294967296);\n x = state / 4294967296;\n vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));\nend\ngap = P.strike_gap_m;\nif gap <= 0\n gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;\nend\ntravel = gap ./ vstr + vstr / (2 * P.strike_accel);\ntlaunch = tstr - travel;\noff = -min(0, min(tlaunch));\ntlaunch = tlaunch + off;\n\nduration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;\nnsteps = floor(duration * P.sample_rate * P.oversample);\n\nKM = P.head_stiffness;\nalM = P.head_alpha;\ntipi = MB.tipi;\nps = MB.ps;\n\n% Hand state machine: the hand accelerates downward, coasts at the stroke\n% speed until the head makes contact (or a miss depth), holds briefly,\n% then retracts at half speed and brakes so it arrives back at rest,\n% ready for the next scheduled stroke. This replaces the paper's\n% constant-acceleration hand, which never reverses. The machine runs one\n% step ahead because the beam coupling needs u_H at k + 1.\naH = P.strike_accel;\naHr = P.retract_accel;\nholdSteps = max(1, round(P.strike_hold_ms * 1e-3 / ht));\nmissDepth = gap + 0.06;\nuHm = 0; uH0 = 0; uHp = 0;\nvH = 0;\nphase = 0;\nstrokeIdx = 1;\nholdLeft = 0;\nvtar = 0;\n\nD1p = G.D1; D2p = G.D2;\nT1t = G.T1t; T2 = G.T2;\nF = G.F; Gi = G.Gi;\nRd = G.Rd; R2d = G.R2d;\nAiB = S.AiB; AiC = S.AiC; Ai = S.Ai; Fi = S.Fi;\nsNL = ht^2 / rhoh;\n\nU = zeros(nrt, nth);\nU1 = zeros(nrt, nth);\nUh3 = zeros(nrt, 1, nth);\nUh13 = zeros(nrt, 1, nth);\nd = zeros(MB.nfree, 1);\nd1 = zeros(MB.nfree, 1);\n\nacc = zeros(1, nsteps);\nfhist = zeros(1, nsteps);\nshellActive = 0;\npeakDisp = 0;\npeakXi = 0;\ninContact = 0;\ncvel = zeros(1, K);\n\nfor k = 1:nsteps\n % ---- mallet: contact force from the state at step k ----\n xi = U(ls, js) - d(tipi) - gap;\n if xi > 0\n fM1 = KM * xi^alM;\n shellActive = 1;\n if inContact == 0 && strokeIdx <= K && cvel(strokeIdx) == 0\n cvel(strokeIdx) = (d(tipi) - d1(tipi)) / ht; % tip speed at contact\n end\n inContact = 1;\n if xi > peakXi\n peakXi = xi;\n end\n else\n fM1 = 0;\n inContact = 0;\n end\n fhist(k) = fM1;\n hA = -(uHp - 2 * uH0 + uHm);\n hB = -(ht / 2) * (uHp - uHm);\n hC = -ht^2 * (ps * uH0 + ((1 - ps) / 2) * (uHp + uHm));\n dnew = MB.ASiB * d + MB.ASiC * d1 + MB.vM * hA + MB.vR * hB + MB.vK * hC ...\n + (ht^2 * fM1) * MB.tipcol;\n\n % ---- advance the hand one step ----\n if phase == 0\n if strokeIdx <= K && k * ht >= tlaunch(strokeIdx)\n phase = 1;\n vtar = vstr(strokeIdx);\n end\n elseif phase == 1\n vH = max(vH - aH * ht, -vtar);\n if inContact ~= 0\n phase = 3;\n holdLeft = holdSteps;\n elseif vH <= -vtar\n phase = 2;\n end\n elseif phase == 2\n if inContact ~= 0\n phase = 3;\n holdLeft = holdSteps;\n elseif uHp <= -missDepth\n phase = 4;\n end\n elseif phase == 3\n holdLeft = holdLeft - 1;\n if holdLeft <= 0\n phase = 4;\n end\n elseif phase == 4\n vH = min(vH + aHr * ht, 0.5 * vtar);\n if vH > 0 && uHp >= -vH^2 / (2 * aHr)\n phase = 5;\n end\n else\n vH = vH - aHr * ht;\n if vH <= 0 || uHp >= 0\n vH = 0;\n phase = 0;\n strokeIdx = strokeIdx + 1;\n end\n end\n uHm = uH0;\n uH0 = uHp;\n uHp = uHp + vH * ht;\n\n if shellActive ~= 0\n % ---- shell derivatives and von Karman brackets at step k ----\n Ur = D1p * U;\n Urr = D2p * U;\n v1 = Rd .* Ur + R2d .* (U * T2);\n v2 = Rd .* (Ur * T1t) - R2d .* (U * T1t);\n LapU = Urr + v1;\n rhsA = EhRs * LapU - Eh2 * 2 * (Urr .* v1 - v2 .* v2);\n rhsA(1:3, :) = 0;\n rhsA(Nr:nrt, :) = 0;\n Ph3 = pagemtimes(Fi, reshape(rhsA * F, nrt, 1, nth));\n Phi = reshape(Ph3, nrt, nth) * Gi;\n Pr = D1p * Phi;\n Prr = D2p * Phi;\n w1 = Rd .* Pr + R2d .* (Phi * T2);\n w2 = Rd .* (Pr * T1t) - R2d .* (Phi * T1t);\n NL = Prr .* v1 + Urr .* w1 - 2 * w2 .* v2 - (Prr + w1) / Rs;\n NL(ls, js) = NL(ls, js) - fM1 * fsc;\n NL(1:nlr0 - 1, :) = 0;\n NL(Nr + 1:nrt, :) = 0;\n NLh = pagemtimes(Ai, reshape((sNL * NL) * F, nrt, 1, nth));\n Uh3n = pagemtimes(AiB, Uh3) + pagemtimes(AiC, Uh13) + NLh;\n Unew = reshape(Uh3n, nrt, nth) * Gi;\n\n acc(k) = (Unew(lp, jp) - 2 * U(lp, jp) + U1(lp, jp)) / ht^2;\n U1 = U;\n U = Unew;\n Uh13 = Uh3;\n Uh3 = Uh3n;\n if mod(k, 512) == 0\n pd = max(abs(U(:, js)));\n if pd > peakDisp\n peakDisp = pd;\n end\n if ~isfinite(pd)\n error('cymbal2: simulation diverged at t = %.4f s (reduce strike_velocity, or raise n_theta / oversample)', k * ht);\n end\n end\n end\n\n d1 = d;\n d = dnew;\nend\n\n% Contact bookkeeping from the force history\ninC = fhist > 0;\ndinfo.peakForceN = max(fhist);\ndinfo.contactSteps = sum(inC);\ndinfo.contactMs = 1000 * ht * dinfo.contactSteps;\ndinfo.contactEvents = sum(diff([0, inC]) > 0);\ndinfo.peakXiMm = peakXi * 1000;\ndinfo.peakDispM = peakDisp;\ndinfo.peakDispThickness = peakDisp / P.thickness_m;\ndinfo.peakAccel = max(abs(acc));\ndinfo.nsteps = nsteps;\ndinfo.simRate = P.sample_rate * P.oversample;\ndinfo.tstr = tstr;\ndinfo.vstr = vstr;\ndinfo.meanV = mean(vstr);\ndinfo.contactVel = cvel; % tip speed at first contact of each stroke\ndinfo.strokesPlayed = min(strokeIdx - 1 + (phase > 0), K);\nend\n"
64- },
65- {
66- "name": "cymbal2_system.m",
67- "content": "function S = cymbal2_system(P)\n% CYMBAL2_SYSTEM Precompute everything the time loop needs. The linear\n% updates block-diagonalize over azimuthal wavenumber, so the implicit\n% shell solve and the Airy (stress function) solve are stored as one small\n% radial matrix per transform column: the time step becomes\n% Uhat+ = (A^-1 B) Uhat + (A^-1 C) Uhat- + A^-1 (ht^2/rho h) NLhat\n% with A, B, C the matrices of Samejima (2021) Eq. (20) restricted to one\n% wavenumber, the boundary-condition rows of Eqs. (22)-(30) built in, and\n% NL the explicit nonlinear + forcing terms evaluated in physical space.\n% Rows are equilibrated before inversion; the stored inverses absorb the\n% scaling so the loop applies them to raw right-hand sides.\n\nG = cymbal2_grid(P);\nnrt = G.nrt;\nNr = G.Nr;\nnth = G.nth;\n\nh = P.thickness_m;\nrhoh = P.density * h;\nEh = P.youngs_modulus * h;\nRs = P.radius_m^2 / (2 * P.rise_m);\nht = 1 / (P.sample_rate * P.oversample);\ng = P.g_param;\nphiw = P.phi_param;\nkap = 0; % set below from cymbal2_perm (needs Dflex)\ns1 = P.sigma1 * P.radius_m^2;\ns0 = P.sigma0;\n\nMW = P.washer_mass;\nRW = P.washer_damping;\nKW = P.washer_stiffness;\n% Spring term split as phiw/2 on the (k+1, k-1) average and (1 - phiw) on\n% step k. Note the paper prints the complementary split, which is unstable\n% for a massless washer (the roots of (1-phi)/2 (z^2+1) + phi z include\n% |z| > 1 whenever phi > 1/2); with the implicit weight on the average,\n% phi = 0.7 is unconditionally stable.\nAW = MW + (ht / 2) * RW + (phiw / 2) * ht^2 * KW;\nBW = 2 * MW - (1 - phiw) * ht^2 * KW;\nCW = -MW + (ht / 2) * RW - (phiw / 2) * ht^2 * KW;\n\nAiB = zeros(nrt, nrt, nth);\nAiC = zeros(nrt, nrt, nth);\nAi = zeros(nrt, nrt, nth);\nFi = zeros(nrt, nrt, nth);\nnDropped = 0;\n\ne1 = zeros(1, nrt); e1(1) = 1;\ne3 = zeros(1, nrt); e3(3) = 1;\nI = eye(nrt);\n\nfor m = 0:nth / 2\n lam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;\n op = cymbal2_perm(P, G, lam2);\n kap = op.Dflex / rhoh;\n\n Am = zeros(nrt, nrt);\n Bm = zeros(nrt, nrt);\n Cm = zeros(nrt, nrt);\n if P.support_clamped ~= 0\n l0 = 4;\n else\n l0 = 3;\n end\n for l = l0:Nr\n Am(l, :) = (1 + s0 * ht) * I(l, :) + g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);\n Bm(l, :) = 2 * I(l, :) - (1 - 2 * g) * kap * ht^2 * op.Bih(l, :);\n Cm(l, :) = (s0 * ht - 1) * I(l, :) - g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);\n end\n if P.support_clamped ~= 0\n Am(1, :) = e1;\n Am(2, :) = e3;\n Am(3, :) = G.D1(3, :);\n else\n % The shear side of the washer balance gets the same implicit split as\n % the spring; the paper evaluates it at step k only, which feeds a\n % spurious mode on the virtual lines when the washer is massless.\n Am(1, :) = op.MomW;\n Am(2, :) = AW * e3 - (phiw / 2) * ht^2 * op.ShW;\n Bm(2, :) = BW * e3 + (1 - phiw) * ht^2 * op.ShW;\n Cm(2, :) = CW * e3 + (phiw / 2) * ht^2 * op.ShW;\n end\n Am(Nr + 1, :) = op.MomE;\n Am(Nr + 2, :) = op.ShE;\n\n % Row equilibration, absorbed into the stored inverses\n sc = 1.0 ./ max(abs(Am), [], 2);\n Asc = sc .* Am;\n Aisc = inv(Asc);\n AiBm = Aisc * (sc .* Bm);\n AiCm = Aisc * (sc .* Cm);\n Aim = Aisc .* sc'; % inv(Am) applied to raw rhs\n\n % On grids too coarse to resolve a wavenumber's inner boundary layer,\n % the discretized boundary conditions can hand that wavenumber a\n % negative stiffness; such columns are truncated rather than integrated.\n % At the default grid none are.\n sr = max(abs(eig([AiBm, AiCm; I, zeros(nrt)])));\n if sr > 1 + 1e-7\n AiBm = zeros(nrt, nrt);\n AiCm = zeros(nrt, nrt);\n Aim = zeros(nrt, nrt);\n nDropped = nDropped + 1 + (m > 0 && m < nth / 2);\n end\n\n % Airy solve: biharmonic rows inside, phi = phi_r = 0 at both circles\n Fm = zeros(nrt, nrt);\n Fm(1, 1) = 1;\n Fm(2, :) = G.D1(3, :);\n Fm(3, 3) = 1;\n for l = 4:Nr - 1\n Fm(l, :) = op.Bih(l, :);\n end\n Fm(Nr, Nr) = 1;\n Fm(Nr + 1, :) = G.D1(Nr, :);\n Fm(Nr + 2, Nr + 2) = 1;\n scf = 1.0 ./ max(abs(Fm), [], 2);\n Fim = inv(scf .* Fm) .* scf';\n\n if m == 0\n cols = 1;\n elseif m == nth / 2\n cols = nth;\n else\n cols = [2 * m, 2 * m + 1];\n end\n for cc = 1:length(cols)\n AiB(:, :, cols(cc)) = AiBm;\n AiC(:, :, cols(cc)) = AiCm;\n Ai(:, :, cols(cc)) = Aim;\n Fi(:, :, cols(cc)) = Fim;\n end\nend\n\nif nDropped > nth / 2\n error(['cymbal2: the radial grid cannot resolve the inner boundary layer ' ...\n '(need roughly washer_radius > 3 (a - c) / (n_radial - 3)); raise n_radial']);\nend\n\n% Strike and pickup nodes (radial fraction 0 at the washer, 1 at the rim)\nls = 3 + round(min(1, max(0, P.strike_r)) * (Nr - 3));\nlp = 3 + round(min(1, max(0, P.pickup_r)) * (Nr - 3));\njs = 1 + mod(round(P.strike_theta_deg / 360 * nth), nth);\njp = 1 + mod(round(P.pickup_theta_deg / 360 * nth), nth);\n\nS.G = G;\nS.AiB = AiB; S.AiC = AiC; S.Ai = Ai; S.Fi = Fi;\nS.rhoh = rhoh; S.Eh = Eh; S.Rs = Rs; S.ht = ht;\nS.ls = ls; S.lp = lp; S.js = js; S.jp = jp;\nS.eomRow0 = 3 + (P.support_clamped ~= 0); % first equation-of-motion row\nS.nDropped = nDropped; % truncated theta columns\nS.strikeRadiusM = G.r(ls);\nS.pickupRadiusM = G.r(lp);\nS.forceScale = 1 / (G.r(ls) * G.hr * G.hth); % point force -> pressure\nend\n"
68- },
69- {
70- "name": "cymbal2_validate.m",
71- "content": "% Validate the cymbal2 model. There is no external reference\n% implementation, so the checks are (a) closed-form limits the\n% discretization must reproduce (the free circular plate, rigid-body\n% motion on the washer, cantilever statics of the stick), (b) agreement\n% with the linear-FEM eigenfrequencies published in Samejima (2021) for\n% the modes that model resolves the same way, (c) structural identities\n% of the transform and of the assembled update, (d) exact reproducibility,\n% and (e) pinned end-to-end diagnostics from a reduced configuration,\n% recorded when the model was built and treated as a regression reference.\n% Run with: numbl run cymbal2_validate.m\n\nfails = 0;\nP = cymbal2_defaults();\n\n% ---- 1. flat limit: the completely free circular plate ----\n% With the curvature removed and the washer stiffness zeroed, the shell\n% must reproduce the free plate's frequency equation (the 2 x 2 J_n / I_n\n% determinant under zero radial moment and zero Kelvin-Kirchhoff shear).\n% m = 0 is excluded: its limit as the free hole closes is logarithmic, so\n% a finite hole shifts it far more than the m >= 1 modes.\nQ = P;\nQ.rise_m = 1e-9;\nQ.washer_stiffness = 0;\nQ.washer_radius_m = 0.010;\nQ.n_radial = 96;\nDp = Q.youngs_modulus * Q.thickness_m^3 / (12 * (1 - Q.poisson_ratio^2));\nscal = sqrt(Dp / (Q.density * Q.thickness_m)) / Q.radius_m^2 / (2 * pi);\nerrP = 0;\nfor m = 1:4\n f = cymbal2_modes(Q, m);\n fe = f(f > 0.5);\n lam = c2v_plate_roots(m, Q.poisson_ratio);\n fc = lam.^2 * scal;\n nc = min(3, min(length(fe), length(fc)));\n errP = max(errP, max(abs(fe(1:nc) - fc(1:nc)') ./ fc(1:nc)'));\nend\nok = errP < 0.04;\nfails = fails + (~ok);\nfprintf('flat plate : max freq err %.1e over m = 1..4 %s\\n', errP, c2v_verdict(ok));\n\n% ---- 2. rigid limit: bounce and tilt on the washer ----\n% With the plate stiffened a hundredfold the lowest m = 0 and m = 1 modes\n% must approach the closed-form rigid-body values on the line spring:\n% bounce f^2 = KW 2 pi c / (rho h pi (a^2 - c^2)) / (2 pi)^2 and tilt\n% f^2 = 4 KW c^3 / (rho h (a^4 - c^4)) / (2 pi)^2.\nQ = P;\nQ.rise_m = 1e-9;\nQ.youngs_modulus = P.youngs_modulus * 100;\nQ.n_radial = 96;\na = Q.radius_m; c = Q.washer_radius_m;\nrhoh = Q.density * Q.thickness_m;\nKW = Q.washer_stiffness;\nf0r = sqrt(KW * 2 * pi * c / (rhoh * pi * (a^2 - c^2))) / (2 * pi);\nf1r = sqrt(4 * KW * c^3 / (rhoh * (a^4 - c^4))) / (2 * pi);\nfa = cymbal2_modes(Q, 0); fa = fa(fa > 0.01);\nfb = cymbal2_modes(Q, 1); fb = fb(fb > 0.01);\ne0 = abs(fa(1) / f0r - 1);\ne1 = abs(fb(1) / f1r - 1);\nok = e0 < 0.03 && e1 < 0.03;\nfails = fails + (~ok);\nfprintf('washer : bounce err %.1e, tilt err %.1e vs rigid closed form %s\\n', e0, e1, c2v_verdict(ok));\n\n% ---- 3. the paper's linear FEM reference, where comparable ----\n% Samejima's Fig. 9 gives linear-FEM eigenfrequencies of the same shell.\n% For the clamped center, the (2,1), (3,1), (4,1) modes must agree to a\n% few percent ((1,1) is excluded: it is controlled by the support model,\n% where this shell's moment-free hole and their 2-D elements differ).\nQ = P;\nQ.support_clamped = 1;\nref = [20.1, 43.3, 77.7];\ngot = zeros(1, 3);\nfor m = 2:4\n f = cymbal2_modes(Q, m);\n f = f(f > 1);\n got(m - 1) = f(1);\nend\nerrF = max(abs(got - ref) ./ ref);\nok = errF < 0.06;\nfails = fails + (~ok);\nfprintf('paper FEM : clamped (2,1)(3,1)(4,1) = %.1f %.1f %.1f Hz (ref %.1f %.1f %.1f), err %.1e %s\\n', ...\n got(1), got(2), got(3), ref(1), ref(2), ref(3), errF, c2v_verdict(ok));\n\n% ---- 4. transform and update-structure identities ----\n% The azimuthal transform must be exactly invertible, and the physical\n% derivative pipeline must agree with the per-wavenumber operators: for a\n% field cos(m theta) g(r), the physical-space Laplacian must equal the\n% per-m radial Laplacian applied to g.\nG = cymbal2_grid(P);\neT = max(max(abs(G.F * G.Gi - eye(G.nth))));\ngr = exp(-((G.r - 0.12) / 0.05).^2);\nm = 7;\nUf = gr * cos(m * G.th');\nLapPhys = G.D2 * Uf + G.Rd .* (G.D1 * Uf) + G.R2d .* (Uf * G.T2);\nlam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;\nop = cymbal2_perm(P, G, lam2);\nLapM = (op.Lap * gr) * cos(m * G.th');\nr2 = 2:G.nrt - 1;\neL = max(max(abs(LapPhys(r2, :) - LapM(r2, :)))) / max(max(abs(LapM(r2, :))));\nok = eT < 1e-10 && eL < 1e-10;\nfails = fails + (~ok);\nfprintf('operators : transform inverse %.1e, physical vs per-m Laplacian %.1e %s\\n', ...\n eT, eL, c2v_verdict(ok));\n\n% ---- 5. the assembled update is not amplifying ----\n% Every azimuthal column's two-step update must have spectral radius at\n% most 1 (the damping makes physical modes strictly contractive; boundary\n% constraints contribute zeros). The system builder truncates columns\n% that fail; at the defaults none may be.\nS = cymbal2_system(P);\nnrt = S.G.nrt;\nsrMax = 0;\nfor col = [1, 2, S.G.nth / 2, S.G.nth]\n T = [S.AiB(:, :, col), S.AiC(:, :, col); eye(nrt), zeros(nrt)];\n srMax = max(srMax, max(abs(eig(T))));\nend\nok = S.nDropped == 0 && srMax <= 1 + 1e-9;\nfails = fails + (~ok);\nfprintf('stability : %d truncated columns, sampled spectral radius %.9f %s\\n', ...\n S.nDropped, srMax, c2v_verdict(ok));\n\n% ---- 6. stick statics ----\n% A static tip load on the clamped-hand stick must deflect by\n% F L^3 / (3 E I); the rotational-spring hand adds F L^2 / K_R times L.\nMBc = cymbal2_mallet(setfield(P, 'grip_clamped', 1));\nMBs = cymbal2_mallet(P);\nEI = P.stick_youngs * P.stick_inertia_m4;\nL = P.stick_length_m;\nFl = 1.0;\n% static balance of the update: at a fixed point, (I - A^-1 B - A^-1 C) d\n% equals A^-1 ht^2 K d, so the deflection under a constant tip load F is\n% (I - ASiB - ASiC) \\ (tipcol * ht^2 * F)\nhtv = MBc.ht;\ndc = ((eye(size(MBc.ASiB, 1)) - MBc.ASiB - MBc.ASiC) \\ (MBc.tipcol * htv^2)) * Fl;\nds = ((eye(size(MBs.ASiB, 1)) - MBs.ASiB - MBs.ASiC) \\ (MBs.tipcol * htv^2)) * Fl;\nwc = dc(MBc.tipi);\nws = ds(MBs.tipi);\nwcRef = Fl * L^3 / (3 * EI);\nwsRef = wcRef + (Fl * L / P.grip_spring) * L; % grip rotation F L / K_R times the arm\nec = abs(wc / wcRef - 1);\nes = abs(ws / wsRef - 1);\nok = ec < 1e-6 && es < 1e-6;\nfails = fails + (~ok);\nfprintf('stick : clamped-hand tip %.3e (ref %.3e), spring-hand %.3e (ref %.3e) %s\\n', ...\n wc, wcRef, ws, wsRef, c2v_verdict(ok));\n\n% ---- 7. reproducibility and pinned diagnostics ----\n% Two identical renders of a reduced configuration must agree exactly,\n% and its diagnostics are pinned. The trajectory of a nonlinear shell is\n% sensitive, so the acoustic pins carry percent-level tolerances while\n% the contact pins are tight.\nQ = P;\nQ.sample_rate = 44100;\nQ.oversample = 1;\nQ.tail_s = 1.5;\n[y1, d1] = cymbal2_render(Q);\n[y2, d2] = cymbal2_render(Q);\nrms1 = sqrt(mean(y1.^2));\nrms2 = sqrt(mean(y2.^2));\nsame = max(abs(y1 - y2)) == 0;\nok = same ...\n && d1.contactEvents == d2.contactEvents ...\n && abs(d1.peakForceN - 311.14) < 4 ...\n && abs(d1.contactVel(1) + 8.223) < 0.08 ...\n && abs(rms1 - 0.0973) < 0.010;\nfails = fails + (~ok);\nfprintf('pinned : repro %d, peakF %.2f N (ref 311.14), cvel %.3f (ref -8.223), rms %.4f (ref 0.0973) %s\\n', ...\n same, d1.peakForceN, d1.contactVel(1), rms1, c2v_verdict(ok));\n\n% ---- 8. the crash decays ----\nn = length(y1);\nfs = Q.sample_rate;\nhead = y1(1:round(0.2 * fs));\ntail = y1(n - round(0.2 * fs):n);\ndropE = 10 * log10(sum(tail.^2) / sum(head.^2));\nok = dropE < -3 && dropE > -50;\nfails = fails + (~ok);\nfprintf('decay : last/first 0.2 s %.1f dB %s\\n', dropE, c2v_verdict(ok));\n\nif fails > 0\n fprintf('%d CHECK(S) FAILED\\n', fails);\nelse\n fprintf('all checks passed\\n');\nend\n\nfunction lam = c2v_plate_roots(n, nu)\n% Roots of the free-edge frequency equation of a circular plate\nlam = [];\nprev = c2v_plate_det(n, nu, 0.30);\nfor x = 0.30:0.01:20\n cur = c2v_plate_det(n, nu, x + 0.01);\n if prev * cur < 0\n lo = x; hi = x + 0.01;\n for it = 1:60\n mid = (lo + hi) / 2;\n if c2v_plate_det(n, nu, lo) * c2v_plate_det(n, nu, mid) < 0\n hi = mid;\n else\n lo = mid;\n end\n end\n lam = [lam, (lo + hi) / 2];\n end\n prev = cur;\nend\nend\n\nfunction d = c2v_plate_det(n, nu, x)\nJ = besselj(n, x); I = besseli(n, x);\nJp = besselj(n - 1, x) - (n / x) * J;\nIp = besseli(n - 1, x) - (n / x) * I;\nJpp = -Jp / x + (n^2 / x^2 - 1) * J;\nIpp = -Ip / x + (n^2 / x^2 + 1) * I;\na11 = Jpp + nu * (Jp / x - n^2 * J / x^2);\na12 = Ipp + nu * (Ip / x - n^2 * I / x^2);\na21 = -Jp - (1 - nu) * (n^2 / x^2) * (Jp - J / x);\na22 = Ip - (1 - nu) * (n^2 / x^2) * (Ip - I / x);\nd = a11 * a22 - a12 * a21;\nend\n\nfunction s = c2v_verdict(ok)\nif ok\n s = 'ok';\nelse\n s = 'MISMATCH';\nend\nend\n"
72- }
73- ]
74-}
instrument/.gitignoreadded+1−0View file
@@ -0,0 +1 @@
1+node_modules/
instrument/core.tsadded+851−0View file
@@ -0,0 +1,851 @@
1+// The cymbal (Samejima 2021 nonlinear FD shell) as a track-render
2+// instrument: pure rendering logic with no DOM or Worker dependency,
3+// shared by the browser worker (worker.ts) and the node tests (test.ts).
4+// Two request shapes arrive here:
5+//
6+// - the track path (midiboot): every note of a MIDI track. Notes are
7+// grouped by cache key (pitch and bucketed velocity; a crash ignores
8+// note duration, the ring is the physics), each distinct stroke is a
9+// full time integration, and copies are summed at the onsets. MIDI
10+// pitch scales the shell geometry around pitch 60, exactly as the
11+// retired per-note export driver did.
12+//
13+// - the explorer path (the cymbal2 web app, wantDiagnostics set): the
14+// full parameter surface and live-edited .m sources drive one render
15+// with the diagnostics payload the app's plots consume.
16+//
17+// Either way the MATLAB implementation in ../matlab remains the source of
18+// truth; what varies is who runs the time loop. When WebGPU is available
19+// the unmodified setup stage runs through numbl, the operators go to the
20+// WGSL port of cymbal2_sim.m (./gpuSim.ts), and the acceleration comes
21+// back through the unmodified cymbal2_output.m; this is 20x or more
22+// faster than interpreting the loop, and the GPU chunk boundaries give
23+// integration progress and cooperative cancellation. Without WebGPU, or
24+// when cymbal2_sim.m has been edited in the app (the WGSL mirrors the
25+// shipped step), the original interpreter path runs instead.
26+
27+import { executeCode, RTV, type RuntimeValue, type RuntimeTensor } from 'numbl';
28+import type { MFile, TrackRenderRequest } from './types.ts';
29+import { CYMBAL2_KNOBS, type Cymbal2Variant } from './knobs.ts';
30+import { requestGpuDevice, runGpuSim, type GpuSetup } from './gpuSim.ts';
31+import { NoteCache } from './noteCache.ts';
32+
33+export { CYMBAL2_KNOBS, type Cymbal2Variant };
34+
35+/** Thrown by env.checkCancel to abandon a render; the worker maps it to a
36+ * cancelled response rather than an error. */
37+export class RenderCancelled extends Error {
38+ constructor() {
39+ super('render cancelled');
40+ this.name = 'RenderCancelled';
41+ }
42+}
43+
44+/** The environment the host provides: sources, progress, cancellation, and
45+ * a macrotask yield so a browser worker can receive a cancel message. */
46+export interface RenderEnv {
47+ mfiles: MFile[];
48+ /** The shipped sources, for the cymbal2_sim.m edited-step check. */
49+ shipped: MFile[];
50+ post: (frac: number, stage: string) => void;
51+ /** Throws RenderCancelled to abandon the render; checked between stages,
52+ * between notes, and between GPU chunks. */
53+ checkCancel: () => void;
54+ yieldToHost: () => Promise<void>;
55+}
56+
57+// --- diagnostics payload for the explorer path (the web app's plots) ---
58+
59+export interface Cymbal2ExploreDiag {
60+ quality: 'preview' | 'full';
61+ forceT: Float32Array; // contact force history (max-pooled), s
62+ forceN: Float32Array; // and its values, N
63+ modeF: Float32Array; // linear mode frequencies, Hz
64+ modeM: Float32Array; // nodal diameters of each
65+ diag: {
66+ strokesPlayed: number;
67+ meanV: number;
68+ contactVel: number;
69+ peakForceN: number;
70+ contactMs: number;
71+ contactEvents: number;
72+ peakXiMm: number;
73+ peakDispMm: number;
74+ peakDispThickness: number;
75+ peakAccel: number;
76+ strikeRadiusM: number;
77+ pickupRadiusM: number;
78+ RsM: number;
79+ nsteps: number;
80+ simRate: number;
81+ nDropped: number;
82+ durationS: number;
83+ renderMs: number;
84+ oversample: number;
85+ nTheta: number;
86+ backend: 'webgpu' | 'cpu';
87+ warnings: string[];
88+ };
89+}
90+
91+// --- shared helpers ---
92+
93+export function hashMFiles(files: MFile[]): string {
94+ let h = 0x811c9dc5;
95+ for (const f of files) {
96+ const s = f.name + ' ' + f.source;
97+ for (let i = 0; i < s.length; i++) {
98+ h ^= s.charCodeAt(i);
99+ h = Math.imul(h, 0x01000193);
100+ }
101+ }
102+ return (h >>> 0).toString(36);
103+}
104+
105+function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
106+ const v = cache.get(key);
107+ if (v !== undefined) {
108+ cache.delete(key);
109+ cache.set(key, v);
110+ }
111+ return v;
112+}
113+
114+function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
115+ cache.delete(key);
116+ cache.set(key, value);
117+ while (cache.size > max) {
118+ const oldest = cache.keys().next().value as string;
119+ cache.delete(oldest);
120+ }
121+}
122+
123+type Params = Record<string, number | boolean>;
124+
125+// MATLAB source lines assigning parameters onto a struct P. bit_depth is a
126+// client-side WAV setting, not a model parameter.
127+function pLiteral(params: Params): string {
128+ const lines: string[] = [];
129+ for (const [k, v] of Object.entries(params)) {
130+ if (k === 'bit_depth') continue;
131+ if (typeof v === 'number') {
132+ if (!Number.isFinite(v)) throw new Error(`parameter ${k} is not finite`);
133+ lines.push(`P.${k} = ${v};`);
134+ } else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
135+ }
136+ return lines.join('\n');
137+}
138+
139+function keyFor(params: Params, keys: readonly string[]): string {
140+ return keys.map((k) => String(params[k] ?? 'd')).join('|');
141+}
142+
143+function knobValue(settings: Record<string, number>, key: string): number {
144+ const spec = CYMBAL2_KNOBS.find((k) => k.key === key)!;
145+ const raw = settings[key];
146+ return raw === undefined ? spec.default : Math.max(spec.min, Math.min(spec.max, raw));
147+}
148+
149+/** Velocity quantized to a few levels so identical strokes share a cached
150+ * render; 8 levels is inaudible next to the model's own dynamics. */
151+export function bucketVelocity(velocity01: number, levels = 8): number {
152+ return Math.ceil(Math.max(0.001, Math.min(1, velocity01)) * levels) / levels;
153+}
154+
155+// Files a render actually executes: a parse failure in one of these stops
156+// the render, while a typo in a CLI entry point must not.
157+const EXECUTED_BY_APP = new Set([
158+ 'cymbal2_grid.m',
159+ 'cymbal2_perm.m',
160+ 'cymbal2_system.m',
161+ 'cymbal2_mallet.m',
162+ 'cymbal2_sim.m',
163+ 'cymbal2_output.m',
164+ 'cymbal2_modes.m',
165+]);
166+
167+const MODES_KEYS = [
168+ 'radius_m', 'rise_m', 'thickness_m',
169+ 'youngs_modulus', 'poisson_ratio', 'density',
170+ 'washer_radius_m', 'washer_stiffness', 'support_clamped',
171+ 'n_radial', 'n_theta',
172+];
173+
174+// An error in the MATLAB source, carrying the file and line where numbl put
175+// it. Reported to the UI, which stops the render and shows nothing playable.
176+class CodeError extends Error {
177+ file?: string;
178+ line?: number;
179+ constructor(message: string, file?: string, line?: number) {
180+ super(message);
181+ this.name = 'CodeError';
182+ this.file = file;
183+ this.line = line;
184+ }
185+}
186+
187+const SKIP_WARNING = /skipping (\S+\.m) \((?:syntax error at line (\d+)|parse error)\)/;
188+
189+// numbl does not abort on a file it cannot parse: it drops that file with a
190+// warning and carries on, so a render could otherwise finish (or fail much
191+// later with a confusing 'undefined function') while the real fault is a
192+// typo in a file the stage never reached.
193+function skipWarningToError(warnings: string[]): CodeError | null {
194+ for (const w of warnings) {
195+ const hit = SKIP_WARNING.exec(w);
196+ if (!hit) continue;
197+ const file = hit[1];
198+ if (!EXECUTED_BY_APP.has(file)) continue;
199+ const line = hit[2] ? Number(hit[2]) : undefined;
200+ return new CodeError(
201+ line !== undefined
202+ ? 'syntax error; the render was stopped'
203+ : 'the file cannot be parsed; the render was stopped',
204+ file,
205+ line,
206+ );
207+ }
208+ return null;
209+}
210+
211+function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
212+ const warnings: string[] = [];
213+ const prevWarn = console.warn;
214+ console.warn = (...args: unknown[]) => {
215+ warnings.push(args.map(String).join(' '));
216+ prevWarn(...(args as []));
217+ };
218+ try {
219+ const r = executeCode(
220+ source,
221+ { initialVariableValues: inputs, displayResults: false },
222+ mfiles,
223+ 'stage.m',
224+ );
225+ const skipped = skipWarningToError(warnings);
226+ if (skipped) throw skipped;
227+ return r;
228+ } catch (err) {
229+ const skipped = skipWarningToError(warnings);
230+ if (skipped) throw skipped;
231+ throw err;
232+ } finally {
233+ console.warn = prevWarn;
234+ }
235+}
236+
237+// Parse every .m file before any expensive stage runs, so a syntax error is
238+// reported immediately instead of after a stage.
239+function checkMFiles(mfiles: MFile[]): void {
240+ runStage(mfiles, '1;');
241+}
242+
243+function asScalar(v: RuntimeValue): number {
244+ if (typeof v === 'number') return v;
245+ const t = v as RuntimeTensor;
246+ if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
247+ if (typeof v === 'boolean') return v ? 1 : 0;
248+ throw new Error('expected scalar result from numbl stage');
249+}
250+
251+function toF32(v: RuntimeValue): Float32Array {
252+ const t = v as RuntimeTensor;
253+ if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
254+ const out = new Float32Array(t.data.length);
255+ for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
256+ return out;
257+}
258+
259+// Max-pool the contact force history to about 2000 points for the plot.
260+function poolForce(fhist: Float32Array, simRate: number): { tpool: Float32Array; fpool: Float32Array } {
261+ const nsim = fhist.length;
262+ const np0 = Math.min(2000, nsim);
263+ const blk = Math.ceil(nsim / np0);
264+ const np = Math.ceil(nsim / blk);
265+ const fpool = new Float32Array(np);
266+ const tpool = new Float32Array(np);
267+ for (let b = 0; b < np; b++) {
268+ let m = 0;
269+ const end = Math.min(nsim, (b + 1) * blk);
270+ for (let i = b * blk; i < end; i++) m = Math.max(m, fhist[i]);
271+ fpool[b] = m;
272+ tpool[b] = (b * blk + blk / 2) / simRate;
273+ }
274+ return { tpool, fpool };
275+}
276+
277+// --- the simulation, on either backend ---
278+
279+// The stroke-schedule block from cymbal2_sim.m, verbatim, so the GPU render
280+// integrates exactly the schedule the interpreter would.
281+const SCHEDULE_SRC = `
282+K = max(1, round(P.strike_count));
283+tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);
284+vstr = zeros(1, K);
285+state = P.seed * 2654435761;
286+for i = 1:K
287+ state = mod(state * 1664525 + 1013904223, 4294967296);
288+ x = state / 4294967296;
289+ vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));
290+end
291+gap = P.strike_gap_m;
292+if gap <= 0
293+ gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;
294+end
295+travel = gap ./ vstr + vstr / (2 * P.strike_accel);
296+tlaunch = tstr - travel;
297+off = -min(0, min(tlaunch));
298+tlaunch = tlaunch + off;
299+duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
300+nsteps = floor(duration * P.sample_rate * P.oversample);
301+`;
302+
303+// Everything the response assembly needs, produced by either backend. The
304+// scalar model facts (rates, grid, thickness) come out of the setup stage
305+// rather than from a TS-side parameter struct, so the note-mapping path and
306+// the full-surface path share one implementation.
307+interface SimOutcome {
308+ y32: Float32Array;
309+ peakAccel: number;
310+ nDropped: number;
311+ tpool: Float32Array;
312+ fpool: Float32Array;
313+ backend: 'webgpu' | 'cpu';
314+ sampleRate: number;
315+ oversample: number;
316+ nTheta: number;
317+ normalize: boolean;
318+ d: {
319+ strokesPlayed: number;
320+ meanV: number;
321+ contactVel: number;
322+ peakForceN: number;
323+ contactMs: number;
324+ contactEvents: number;
325+ peakXiMm: number;
326+ peakDispM: number;
327+ peakDispThickness: number;
328+ nsteps: number;
329+ simRate: number;
330+ strikeRadiusM: number;
331+ pickupRadiusM: number;
332+ RsM: number;
333+ };
334+}
335+
336+// Scalars every backend extracts from P after the base source has run.
337+const P_FACTS = `
338+c_sr = P.sample_rate; c_os = P.oversample; c_nTheta = P.n_theta;
339+c_thick = P.thickness_m; c_norm = P.normalize;
340+`;
341+
342+async function renderOnGpu(
343+ pBase: string,
344+ mfiles: MFile[],
345+ env: RenderEnv,
346+ onProgress: (frac: number) => void,
347+): Promise<SimOutcome> {
348+ // A WebGPU device first: if there is none, fail before spending seconds on
349+ // the setup stage.
350+ const device = await requestGpuDevice();
351+ try {
352+ const r = runStage(mfiles, `
353+${pBase}
354+S = cymbal2_system(P);
355+MB = cymbal2_mallet(P);
356+G = S.G;
357+o_D1 = G.D1; o_D2 = G.D2;
358+o_F = G.F; o_Gi = G.Gi;
359+o_Rd = G.Rd; o_R2d = G.R2d;
360+o_AiB = S.AiB; o_AiC = S.AiC; o_Ai = S.Ai; o_Fi = S.Fi;
361+o_ASiB = MB.ASiB; o_ASiC = MB.ASiC;
362+o_vM = MB.vM; o_vR = MB.vR; o_vK = MB.vK; o_tipcol = MB.tipcol;
363+${SCHEDULE_SRC}
364+o_tlaunch = tlaunch; o_vstr = vstr;
365+c_nrt = G.nrt; c_Nr = G.Nr; c_nth = G.nth; c_hr = G.hr; c_hth = G.hth;
366+c_ht = S.ht; c_Eh = S.Eh; c_Rs = S.Rs; c_rhoh = S.rhoh;
367+c_fsc = S.forceScale; c_nlr0 = S.eomRow0; c_nDropped = S.nDropped;
368+c_ls = S.ls; c_js = S.js; c_lp = S.lp; c_jp = S.jp;
369+c_strikeR = S.strikeRadiusM; c_pickupR = S.pickupRadiusM;
370+c_tipi = MB.tipi; c_nfree = MB.nfree; c_ps = MB.ps;
371+c_holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / S.ht));
372+c_missDepth = gap + 0.06;
373+c_gap = gap; c_K = K; c_nsteps = nsteps;
374+c_simRate = P.sample_rate * P.oversample;
375+c_KM = P.head_stiffness; c_alM = P.head_alpha;
376+c_aH = P.strike_accel; c_aHr = P.retract_accel;
377+${P_FACTS}
378+`);
379+ const v = r.variableValues;
380+ const s = (n: string) => asScalar(v[n]);
381+ const t = (n: string) => toF32(v[n]);
382+ const setup: GpuSetup = {
383+ ops: {
384+ D1: t('o_D1'), D2: t('o_D2'),
385+ F: t('o_F'), Gi: t('o_Gi'),
386+ Rd: t('o_Rd'), R2d: t('o_R2d'),
387+ AiB: t('o_AiB'), AiC: t('o_AiC'), Ai: t('o_Ai'), Fi: t('o_Fi'),
388+ ASiB: t('o_ASiB'), ASiC: t('o_ASiC'),
389+ vM: t('o_vM'), vR: t('o_vR'), vK: t('o_vK'), tipcol: t('o_tipcol'),
390+ },
391+ tlaunch: Array.from((v.o_tlaunch as RuntimeTensor).data),
392+ vstr: Array.from((v.o_vstr as RuntimeTensor).data),
393+ nrt: s('c_nrt'), Nr: s('c_Nr'), nth: s('c_nth'),
394+ hr: s('c_hr'), hth: s('c_hth'),
395+ ht: s('c_ht'), Eh: s('c_Eh'), Rs: s('c_Rs'), rhoh: s('c_rhoh'),
396+ fsc: s('c_fsc'), nlr0: s('c_nlr0'),
397+ ls: s('c_ls'), js: s('c_js'), lp: s('c_lp'), jp: s('c_jp'),
398+ tipi: s('c_tipi'), nfree: s('c_nfree'), ps: s('c_ps'),
399+ KM: s('c_KM'), alM: s('c_alM'),
400+ aH: s('c_aH'), aHr: s('c_aHr'),
401+ holdSteps: s('c_holdSteps'), missDepth: s('c_missDepth'),
402+ gap: s('c_gap'), K: s('c_K'), nsteps: s('c_nsteps'),
403+ };
404+ const simRate = s('c_simRate');
405+ const thickness = s('c_thick');
406+
407+ const sim = await runGpuSim(device, setup, (done, n) => {
408+ env.checkCancel();
409+ onProgress(done / n);
410+ });
411+
412+ // The output stage through the unmodified cymbal2_output.m.
413+ const accT = RTV.tensor(Float64Array.from(sim.acc), [1, sim.acc.length]);
414+ const ro = runStage(mfiles, `
415+${pBase}
416+[y, peakAccel] = cymbal2_output(P, acc);
417+`, { acc: accT });
418+
419+ // Contact bookkeeping from the force history, as cymbal2_sim does it.
420+ let peakF = 0, contactSteps = 0, events = 0, prev = 0;
421+ for (let i = 0; i < sim.fhist.length; i++) {
422+ const f = sim.fhist[i];
423+ if (f > peakF) peakF = f;
424+ if (f > 0) { contactSteps++; if (prev === 0) events++; }
425+ prev = f;
426+ }
427+ const { tpool, fpool } = poolForce(sim.fhist, simRate);
428+ const meanV = setup.vstr.reduce((a, b) => a + b, 0) / setup.vstr.length;
429+
430+ return {
431+ y32: toF32(ro.variableValues.y),
432+ peakAccel: asScalar(ro.variableValues.peakAccel),
433+ nDropped: s('c_nDropped'),
434+ tpool,
435+ fpool,
436+ backend: 'webgpu',
437+ sampleRate: s('c_sr'),
438+ oversample: s('c_os'),
439+ nTheta: s('c_nTheta'),
440+ normalize: s('c_norm') !== 0,
441+ d: {
442+ strokesPlayed: sim.strokesPlayed,
443+ meanV,
444+ contactVel: sim.contactVel[0] ?? 0,
445+ peakForceN: peakF,
446+ contactMs: (1000 * contactSteps) / simRate,
447+ contactEvents: events,
448+ peakXiMm: sim.peakXi * 1000,
449+ peakDispM: sim.peakDisp,
450+ peakDispThickness: sim.peakDisp / thickness,
451+ nsteps: setup.nsteps,
452+ simRate,
453+ strikeRadiusM: s('c_strikeR'),
454+ pickupRadiusM: s('c_pickupR'),
455+ RsM: setup.Rs,
456+ },
457+ };
458+ } finally {
459+ device.destroy();
460+ }
461+}
462+
463+function renderOnCpu(pBase: string, mfiles: MFile[]): SimOutcome {
464+ const r = runStage(mfiles, `
465+${pBase}
466+S = cymbal2_system(P);
467+MB = cymbal2_mallet(P);
468+[acc, dinfo, fhist] = cymbal2_sim(P, S, MB);
469+[y, peakAccel] = cymbal2_output(P, acc);
470+% max-pool the contact force history to about 2000 points for the plot
471+nsim = length(fhist);
472+np = min(2000, nsim);
473+blk = ceil(nsim / np);
474+np = ceil(nsim / blk);
475+fpool = max(reshape([fhist, zeros(1, np * blk - nsim)], blk, np), [], 1);
476+tpool = ((0:np - 1) * blk + blk / 2) / dinfo.simRate;
477+d_strokes = dinfo.strokesPlayed;
478+d_meanV = dinfo.meanV;
479+d_cvel = dinfo.contactVel(1);
480+d_peakF = dinfo.peakForceN;
481+d_contactMs = dinfo.contactMs;
482+d_events = dinfo.contactEvents;
483+d_xi = dinfo.peakXiMm;
484+d_disp = dinfo.peakDispM;
485+d_dispH = dinfo.peakDispThickness;
486+d_nsteps = dinfo.nsteps;
487+d_simRate = dinfo.simRate;
488+d_dropped = S.nDropped;
489+d_strikeR = S.strikeRadiusM;
490+d_pickupR = S.pickupRadiusM;
491+d_Rs = S.Rs;
492+${P_FACTS}
493+`);
494+ const v = r.variableValues;
495+ return {
496+ y32: toF32(v.y),
497+ peakAccel: asScalar(v.peakAccel),
498+ nDropped: asScalar(v.d_dropped),
499+ tpool: toF32(v.tpool),
500+ fpool: toF32(v.fpool),
501+ backend: 'cpu',
502+ sampleRate: asScalar(v.c_sr),
503+ oversample: asScalar(v.c_os),
504+ nTheta: asScalar(v.c_nTheta),
505+ normalize: asScalar(v.c_norm) !== 0,
506+ d: {
507+ strokesPlayed: asScalar(v.d_strokes),
508+ meanV: asScalar(v.d_meanV),
509+ contactVel: asScalar(v.d_cvel),
510+ peakForceN: asScalar(v.d_peakF),
511+ contactMs: asScalar(v.d_contactMs),
512+ contactEvents: asScalar(v.d_events),
513+ peakXiMm: asScalar(v.d_xi),
514+ peakDispM: asScalar(v.d_disp),
515+ peakDispThickness: asScalar(v.d_dispH),
516+ nsteps: asScalar(v.d_nsteps),
517+ simRate: asScalar(v.d_simRate),
518+ strikeRadiusM: asScalar(v.d_strikeR),
519+ pickupRadiusM: asScalar(v.d_pickupR),
520+ RsM: asScalar(v.d_Rs),
521+ },
522+ };
523+}
524+
525+// The WGSL mirrors the shipped cymbal2_sim.m, so an edited step falls back
526+// to the interpreter, where the edit is live. Edits to the setup and output
527+// stages flow through numbl on both paths.
528+function simEdited(env: RenderEnv): boolean {
529+ const current = env.mfiles.find((f) => f.name === 'cymbal2_sim.m')?.source;
530+ const shipped = env.shipped.find((f) => f.name === 'cymbal2_sim.m')?.source;
531+ return current !== shipped;
532+}
533+
534+async function renderSim(
535+ pBase: string,
536+ env: RenderEnv,
537+ edited: boolean,
538+ onIntegrateProgress: (frac: number) => void,
539+ warnings: string[],
540+): Promise<SimOutcome> {
541+ if (!edited) {
542+ try {
543+ return await renderOnGpu(pBase, env.mfiles, env, onIntegrateProgress);
544+ } catch (err) {
545+ if (err instanceof RenderCancelled) throw err;
546+ const detail = err instanceof Error ? err.message : String(err);
547+ // An error in the MATLAB source (it carries a file, or names the
548+ // divergence guard) is a fault to report, not a reason to retry the
549+ // same code on the CPU. Only a genuine GPU problem falls back.
550+ const fromCode =
551+ err instanceof CodeError ||
552+ typeof (err as { file?: unknown }).file === 'string' ||
553+ (err as { name?: string }).name === 'RuntimeError';
554+ if (fromCode || /diverged/.test(detail)) throw err;
555+ // console.warn is intercepted during a render; use info so the reason
556+ // is visible in the devtools console.
557+ console.info(`cymbal2: WebGPU unavailable (${detail}); rendering on the CPU`);
558+ }
559+ }
560+ const out = renderOnCpu(pBase, env.mfiles);
561+ if (edited) {
562+ // Not a fault: just say why the render took the slow path.
563+ warnings.push('cymbal2_sim.m differs from the shipped source, so the render used the interpreter (the GPU path mirrors the shipped step)');
564+ }
565+ return out;
566+}
567+
568+// --- the explorer path (full params, with diagnostics) ---
569+
570+interface ModesEntry {
571+ modeF: Float32Array;
572+ modeM: Float32Array;
573+}
574+
575+const modesCache = new Map<string, ModesEntry>();
576+const MODES_CACHE_MAX = 4;
577+
578+// Finished per-note audio: an in-memory LRU in front of an IndexedDB
579+// store, so renders survive page reloads (budgets and eviction live in
580+// noteCache.ts).
581+const noteCache = new NoteCache('cymbal2-instrument-notes');
582+
583+async function renderExplore(
584+ req: TrackRenderRequest,
585+ env: RenderEnv,
586+): Promise<{ audio: Float32Array; diag: Cymbal2ExploreDiag }> {
587+ const quality = req.quality;
588+ const params: Params = { ...(req.params ?? {}) };
589+ if (quality === 'preview') {
590+ // oversample 1, everything else as given: the render covers the full
591+ // duration, so the decay can be judged, and the azimuthal grid is kept
592+ // because hard strokes need the full wavenumber band to carry the
593+ // cascade, and coarsening it destabilizes the integration.
594+ params.oversample = 1;
595+ }
596+ const pBase = `P = cymbal2_defaults();\n${pLiteral(params)}`;
597+
598+ // Cheap, and it must come first: a file that does not parse stops the
599+ // render here rather than partway through a stage.
600+ env.post(0.01, 'parsing');
601+ await env.yieldToHost();
602+ env.checkCancel();
603+ checkMFiles(env.mfiles);
604+
605+ const codeKey = hashMFiles(env.mfiles);
606+
607+ // ---- linear modes, cheap, cached on geometry ----
608+ const modesKey = codeKey + '|' + keyFor(params, MODES_KEYS);
609+ let modesEntry = lruGet(modesCache, modesKey);
610+ if (!modesEntry) {
611+ env.post(0.05, 'modes');
612+ await env.yieldToHost();
613+ env.checkCancel();
614+ const r = runStage(env.mfiles, `
615+${pBase}
616+mF = [];
617+mM = [];
618+for mw = 0:min(10, floor(P.n_theta / 2))
619+ f = cymbal2_modes(P, mw);
620+ f = f(f > 0.5 & f < 4000);
621+ nk = min(4, length(f));
622+ mF = [mF; f(1:nk)];
623+ mM = [mM; mw * ones(nk, 1)];
624+end
625+`);
626+ modesEntry = {
627+ modeF: toF32(r.variableValues.mF),
628+ modeM: toF32(r.variableValues.mM),
629+ };
630+ lruSet(modesCache, modesKey, modesEntry, MODES_CACHE_MAX);
631+ }
632+
633+ // ---- the time integration ----
634+ env.post(0.1, 'integrate');
635+ await env.yieldToHost();
636+ env.checkCancel();
637+ const t0 = performance.now();
638+ const warnings: string[] = [];
639+ const out = await renderSim(pBase, env, simEdited(env), (frac) =>
640+ env.post(0.1 + 0.85 * frac, 'integrate'), warnings);
641+ const renderMs = performance.now() - t0;
642+
643+ env.post(0.97, 'output');
644+ const { y32, peakAccel, nDropped } = out;
645+
646+ if (nDropped > 0) {
647+ warnings.push(
648+ `${nDropped} azimuthal columns were truncated: the radial grid cannot resolve their inner boundary layer; raise n_radial`,
649+ );
650+ }
651+ if (!out.normalize) {
652+ let peak = 0;
653+ for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
654+ if (peak > 1.0) {
655+ warnings.push(
656+ `clipping: peak acceleration ${peakAccel.toExponential(2)} m/s² exceeds full_scale_accel`,
657+ );
658+ }
659+ }
660+
661+ const diag: Cymbal2ExploreDiag = {
662+ quality,
663+ forceT: out.tpool,
664+ forceN: out.fpool,
665+ modeF: modesEntry.modeF.slice(),
666+ modeM: modesEntry.modeM.slice(),
667+ diag: {
668+ strokesPlayed: out.d.strokesPlayed,
669+ meanV: out.d.meanV,
670+ contactVel: out.d.contactVel,
671+ peakForceN: out.d.peakForceN,
672+ contactMs: out.d.contactMs,
673+ contactEvents: out.d.contactEvents,
674+ peakXiMm: out.d.peakXiMm,
675+ peakDispMm: out.d.peakDispM * 1000,
676+ peakDispThickness: out.d.peakDispThickness,
677+ peakAccel,
678+ strikeRadiusM: out.d.strikeRadiusM,
679+ pickupRadiusM: out.d.pickupRadiusM,
680+ RsM: out.d.RsM,
681+ nsteps: out.d.nsteps,
682+ simRate: out.d.simRate,
683+ nDropped,
684+ durationS: y32.length / out.sampleRate,
685+ renderMs,
686+ oversample: out.oversample,
687+ nTheta: out.nTheta,
688+ backend: out.backend,
689+ warnings,
690+ },
691+ };
692+ return { audio: y32, diag };
693+}
694+
695+// --- the track path (midiboot): per-stroke rendering with caching ---
696+
697+function fadeTail(buf: Float32Array, sampleRate: number): void {
698+ const n = Math.min(buf.length, Math.round(0.25 * sampleRate));
699+ if (n < 2) return;
700+ for (let i = 0; i < n; i++) {
701+ buf[buf.length - n + i] *= 0.5 * (1 + Math.cos((Math.PI * (i + 1)) / n));
702+ }
703+}
704+
705+function notePLines(note: { pitch: number; velocity: number }, req: TrackRenderRequest, previewTailS: number): string {
706+ // The note-to-parameter mapping: MIDI pitch scales the shell geometry,
707+ // every length by s = 2^-((pitch - 60)/12), which shifts the whole
708+ // spectrum by the semitone ratio exactly (shell modes go as
709+ // thickness / radius^2, and the curvature and washer frequencies scale
710+ // the same way when the washer circle scales too). The mallet stays
711+ // fixed, so this is one player striking a family of cymbal sizes; pitch
712+ // 60 is the default 18-inch shell, playable 48..72. The sound is
713+ // inharmonic, so the note number is a handle on register, not a tuned
714+ // fundamental. Velocity drives the hand speed of the stroke, and the
715+ // crash character follows from the physics: soft notes stay nearly
716+ // linear, hard notes wash. Fixed: a single stroke, a fixed seed,
717+ // normalize off so velocity survives. The preview tier drops the
718+ // oversampling and shortens the tail but keeps the azimuthal grid; how
719+ // short the tail is depends on the backend (the caller probes for
720+ // WebGPU), since a longer ring costs seconds on the GPU but minutes on
721+ // the interpreter. Optional request params land after the mapping,
722+ // overriding it, which is what the node tests use to shrink the render.
723+ const settings = req.settings ?? {};
724+ const s = Math.pow(2, -(Math.max(48, Math.min(72, Math.round(note.pitch))) - 60) / 12);
725+ const lines = [
726+ 'P = cymbal2_defaults();',
727+ `P.sample_rate = ${req.sampleRate};`,
728+ 'P.strike_count = 1;',
729+ 'P.seed = 1;',
730+ 'P.pre_delay_ms = 0;',
731+ 'P.normalize = 0;',
732+ ];
733+ if (req.quality === 'preview') lines.push('P.oversample = 1;', `P.tail_s = ${previewTailS};`);
734+ lines.push(
735+ `P.radius_m = P.radius_m * ${s};`,
736+ `P.rise_m = P.rise_m * ${s};`,
737+ `P.thickness_m = P.thickness_m * ${s};`,
738+ `P.washer_radius_m = P.washer_radius_m * ${s};`,
739+ `P.strike_velocity = ${0.3 + (6.5 - 0.3) * note.velocity};`,
740+ `P.strike_r = ${knobValue(settings, 'strike_r')};`,
741+ `P.grip_spring = ${knobValue(settings, 'grip_spring')};`,
742+ );
743+ if (req.params) lines.push(pLiteral(req.params));
744+ return lines.join('\n');
745+}
746+
747+async function renderTrackNotes(
748+ req: TrackRenderRequest,
749+ env: RenderEnv,
750+): Promise<{ audio: Float32Array }> {
751+ if (req.notes.length === 0) return { audio: new Float32Array(req.lengthSamples) };
752+
753+ env.post(0.005, 'parsing');
754+ await env.yieldToHost();
755+ env.checkCancel();
756+ checkMFiles(env.mfiles);
757+
758+ const codeKey = hashMFiles(env.mfiles);
759+ const settings = req.settings ?? {};
760+ const knobsKey = CYMBAL2_KNOBS.map((k) => `${k.key}=${knobValue(settings, k.key)}`).join(',');
761+ const overridesKey = req.params ? JSON.stringify(req.params) : '';
762+ const edited = simEdited(env);
763+
764+ // The preview tail is a compromise with the backend: the cymbal's low
765+ // modes carry a T60 of about 5 s, and on the GPU a 4 s ring costs a few
766+ // seconds, but on the interpreter it costs minutes. Probe once per
767+ // request and choose accordingly; the tail is part of the cache key, so
768+ // a browser gaining or losing WebGPU re-renders rather than reusing the
769+ // other tier's tail.
770+ let previewTailS = 1.5;
771+ if (req.quality === 'preview' && !edited) {
772+ try {
773+ const probe = await requestGpuDevice();
774+ probe.destroy();
775+ previewTailS = 4;
776+ } catch { /* no WebGPU; keep the short tail */ }
777+ }
778+
779+ // Group notes by everything that affects their sound. A crash ignores
780+ // note duration (the ring is the physics), so a repeated stroke costs
781+ // one integration per distinct pitch and velocity level.
782+ const groups = new Map<string, { pitch: number; velocity: number; onsets: number[] }>();
783+ for (const n of req.notes) {
784+ const pitch = Math.round(n.pitch);
785+ const vel = bucketVelocity(n.velocity);
786+ const key = `${codeKey}|${req.quality}|${previewTailS}|${req.sampleRate}|${knobsKey}|${overridesKey}|${pitch}|${vel}`;
787+ let g = groups.get(key);
788+ if (!g) groups.set(key, (g = { pitch, velocity: vel, onsets: [] }));
789+ g.onsets.push(n.timeSec);
790+ }
791+
792+ const placed: { buf: Float32Array; onsets: number[] }[] = [];
793+ // Resolve cache hits first (memory, then IndexedDB), so the progress
794+ // denominator counts only the sounds that actually need rendering.
795+ const resolved = new Map<string, Float32Array>();
796+ for (const key of groups.keys()) {
797+ const hit = await noteCache.get(key);
798+ if (hit) resolved.set(key, hit);
799+ }
800+ const misses = groups.size - resolved.size;
801+ let done = 0;
802+ for (const [key, g] of groups) {
803+ let buf = resolved.get(key);
804+ if (!buf) {
805+ const i = done;
806+ const label = `note ${i + 1}/${misses}: integrate`;
807+ env.post(i / Math.max(1, misses), label);
808+ await env.yieldToHost();
809+ env.checkCancel();
810+ const warnings: string[] = [];
811+ const sim = await renderSim(notePLines(g, req, previewTailS), env, edited, (frac) =>
812+ env.post((i + 0.05 + 0.9 * frac) / Math.max(1, misses), label), warnings);
813+ buf = sim.y32;
814+ // The integration stops at tail_s, which for a cymbal is well before
815+ // silence; a short raised-cosine fade removes the cliff.
816+ fadeTail(buf, req.sampleRate);
817+ noteCache.put(key, buf);
818+ done++;
819+ }
820+ placed.push({ buf, onsets: g.onsets });
821+ }
822+
823+ // lengthSamples is a minimum: size the output to the longest tail so a
824+ // ring that outlives the last note is returned in full, and place every
825+ // copy.
826+ let length = req.lengthSamples;
827+ for (const p of placed) {
828+ for (const t of p.onsets) length = Math.max(length, Math.round(t * req.sampleRate) + p.buf.length);
829+ }
830+ const out = new Float32Array(length);
831+ for (const p of placed) {
832+ for (const t of p.onsets) {
833+ const off = Math.round(t * req.sampleRate);
834+ const end = Math.min(out.length, off + p.buf.length);
835+ for (let j = Math.max(0, off); j < end; j++) out[j] += p.buf[j - off];
836+ }
837+ }
838+ return { audio: out };
839+}
840+
841+/** The single entry point: explorer path when diagnostics are requested,
842+ * the per-stroke track path otherwise. */
843+export async function renderRequest(
844+ variant: Cymbal2Variant,
845+ req: TrackRenderRequest,
846+ env: RenderEnv,
847+): Promise<{ audio: Float32Array; diag?: Cymbal2ExploreDiag }> {
848+ void variant;
849+ if (req.wantDiagnostics) return renderExplore(req, env);
850+ return renderTrackNotes(req, env);
851+}
webapp/src/worker/gpuSim.ts →instrument/gpuSim.tsrenamed+0−0View file
No changes to the file's content.
instrument/index.tsadded+116−0View file
@@ -0,0 +1,116 @@
1+// The cymbal (Samejima 2021 nonlinear FD shell) behind the generic
2+// track-render instrument interface, for direct import by midiboot and by
3+// the cymbal2 web app. The MATLAB model in ../matlab is executed in a
4+// lazily spawned worker that this module owns: setup and output run
5+// through numbl, and the time loop runs on WebGPU when the browser
6+// exposes it (a WGSL port of cymbal2_sim.m), falling back to the numbl
7+// interpreter otherwise. The worker is kept alive across renders so its
8+// caches (the mode tables and finished stroke renders) persist. Aborting
9+// a request is cooperative, checked between stages and between GPU
10+// chunks; a CPU integration cannot be interrupted mid-call, so there the
11+// grace period elapses and the worker is torn down (and the caches with
12+// it) and respawned.
13+import type { Instrument, TrackRenderRequest, TrackRenderResult } from './types.ts';
14+import { CYMBAL2_KNOBS, type Cymbal2Variant } from './knobs.ts';
15+import type { WorkerCancelMsg, WorkerRenderMsg, WorkerResponse } from './worker.ts';
16+
17+const CANCEL_GRACE_MS = 5000;
18+
19+interface Pending {
20+ resolve: (r: TrackRenderResult) => void;
21+ reject: (e: Error) => void;
22+ onProgress?: (p: { frac: number; stage: string }) => void;
23+}
24+
25+let worker: Worker | null = null;
26+let nextId = 1;
27+const pending = new Map<number, Pending>();
28+
29+function ensureWorker(): Worker {
30+ if (worker) return worker;
31+ worker = new Worker(new URL('./worker.ts', import.meta.url), { type: 'module' });
32+ worker.onmessage = (ev: MessageEvent<WorkerResponse>) => {
33+ const msg = ev.data;
34+ const p = pending.get(msg.id);
35+ if (!p) return;
36+ if (msg.kind === 'progress') {
37+ p.onProgress?.({ frac: msg.frac, stage: msg.stage });
38+ return;
39+ }
40+ pending.delete(msg.id);
41+ if (msg.kind === 'result') p.resolve({ audio: msg.audio, diag: msg.diag });
42+ else if (msg.kind === 'cancelled') p.reject(new DOMException('render cancelled', 'AbortError'));
43+ else p.reject(Object.assign(new Error(msg.message), { file: msg.file, line: msg.line }));
44+ };
45+ worker.onerror = (ev) => {
46+ const err = new Error(`cymbal2 instrument worker failed: ${ev.message}`);
47+ for (const p of pending.values()) p.reject(err);
48+ pending.clear();
49+ };
50+ return worker;
51+}
52+
53+function teardownWorker(reason: Error): void {
54+ if (worker) worker.terminate();
55+ worker = null;
56+ for (const p of pending.values()) p.reject(reason);
57+ pending.clear();
58+}
59+
60+function renderOnWorker(variant: Cymbal2Variant, req: TrackRenderRequest): Promise<TrackRenderResult> {
61+ if (req.signal?.aborted) return Promise.reject(new DOMException('render cancelled', 'AbortError'));
62+ const w = ensureWorker();
63+ const id = nextId++;
64+ const msg: WorkerRenderMsg = {
65+ kind: 'render',
66+ id,
67+ variant,
68+ notes: req.notes.map((n) => ({ timeSec: n.timeSec, durationSec: n.durationSec, pitch: n.pitch, velocity: n.velocity })),
69+ sampleRate: req.sampleRate,
70+ lengthSamples: req.lengthSamples,
71+ quality: req.quality,
72+ settings: req.settings ?? {},
73+ params: req.params,
74+ mfiles: req.mfiles,
75+ wantDiagnostics: req.wantDiagnostics,
76+ };
77+ return new Promise<TrackRenderResult>((resolve, reject) => {
78+ let graceTimer: ReturnType<typeof setTimeout> | null = null;
79+ const onAbort = () => {
80+ const cancel: WorkerCancelMsg = { kind: 'cancel', id };
81+ worker?.postMessage(cancel);
82+ // Escalation: a stage that never yields (or a hung worker) is torn
83+ // down after the grace period, losing the caches.
84+ graceTimer = setTimeout(() => {
85+ if (pending.has(id)) teardownWorker(new DOMException('render cancelled', 'AbortError'));
86+ }, CANCEL_GRACE_MS);
87+ };
88+ req.signal?.addEventListener('abort', onAbort, { once: true });
89+ pending.set(id, {
90+ resolve: (r) => {
91+ req.signal?.removeEventListener('abort', onAbort);
92+ if (graceTimer !== null) clearTimeout(graceTimer);
93+ resolve(r);
94+ },
95+ reject: (e) => {
96+ req.signal?.removeEventListener('abort', onAbort);
97+ if (graceTimer !== null) clearTimeout(graceTimer);
98+ reject(e);
99+ },
100+ onProgress: req.onProgress,
101+ });
102+ w.postMessage(msg);
103+ });
104+}
105+
106+class Cymbal2Instrument implements Instrument {
107+ readonly id = 'cymbal2';
108+ readonly name = 'Cymbal2';
109+ readonly knobs = CYMBAL2_KNOBS;
110+ renderTrack(req: TrackRenderRequest): Promise<TrackRenderResult> {
111+ return renderOnWorker('cymbal2', req);
112+ }
113+}
114+
115+export const cymbal2Instrument: Instrument = new Cymbal2Instrument();
116+export const cymbal2Instruments: Instrument[] = [cymbal2Instrument];
instrument/knobs.tsadded+10−0View file
@@ -0,0 +1,10 @@
1+// Curated knob declarations, in their own module so the app-facing
2+// index.ts can import them without pulling numbl into the main bundle.
3+import type { KnobSpec } from './types.ts';
4+
5+export const CYMBAL2_KNOBS: KnobSpec[] = [
6+ { key: 'strike_r', label: 'strike radius', min: 0.05, max: 0.99, default: 0.6 },
7+ { key: 'grip_spring', label: 'grip stiffness (N m/rad)', min: 20, max: 2000, default: 140, scale: 'log' },
8+];
9+
10+export type Cymbal2Variant = 'cymbal2';
instrument/mfiles.tsadded+14−0View file
@@ -0,0 +1,14 @@
1+// The shipped MATLAB sources, loaded from ../matlab at build time. Browser
2+// only (import.meta.glob is a Vite feature); node tests read the same files
3+// through fs instead.
4+import type { MFile } from './types.ts';
5+
6+const modules = import.meta.glob('../matlab/*.m', {
7+ eager: true,
8+ query: '?raw',
9+ import: 'default',
10+}) as Record<string, string>;
11+
12+export const DEFAULT_MFILES: MFile[] = Object.entries(modules)
13+ .map(([path, source]) => ({ name: path.split('/').pop()!, source }))
14+ .sort((a, b) => a.name.localeCompare(b.name));
instrument/noteCache.tsadded+180−0View file
@@ -0,0 +1,180 @@
1+// Persistent cache for finished note renders: an in-memory LRU in front of
2+// an IndexedDB store, so a rendered note survives page reloads, worker
3+// teardowns, and browser restarts. Keys already carry everything that
4+// identifies a sound (the .m code hash, quality tier, sample rate, knob
5+// values, pitch, bucketed velocity/duration), so a stale entry from older
6+// code is simply never asked for again and ages out.
7+//
8+// Eviction is by budget and by age. The memory tier is bounded by total
9+// samples and evicts least-recently-used, as before. The IndexedDB tier
10+// has its own larger sample budget, trimmed least-recently-used after
11+// each write, and entries unused for 30 days are dropped when the store
12+// is opened. Sample data and bookkeeping live in separate object stores
13+// ('data' and 'meta'), so the eviction scan never loads audio.
14+//
15+// Everything degrades to memory-only: no indexedDB global (the node
16+// tests), a blocked or failing database, or a private-mode quota all just
17+// disable persistence.
18+
19+interface MetaRec {
20+ key: string;
21+ sampleCount: number;
22+ lastUsedAt: number;
23+}
24+
25+interface DataRec {
26+ key: string;
27+ samples: Float32Array;
28+}
29+
30+const MEM_MAX_SAMPLES = 48e6; // ~192 MB of Float32 renders in memory
31+const DB_MAX_SAMPLES = 96e6; // ~384 MB on disk
32+const MAX_AGE_MS = 30 * 24 * 3600e3;
33+
34+function reqDone<T>(req: IDBRequest<T>): Promise<T> {
35+ return new Promise((resolve, reject) => {
36+ req.onsuccess = () => resolve(req.result);
37+ req.onerror = () => reject(req.error);
38+ });
39+}
40+
41+export class NoteCache {
42+ private mem = new Map<string, Float32Array>();
43+ private memSamples = 0;
44+ private dbPromise: Promise<IDBDatabase | null> | null = null;
45+ /** Bookkeeping mirror of the 'meta' store, loaded once at open. */
46+ private meta = new Map<string, MetaRec>();
47+ private dbSamples = 0;
48+
49+ constructor(private readonly dbName: string) {}
50+
51+ private open(): Promise<IDBDatabase | null> {
52+ if (this.dbPromise) return this.dbPromise;
53+ this.dbPromise = (async () => {
54+ if (typeof indexedDB === 'undefined') return null;
55+ try {
56+ const req = indexedDB.open(this.dbName, 1);
57+ req.onupgradeneeded = () => {
58+ const db = req.result;
59+ if (!db.objectStoreNames.contains('data')) db.createObjectStore('data', { keyPath: 'key' });
60+ if (!db.objectStoreNames.contains('meta')) db.createObjectStore('meta', { keyPath: 'key' });
61+ };
62+ const db = await new Promise<IDBDatabase>((resolve, reject) => {
63+ req.onsuccess = () => resolve(req.result);
64+ req.onerror = () => reject(req.error);
65+ req.onblocked = () => reject(new Error('blocked'));
66+ });
67+ db.onversionchange = () => db.close();
68+
69+ // Load the bookkeeping and drop entries unused for too long.
70+ const metas = (await reqDone(
71+ db.transaction('meta', 'readonly').objectStore('meta').getAll(),
72+ )) as MetaRec[];
73+ const now = Date.now();
74+ const stale: string[] = [];
75+ for (const m of metas) {
76+ if (now - m.lastUsedAt > MAX_AGE_MS) stale.push(m.key);
77+ else {
78+ this.meta.set(m.key, m);
79+ this.dbSamples += m.sampleCount;
80+ }
81+ }
82+ if (stale.length > 0) {
83+ const tx = db.transaction(['data', 'meta'], 'readwrite');
84+ for (const key of stale) {
85+ tx.objectStore('data').delete(key);
86+ tx.objectStore('meta').delete(key);
87+ }
88+ }
89+ return db;
90+ } catch {
91+ return null; // persistence unavailable; memory tier still works
92+ }
93+ })();
94+ return this.dbPromise;
95+ }
96+
97+ private memPut(key: string, samples: Float32Array): void {
98+ const old = this.mem.get(key);
99+ if (old) this.memSamples -= old.length;
100+ this.mem.delete(key);
101+ this.mem.set(key, samples);
102+ this.memSamples += samples.length;
103+ while (this.memSamples > MEM_MAX_SAMPLES && this.mem.size > 1) {
104+ const oldest = this.mem.keys().next().value as string;
105+ this.memSamples -= this.mem.get(oldest)!.length;
106+ this.mem.delete(oldest);
107+ }
108+ }
109+
110+ /** Record a use; fire-and-forget so hits never wait on a write. */
111+ private touch(key: string): void {
112+ const m = this.meta.get(key);
113+ if (!m) return;
114+ m.lastUsedAt = Date.now();
115+ void this.open().then((db) => {
116+ if (!db) return;
117+ try {
118+ db.transaction('meta', 'readwrite').objectStore('meta').put(m);
119+ } catch { /* persistence unavailable */ }
120+ });
121+ }
122+
123+ async get(key: string): Promise<Float32Array | undefined> {
124+ const hit = this.mem.get(key);
125+ if (hit) {
126+ this.mem.delete(key);
127+ this.mem.set(key, hit); // LRU refresh
128+ this.touch(key);
129+ return hit;
130+ }
131+ const db = await this.open();
132+ if (!db) return undefined;
133+ try {
134+ const rec = (await reqDone(
135+ db.transaction('data', 'readonly').objectStore('data').get(key),
136+ )) as DataRec | undefined;
137+ if (!rec) return undefined;
138+ const samples = rec.samples instanceof Float32Array ? rec.samples : new Float32Array(rec.samples as ArrayLike<number>);
139+ this.memPut(key, samples);
140+ this.touch(key);
141+ return samples;
142+ } catch {
143+ return undefined;
144+ }
145+ }
146+
147+ /** Store a fresh render: memory now, disk in the background. */
148+ put(key: string, samples: Float32Array): void {
149+ this.memPut(key, samples);
150+ void this.persist(key, samples);
151+ }
152+
153+ private async persist(key: string, samples: Float32Array): Promise<void> {
154+ const db = await this.open();
155+ if (!db) return;
156+ try {
157+ const prev = this.meta.get(key);
158+ if (prev) this.dbSamples -= prev.sampleCount;
159+ const m: MetaRec = { key, sampleCount: samples.length, lastUsedAt: Date.now() };
160+ const tx = db.transaction(['data', 'meta'], 'readwrite');
161+ tx.objectStore('data').put({ key, samples } satisfies DataRec);
162+ tx.objectStore('meta').put(m);
163+ this.meta.set(key, m);
164+ this.dbSamples += m.sampleCount;
165+
166+ // Trim least-recently-used until the disk tier is under budget.
167+ if (this.dbSamples > DB_MAX_SAMPLES) {
168+ const byAge = [...this.meta.values()].sort((a, b) => a.lastUsedAt - b.lastUsedAt);
169+ for (const victim of byAge) {
170+ if (this.dbSamples <= DB_MAX_SAMPLES) break;
171+ if (victim.key === key) continue; // never evict what was just written
172+ tx.objectStore('data').delete(victim.key);
173+ tx.objectStore('meta').delete(victim.key);
174+ this.meta.delete(victim.key);
175+ this.dbSamples -= victim.sampleCount;
176+ }
177+ }
178+ } catch { /* persistence unavailable; the memory tier holds it */ }
179+ }
180+}
instrument/package-lock.jsonadded+2309−0View file
This diff is 2,314 lines long and is not shown.
instrument/package.jsonadded+17−0View file
@@ -0,0 +1,17 @@
1+{
2+ "name": "cymbal2-instrument",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "description": "The cymbal2 model (Samejima 2021 FD shell) behind the track-render instrument interface",
7+ "scripts": {
8+ "test": "tsx test.ts"
9+ },
10+ "dependencies": {
11+ "numbl": "^0.4.19"
12+ },
13+ "devDependencies": {
14+ "@webgpu/types": "^0.1.72",
15+ "tsx": "^4.19.0"
16+ }
17+}
instrument/test.tsadded+104−0View file
@@ -0,0 +1,104 @@
1+// Node-side test of the cymbal2 track-render instrument: executes the
2+// actual MATLAB model through numbl on the CPU path (no browser, no
3+// worker, no GPU) and checks the audio. The request's params override
4+// shrinks the render (low rate, short tail) so the full time integration
5+// stays test-sized; the GPU path has its own harness in
6+// webapp/scripts/gpu-selftest.ts. Run with: npx tsx test.ts
7+import { readdirSync, readFileSync } from 'node:fs';
8+import { join, dirname } from 'node:path';
9+import { fileURLToPath } from 'node:url';
10+import { renderRequest, type RenderEnv } from './core.ts';
11+import type { NoteSpec, TrackRenderRequest } from './types.ts';
12+
13+const matlabDir = join(dirname(fileURLToPath(import.meta.url)), '..', 'matlab');
14+const mfiles = readdirSync(matlabDir)
15+ .filter((f) => f.endsWith('.m'))
16+ .sort()
17+ .map((f) => ({ name: f, source: readFileSync(join(matlabDir, f), 'utf8') }));
18+
19+const env: RenderEnv = {
20+ mfiles,
21+ shipped: mfiles,
22+ post: () => {},
23+ checkCancel: () => {},
24+ yieldToHost: async () => {},
25+};
26+
27+let failures = 0;
28+function check(label: string, ok: boolean, detail = ''): void {
29+ console.log(`${ok ? 'ok ' : 'FAIL'} ${label}${detail ? ` (${detail})` : ''}`);
30+ if (!ok) failures++;
31+}
32+
33+function stats(y: Float32Array) {
34+ let peak = 0;
35+ let finite = true;
36+ for (let i = 0; i < y.length; i++) {
37+ if (!Number.isFinite(y[i])) finite = false;
38+ peak = Math.max(peak, Math.abs(y[i]));
39+ }
40+ return { peak, finite };
41+}
42+
43+// The explicit nonlinearity needs the full simulation rate: dropping the
44+// sample rate raises the time step and the integration diverges. Shrink
45+// the render by shortening the tail instead.
46+const SR = 44100;
47+
48+function req(notes: NoteSpec[], overrides: Partial<TrackRenderRequest> = {}): TrackRenderRequest {
49+ const end = notes.reduce((m, n) => Math.max(m, n.timeSec + n.durationSec), 0);
50+ return {
51+ notes,
52+ sampleRate: SR,
53+ lengthSamples: Math.ceil((end + 0.4) * SR),
54+ quality: 'preview',
55+ settings: {},
56+ params: { tail_s: 0.25 }, // test-sized ring; the physics is the same
57+ ...overrides,
58+ };
59+}
60+
61+const note = (timeSec: number, pitch: number, velocity: number): NoteSpec =>
62+ ({ timeSec, durationSec: 0.1, pitch, velocity });
63+
64+const t0 = performance.now();
65+
66+// One stroke: finite and audible. (The renders are cached by note key, so
67+// an identical repeat returns the same buffer; determinism of the model
68+// itself is pinned in matlab/cymbal2_validate.m.)
69+const a = await renderRequest('cymbal2', req([note(0, 60, 0.4)]), env);
70+{
71+ const s = stats(a.audio);
72+ check('single stroke is finite', s.finite);
73+ check('single stroke is audible', s.peak > 1e-6, `peak ${s.peak.toExponential(2)}`);
74+}
75+
76+// The note cache.
77+{
78+ const t = performance.now();
79+ await renderRequest('cymbal2', req([note(0, 60, 0.4)]), env);
80+ const ms = performance.now() - t;
81+ check('repeat render hits the note cache', ms < 500, `${ms.toFixed(0)} ms`);
82+}
83+
84+// Velocity changes the sound (normalize is off, so amplitude scales).
85+const soft = await renderRequest('cymbal2', req([note(0, 60, 0.15)]), env);
86+{
87+ const pa = stats(a.audio).peak;
88+ const ps = stats(soft.audio).peak;
89+ check('velocity changes the sound', ps < pa, `soft ${ps.toExponential(2)} < hard ${pa.toExponential(2)}`);
90+}
91+
92+// A pattern: both strokes come from the cache and land at their onsets.
93+const tr = await renderRequest('cymbal2', req([note(0, 60, 0.4), note(0.5, 60, 0.4), note(1.0, 60, 0.15)]), env);
94+{
95+ const s = stats(tr.audio);
96+ check('track render is finite', s.finite);
97+ const at = (sec: number) => stats(tr.audio.slice(Math.floor(sec * SR), Math.floor((sec + 0.25) * SR))).peak;
98+ check('energy at first onset', at(0.0) > 1e-6);
99+ check('energy at second onset', at(0.5) > 1e-6);
100+ check('energy at third onset', at(1.0) > 1e-6);
101+}
102+
103+console.log(`total time: ${((performance.now() - t0) / 1000).toFixed(1)} s`);
104+if (failures > 0) { console.error(`${failures} check(s) failed`); process.exit(1); }
instrument/types.tsadded+91−0View file
@@ -0,0 +1,91 @@
1+// The track-render instrument contract shared between midiboot and the
2+// instrument web apps. Each simulated-instrument repo carries its own
3+// structurally identical copy of these types (there is no shared package;
4+// TypeScript's structural typing keeps the copies interchangeable), and
5+// midiboot's src/instruments/types.ts is the reference.
6+//
7+// An instrument receives an entire MIDI track per call: every note with
8+// times already in seconds. It owns its rendering strategy end to end,
9+// including any internal caching (per-note, per-stage, or none), progress
10+// reporting, and cancellation. The optional fields at the bottom exist for
11+// the instrument's own web app, which drives the full parameter surface
12+// and live-edited MATLAB sources; midiboot never sets them.
13+
14+export interface NoteSpec {
15+ timeSec: number;
16+ durationSec: number;
17+ pitch: number; // MIDI note number
18+ velocity: number; // 0..1
19+}
20+
21+/** A user-editable instrument parameter, declared by the instrument. */
22+export interface KnobSpec {
23+ key: string;
24+ label: string;
25+ min: number;
26+ max: number;
27+ default: number;
28+ step?: number;
29+ scale?: 'log';
30+}
31+
32+/** One MATLAB source file, as numbl consumes it. */
33+export interface MFile {
34+ name: string;
35+ source: string;
36+}
37+
38+export interface RenderProgress {
39+ /** 0..1 through the whole request. */
40+ frac: number;
41+ /** Short human-readable stage label, e.g. "note 3/8: solve". */
42+ stage: string;
43+}
44+
45+export interface TrackRenderRequest {
46+ /** Every note of the track, sorted by timeSec. */
47+ notes: NoteSpec[];
48+ sampleRate: number;
49+ /** Song length incl. tail, as a minimum: the returned buffer may be
50+ * longer, so a ring that outlives the last note (a cymbal crash, a
51+ * long bass tail) is kept in full; the host mixes buffers of unequal
52+ * lengths. */
53+ lengthSamples: number;
54+ /** Fast preview for playback, full quality for export. */
55+ quality: 'preview' | 'full';
56+ /** Values for the instrument's knobs; missing keys mean the default. */
57+ settings?: Record<string, number>;
58+ /** Cooperative cancellation; checked between stages and notes. */
59+ signal?: AbortSignal;
60+ onProgress?: (p: RenderProgress) => void;
61+
62+ // --- fields used only by the instrument's own web app ---
63+
64+ /** Full model parameter overrides, applied after the defaults (and after
65+ * any note mapping). midiboot never sets this. */
66+ params?: Record<string, number | boolean>;
67+ /** Live-edited MATLAB sources replacing the shipped ones. */
68+ mfiles?: MFile[];
69+ /** Also compute the instrument-specific diagnostics payload. */
70+ wantDiagnostics?: boolean;
71+}
72+
73+export interface TrackRenderResult {
74+ audio: Float32Array;
75+ /** Instrument-specific diagnostics; present only when requested. */
76+ diag?: unknown;
77+}
78+
79+export interface LiveVoice {
80+ noteOff(): void; // begin release
81+ cut(): void; // stop immediately
82+}
83+
84+export interface Instrument {
85+ id: string;
86+ name: string;
87+ /** Curated editable parameters; tracks store their values in settings. */
88+ knobs?: KnobSpec[];
89+ renderTrack(req: TrackRenderRequest): Promise<TrackRenderResult>;
90+ noteOn?(ctx: BaseAudioContext, dest: AudioNode, pitch: number, velocity: number, whenSec: number): LiveVoice;
91+}
instrument/worker.tsadded+106−0View file
@@ -0,0 +1,106 @@
1+// The render worker: numbl (and, when available, WebGPU) runs here so a
2+// minutes-long time integration never blocks the UI. Requests are
3+// processed one at a time; the core yields to the macrotask queue between
4+// stages and awaits the GPU between chunks, which is when queued cancel
5+// messages are delivered, so cancellation is cooperative and the caches
6+// survive it.
7+import { RenderCancelled, renderRequest, type Cymbal2Variant } from './core.ts';
8+import { DEFAULT_MFILES } from './mfiles.ts';
9+import type { MFile, NoteSpec } from './types.ts';
10+
11+export interface WorkerRenderMsg {
12+ kind: 'render';
13+ id: number;
14+ variant: Cymbal2Variant;
15+ notes: NoteSpec[];
16+ sampleRate: number;
17+ lengthSamples: number;
18+ quality: 'preview' | 'full';
19+ settings?: Record<string, number>;
20+ params?: Record<string, number | boolean>;
21+ mfiles?: MFile[];
22+ wantDiagnostics?: boolean;
23+}
24+
25+export interface WorkerCancelMsg {
26+ kind: 'cancel';
27+ id: number;
28+}
29+
30+export type WorkerResponse =
31+ | { kind: 'progress'; id: number; frac: number; stage: string }
32+ | { kind: 'result'; id: number; audio: Float32Array; diag?: unknown }
33+ | { kind: 'cancelled'; id: number }
34+ | { kind: 'error'; id: number; message: string; file?: string; line?: number };
35+
36+const cancelled = new Set<number>();
37+const queue: WorkerRenderMsg[] = [];
38+let running = false;
39+
40+const post = (m: WorkerResponse, transfer?: Transferable[]) =>
41+ (self as unknown as Worker).postMessage(m, transfer ?? []);
42+
43+async function run(msg: WorkerRenderMsg): Promise<void> {
44+ try {
45+ const result = await renderRequest(msg.variant, {
46+ notes: msg.notes,
47+ sampleRate: msg.sampleRate,
48+ lengthSamples: msg.lengthSamples,
49+ quality: msg.quality,
50+ settings: msg.settings,
51+ params: msg.params,
52+ wantDiagnostics: msg.wantDiagnostics,
53+ }, {
54+ mfiles: msg.mfiles ?? DEFAULT_MFILES,
55+ shipped: DEFAULT_MFILES,
56+ post: (frac, stage) => post({ kind: 'progress', id: msg.id, frac, stage }),
57+ checkCancel: () => { if (cancelled.has(msg.id)) throw new RenderCancelled(); },
58+ yieldToHost: () => new Promise((res) => setTimeout(res, 0)),
59+ });
60+ const transfer: Transferable[] = [result.audio.buffer];
61+ if (result.diag) {
62+ for (const v of Object.values(result.diag as unknown as Record<string, unknown>)) {
63+ if (v instanceof Float32Array && v.buffer !== result.audio.buffer) transfer.push(v.buffer);
64+ }
65+ }
66+ post({ kind: 'result', id: msg.id, audio: result.audio, diag: result.diag }, transfer);
67+ } catch (err) {
68+ if (err instanceof RenderCancelled) {
69+ post({ kind: 'cancelled', id: msg.id });
70+ } else {
71+ // Faults in the (user-editable) .m code arrive as errors carrying a
72+ // file and line, which the core's parse check has already attributed.
73+ const e = err as { message?: string; file?: string; line?: number };
74+ post({
75+ kind: 'error',
76+ id: msg.id,
77+ message: e?.message ? String(e.message) : String(err),
78+ file: typeof e?.file === 'string' ? e.file : undefined,
79+ line: typeof e?.line === 'number' ? e.line : undefined,
80+ });
81+ }
82+ } finally {
83+ cancelled.delete(msg.id);
84+ }
85+}
86+
87+async function pump(): Promise<void> {
88+ if (running) return;
89+ running = true;
90+ while (queue.length > 0) {
91+ const msg = queue.shift()!;
92+ if (cancelled.has(msg.id)) { post({ kind: 'cancelled', id: msg.id }); cancelled.delete(msg.id); continue; }
93+ await run(msg);
94+ }
95+ running = false;
96+}
97+
98+self.onmessage = (ev: MessageEvent<WorkerRenderMsg | WorkerCancelMsg>) => {
99+ const msg = ev.data;
100+ if (msg.kind === 'render') {
101+ queue.push(msg);
102+ void pump();
103+ } else if (msg.kind === 'cancel') {
104+ cancelled.add(msg.id);
105+ }
106+};
scripts/build-exports.mjsdeleted+0−79View file
@@ -1,79 +0,0 @@
1-#!/usr/bin/env node
2-// Build exported instrument bundles (format numbl-instrument-v1) for one
3-// instrument repository. Run from the repo root:
4-//
5-// node scripts/build-exports.mjs
6-//
7-// For each export-src/<name>/export.json, the bundle's files are that
8-// directory's .m sources (the drivers, render_note.m required) plus every
9-// matlab/*.m of the repo, and the result is written to exports/<id>.json.
10-// The canonical copy of this script lives in
11-// simulated_instruments/instrument-format; instrument repos carry a copy so
12-// they build standalone.
13-
14-import { mkdirSync, readdirSync, readFileSync, writeFileSync, existsSync } from 'node:fs';
15-import { join } from 'node:path';
16-
17-const ID_RE = /^[a-z0-9][a-z0-9-]*$/;
18-const SRC = 'export-src';
19-const OUT = 'exports';
20-
21-function fail(msg) {
22- console.error(`build-exports: ${msg}`);
23- process.exit(1);
24-}
25-
26-function checkKnobs(id, knobs) {
27- if (!Array.isArray(knobs)) fail(`${id}: 'knobs' must be an array`);
28- for (const k of knobs) {
29- if (typeof k.key !== 'string' || typeof k.label !== 'string') fail(`${id}: knobs need key and label`);
30- if (typeof k.min !== 'number' || typeof k.max !== 'number' || typeof k.default !== 'number' || !(k.min < k.max)) {
31- fail(`${id}: knob '${k.key}' needs min < max and a default`);
32- }
33- if (k.default < k.min || k.default > k.max) fail(`${id}: knob '${k.key}' default outside its range`);
34- if (k.step !== undefined && !(k.step > 0)) fail(`${id}: knob '${k.key}' step must be > 0`);
35- if (k.scale !== undefined && k.scale !== 'log') fail(`${id}: knob '${k.key}' has unknown scale`);
36- if (k.scale === 'log' && k.min <= 0) fail(`${id}: knob '${k.key}' log scale needs min > 0`);
37- }
38-}
39-
40-function mfilesIn(dir) {
41- if (!existsSync(dir)) return [];
42- return readdirSync(dir)
43- .filter((f) => f.endsWith('.m'))
44- .sort()
45- .map((name) => ({ name, content: readFileSync(join(dir, name), 'utf8') }));
46-}
47-
48-if (!existsSync(SRC)) fail(`no ${SRC}/ directory here; run from the instrument repo root`);
49-const shared = mfilesIn('matlab');
50-const dirs = readdirSync(SRC, { withFileTypes: true }).filter((d) => d.isDirectory()).map((d) => d.name).sort();
51-if (dirs.length === 0) fail(`${SRC}/ has no instrument directories`);
52-mkdirSync(OUT, { recursive: true });
53-
54-for (const dir of dirs) {
55- const metaPath = join(SRC, dir, 'export.json');
56- if (!existsSync(metaPath)) fail(`${metaPath} is missing`);
57- const meta = JSON.parse(readFileSync(metaPath, 'utf8'));
58- const id = meta.id;
59- if (typeof id !== 'string' || !ID_RE.test(id)) fail(`${metaPath}: 'id' must match ${ID_RE}`);
60- if (typeof meta.name !== 'string' || meta.name === '') fail(`${id}: missing 'name'`);
61- checkKnobs(id, meta.knobs ?? []);
62- const drivers = mfilesIn(join(SRC, dir));
63- if (!drivers.some((f) => f.name === 'render_note.m')) fail(`${id}: no render_note.m in ${SRC}/${dir}/`);
64- const names = new Set(drivers.map((f) => f.name));
65- const clash = shared.find((f) => names.has(f.name));
66- if (clash) fail(`${id}: driver file '${clash.name}' collides with matlab/${clash.name}`);
67- const bundle = {
68- format: 'numbl-instrument-v1',
69- id,
70- name: meta.name,
71- ...(meta.description ? { description: meta.description } : {}),
72- ...(meta.source ? { source: meta.source } : {}),
73- knobs: meta.knobs ?? [],
74- files: [...drivers, ...shared],
75- };
76- const outPath = join(OUT, `${id}.json`);
77- writeFileSync(outPath, JSON.stringify(bundle, null, 2) + '\n');
78- console.log(`wrote ${outPath} (${bundle.files.length} files, ${bundle.knobs.length} knobs)`);
79-}
webapp/README.mdmodified+1−1View file
@@ -21,7 +21,7 @@ The self-test runs the closed-form and reproducibility checks and the pinned dia
2121
2222 Everywhere. Unlike the sibling `cymbal` app, whose expensive stage (the eigensolves) could be cached and whose synthesis was an oscillator bank, here the sound is a nonlinear time-domain integration of the whole shell, several hundred thousand steps of coupled solves, and every parameter feeds it. The Render preview button runs a cheap pass (no oversampling; the duration and the azimuthal grid are those set, the latter because hard strokes diverge without the full wavenumber ladder), Render full runs the real thing, and either can be cancelled while it runs; the status shows the stage, a progress bar and the elapsed time. This is the same control set the sibling instrument apps carry, and there as here nothing renders on its own. The only piece the worker caches is the cheap linear mode table that the spectrum plot overlays.
2323
24-When the browser exposes WebGPU, the time loop runs on the GPU instead of through the interpreter: the worker still executes the unmodified setup stage (grid, per-wavenumber inverses, mallet, stroke schedule) and output stage through numbl, but hands the operators to a WGSL port of `cymbal2_sim.m` (`src/worker/gpuSim.ts`, eight compute dispatches per time step; the design and its measurements are documented in the [cymbal2-webgpu](https://vault1.magland.org/proof-of-concept/cymbal2-webgpu) exploration). This is 20x or more faster than interpreting the loop, even on an integrated GPU, and the diagnostics panel says which path a render took. The GPU integrates in single precision; the contact diagnostics match the double-precision interpreter to a fraction of a percent, and the trajectory of a crash is chaotic either way. Without WebGPU the interpreter path runs as before, and expect a minute or two for a preview and many minutes for a full render. One interaction with the code editor: the WGSL mirrors the shipped `cymbal2_sim.m`, so editing that one file falls back to the interpreter, where the edit is live (a warning in the diagnostics says so); edits to every other `.m` file flow through numbl on both paths. `npm i --no-save webgpu && npx tsx scripts/gpu-selftest.ts` runs the GPU loop from Node on desktop WebGPU and checks the pinned diagnostics.
24+When the browser exposes WebGPU, the time loop runs on the GPU instead of through the interpreter: the worker still executes the unmodified setup stage (grid, per-wavenumber inverses, mallet, stroke schedule) and output stage through numbl, but hands the operators to a WGSL port of `cymbal2_sim.m` (`../instrument/gpuSim.ts`, eight compute dispatches per time step; the design and its measurements are documented in the [cymbal2-webgpu](https://vault1.magland.org/proof-of-concept/cymbal2-webgpu) exploration). This is 20x or more faster than interpreting the loop, even on an integrated GPU, and the diagnostics panel says which path a render took. The GPU integrates in single precision; the contact diagnostics match the double-precision interpreter to a fraction of a percent, and the trajectory of a crash is chaotic either way. Without WebGPU the interpreter path runs as before, and expect a minute or two for a preview and many minutes for a full render. One interaction with the code editor: the WGSL mirrors the shipped `cymbal2_sim.m`, so editing that one file falls back to the interpreter, where the edit is live (a warning in the diagnostics says so); edits to every other `.m` file flow through numbl on both paths. `npm i --no-save webgpu && npx tsx scripts/gpu-selftest.ts` runs the GPU loop from Node on desktop WebGPU and checks the pinned diagnostics.
2525
2626 ## Presets
2727
webapp/scripts/gpu-selftest.tsmodified+2−2View file
@@ -1,4 +1,4 @@
1-// Runs the webapp's WebGPU time loop (src/worker/gpuSim.ts) from Node against
1+// Runs the WebGPU time loop (../instrument/gpuSim.ts) from Node against
22 // desktop WebGPU (the optional `webgpu` package, prebuilt Google Dawn) on the
33 // validation suite's reduced configuration, and checks the pinned contact
44 // diagnostics and output rms. This exercises exactly the code the browser
@@ -50,7 +50,7 @@ async function main() {
5050 console.log('gpu-selftest: `webgpu` package not installed; skipping (npm i --no-save webgpu)');
5151 return;
5252 }
53- const { requestGpuDevice, runGpuSim } = await import('../src/worker/gpuSim');
53+ const { requestGpuDevice, runGpuSim } = await import('../../instrument/gpuSim.ts');
5454
5555 const files = mfiles();
5656 const P = `
webapp/src/components/DiagnosticsPanel.tsxmodified+2−2View file
@@ -1,6 +1,6 @@
1-import type { RenderResponse } from '../worker/renderWorker';
1+import type { RenderResult } from '../hooks/useRenderer';
22
3-export default function DiagnosticsPanel({ result }: { result: RenderResponse }) {
3+export default function DiagnosticsPanel({ result }: { result: RenderResult }) {
44 const d = result.diag;
55 return (
66 <div className="diag">
webapp/src/hooks/useRenderer.tsmodified+72−66View file
@@ -1,13 +1,24 @@
1-// Hook owning the render worker. Unlike the sibling cymbal app there is no
2-// debounced auto-preview: a cymbal2 render is a full time-domain integration
3-// and takes from under a minute (preview) to minutes (full quality) in the
4-// browser, so renders are explicit. A request in flight can be cancelled,
5-// which tears the worker down and builds a fresh one.
1+// Hook owning the render requests. Renders are explicit: nothing starts
2+// until Render preview or Render full is pressed, so changing a parameter
3+// or a .m source never kicks off work on its own. Rendering goes through
4+// the repo's instrument/ module (the same track-render interface midiboot
5+// uses), with the full parameter surface, the edited .m sources, and the
6+// diagnostics payload passed through the request's app-facing fields.
7+// Cancelling is cooperative: the worker checks between render stages, and
8+// the instrument module tears the worker down only if it never responds.
69
710 import { useCallback, useEffect, useRef, useState } from 'react';
811 import type { Cymbal2Params } from '../params';
912 import type { MFile } from '../mfiles';
10-import type { ErrorMessage, ProgressMessage, RenderResponse } from '../worker/renderWorker';
13+import { cymbal2Instrument } from '../../../instrument/index.ts';
14+import type { Cymbal2ExploreDiag } from '../../../instrument/core.ts';
15+
16+/** The old worker response shape, reassembled from the instrument result
17+ * so the plot components are unchanged. */
18+export type RenderResult = Cymbal2ExploreDiag & {
19+ y: Float32Array;
20+ sampleRate: number;
21+};
1122
1223 export interface RenderError {
1324 message: string;
@@ -16,7 +27,7 @@ export interface RenderError {
1627 }
1728
1829 export interface RenderState {
19- result: RenderResponse | null;
30+ result: RenderResult | null;
2031 error: RenderError | null;
2132 busy: boolean;
2233 stage: string;
@@ -25,9 +36,7 @@ export interface RenderState {
2536 }
2637
2738 export function useRenderer(params: Cymbal2Params, mfiles: MFile[]) {
28- const workerRef = useRef<Worker | null>(null);
29- const nextId = useRef(1);
30- const inFlight = useRef<number | null>(null);
39+ const abortRef = useRef<AbortController | null>(null);
3140 const startedAt = useRef(0);
3241 const [state, setState] = useState<RenderState>({
3342 result: null,
@@ -38,49 +47,7 @@ export function useRenderer(params: Cymbal2Params, mfiles: MFile[]) {
3847 elapsedS: 0,
3948 });
4049
41- const attachHandlers = useCallback((worker: Worker) => {
42- worker.onmessage = (ev: MessageEvent<RenderResponse | ProgressMessage | ErrorMessage>) => {
43- const msg = ev.data;
44- if (msg.type === 'progress') {
45- setState((s) => ({ ...s, stage: msg.stage, progress: msg.frac }));
46- return;
47- }
48- inFlight.current = null;
49- if (msg.type === 'error') {
50- // A failed render leaves nothing playable: the previous result is
51- // dropped, so Play and Download cannot hand back audio that no
52- // longer corresponds to the code or parameters on screen.
53- setState((s) => ({
54- ...s,
55- result: null,
56- error: { message: msg.message, file: msg.file, line: msg.line },
57- busy: false,
58- stage: '',
59- progress: 0,
60- }));
61- return;
62- }
63- setState((s) => ({ ...s, result: msg, error: null, busy: false, stage: '', progress: 0 }));
64- };
65- }, []);
66-
67- const makeWorker = useCallback(() => {
68- const worker = new Worker(new URL('../worker/renderWorker.ts', import.meta.url), {
69- type: 'module',
70- });
71- attachHandlers(worker);
72- workerRef.current = worker;
73- return worker;
74- }, [attachHandlers]);
75-
76- useEffect(() => {
77- const worker = makeWorker();
78- return () => {
79- worker.terminate();
80- workerRef.current = null;
81- inFlight.current = null;
82- };
83- }, [makeWorker]);
50+ useEffect(() => () => abortRef.current?.abort(), []);
8451
8552 // elapsed-time ticker while a render is in flight
8653 useEffect(() => {
@@ -94,27 +61,66 @@ export function useRenderer(params: Cymbal2Params, mfiles: MFile[]) {
9461
9562 const render = useCallback(
9663 (quality: 'preview' | 'full') => {
97- const worker = workerRef.current;
98- if (!worker || inFlight.current !== null) return;
99- const id = nextId.current++;
100- inFlight.current = id;
64+ if (abortRef.current) return; // a request is already in flight
65+ const ctrl = new AbortController();
66+ abortRef.current = ctrl;
10167 startedAt.current = performance.now();
10268 setState((s) => ({ ...s, busy: true, stage: 'starting', progress: 0, elapsedS: 0 }));
103- worker.postMessage({ type: 'render', id, params, mfiles, quality });
69+ const sampleRate = Number(params.sample_rate);
70+ void cymbal2Instrument.renderTrack({
71+ notes: [],
72+ sampleRate,
73+ lengthSamples: 1,
74+ quality,
75+ params: params as unknown as Record<string, number | boolean>,
76+ mfiles,
77+ wantDiagnostics: true,
78+ signal: ctrl.signal,
79+ onProgress: (p) => setState((s) => ({ ...s, stage: p.stage, progress: p.frac })),
80+ }).then((r) => {
81+ if (abortRef.current === ctrl) abortRef.current = null;
82+ if (ctrl.signal.aborted) return;
83+ const diag = r.diag as Cymbal2ExploreDiag;
84+ setState((s) => ({
85+ ...s,
86+ result: { ...diag, y: r.audio, sampleRate },
87+ error: null,
88+ busy: false,
89+ stage: '',
90+ progress: 0,
91+ }));
92+ }).catch((err: unknown) => {
93+ if (abortRef.current === ctrl) abortRef.current = null;
94+ if (err instanceof DOMException && err.name === 'AbortError') {
95+ setState((s) => ({ ...s, busy: false, stage: '', progress: 0 }));
96+ return;
97+ }
98+ // A failed render leaves nothing playable: the previous result is
99+ // dropped, so Play and Download cannot hand back audio that no
100+ // longer corresponds to the code or parameters on screen.
101+ const e = err as { message?: string; file?: string; line?: number };
102+ setState((s) => ({
103+ ...s,
104+ result: null,
105+ error: {
106+ message: e?.message ? String(e.message) : String(err),
107+ file: typeof e?.file === 'string' ? e.file : undefined,
108+ line: typeof e?.line === 'number' ? e.line : undefined,
109+ },
110+ busy: false,
111+ stage: '',
112+ progress: 0,
113+ }));
114+ });
104115 },
105116 [params, mfiles],
106117 );
107118
108- // Cancelling terminates the worker outright: the .m interpreter runs a tight
109- // loop that never yields, so there is nothing to poll a flag from. The
110- // replacement worker starts with empty stage caches.
111119 const cancel = useCallback(() => {
112- const worker = workerRef.current;
113- if (worker) worker.terminate();
114- inFlight.current = null;
115- makeWorker();
120+ abortRef.current?.abort();
121+ abortRef.current = null;
116122 setState((s) => ({ ...s, busy: false, stage: '', progress: 0, elapsedS: 0 }));
117- }, [makeWorker]);
123+ }, []);
118124
119125 const renderPreview = useCallback(() => render('preview'), [render]);
120126 const renderFull = useCallback(() => render('full'), [render]);
webapp/src/worker/renderWorker.tsdeleted+0−628View file
@@ -1,628 +0,0 @@
1-// Web Worker that renders the cymbal2 stroke. The MATLAB implementation in
2-// ../matlab remains the source of truth; what varies is who runs the time
3-// loop. When the browser exposes WebGPU, the worker executes the unmodified
4-// setup stage (grid, per-wavenumber inverses, mallet, stroke schedule)
5-// through numbl, hands the operators to a WGSL port of cymbal2_sim.m
6-// (./gpuSim.ts), and feeds the returned acceleration back through the
7-// unmodified cymbal2_output.m; this is 20x or more faster than interpreting
8-// the loop. Without WebGPU, or when cymbal2_sim.m has been edited in the app
9-// (the WGSL mirrors the shipped step, so edits there must fall back), the
10-// original interpreter path runs instead. Requests carry a quality tier
11-// preview -> oversample 1, everything else as given: the render covers
12-// the full duration, so the decay can be judged, and the
13-// azimuthal grid is kept because hard strokes need the full
14-// wavenumber band to carry the cascade, and coarsening it
15-// destabilizes the integration
16-// full -> the parameters as given
17-// and the UI shows elapsed time, plus integration progress on the GPU path.
18-// A cheap linear-modes stage (cymbal2_modes per wavenumber) is cached on the
19-// geometry so the spectrum plot can overlay the mode table immediately.
20-
21-import { executeCode, RTV, type RuntimeValue, type RuntimeTensor } from 'numbl';
22-import { hashMFiles, DEFAULT_MFILES, EXECUTED_BY_APP, type MFile } from '../mfiles';
23-import { MODES_KEYS, type Cymbal2Params } from '../params';
24-import { requestGpuDevice, runGpuSim, type GpuSetup } from './gpuSim';
25-
26-export interface RenderRequest {
27- type: 'render';
28- id: number;
29- params: Cymbal2Params;
30- quality: 'preview' | 'full';
31- mfiles: MFile[];
32-}
33-
34-export interface RenderResponse {
35- type: 'result';
36- id: number;
37- quality: 'preview' | 'full';
38- y: Float32Array; // final signal, [-1, 1]
39- sampleRate: number;
40- forceT: Float32Array; // contact force history (max-pooled), s
41- forceN: Float32Array; // and its values, N
42- modeF: Float32Array; // linear mode frequencies, Hz
43- modeM: Float32Array; // nodal diameters of each
44- diag: {
45- strokesPlayed: number;
46- meanV: number;
47- contactVel: number;
48- peakForceN: number;
49- contactMs: number;
50- contactEvents: number;
51- peakXiMm: number;
52- peakDispMm: number;
53- peakDispThickness: number;
54- peakAccel: number;
55- strikeRadiusM: number;
56- pickupRadiusM: number;
57- RsM: number;
58- nsteps: number;
59- simRate: number;
60- nDropped: number;
61- durationS: number;
62- renderMs: number;
63- oversample: number;
64- nTheta: number;
65- backend: 'webgpu' | 'cpu';
66- warnings: string[];
67- };
68-}
69-
70-export interface ProgressMessage {
71- type: 'progress';
72- id: number;
73- stage: 'parsing' | 'modes' | 'integrate' | 'output';
74- frac: number;
75-}
76-
77-export interface ErrorMessage {
78- type: 'error';
79- id: number;
80- message: string;
81- file?: string;
82- line?: number;
83-}
84-
85-interface ModesEntry {
86- modeF: Float32Array;
87- modeM: Float32Array;
88-}
89-
90-const modesCache = new Map<string, ModesEntry>();
91-const MODES_CACHE_MAX = 4;
92-
93-function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
94- const v = cache.get(key);
95- if (v !== undefined) {
96- cache.delete(key);
97- cache.set(key, v);
98- }
99- return v;
100-}
101-
102-function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
103- cache.delete(key);
104- cache.set(key, value);
105- while (cache.size > max) {
106- const oldest = cache.keys().next().value as string;
107- cache.delete(oldest);
108- }
109-}
110-
111-function keyFor(params: Cymbal2Params, keys: readonly (keyof Cymbal2Params)[]): string {
112- return keys.map((k) => String(params[k])).join('|');
113-}
114-
115-// MATLAB source lines assigning every parameter onto a struct P.
116-function pLiteral(params: Cymbal2Params): string {
117- const lines: string[] = [];
118- for (const [k, v] of Object.entries(params)) {
119- if (k === 'bit_depth') continue;
120- if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
121- else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
122- }
123- return lines.join('\n');
124-}
125-
126-// An error in the MATLAB source, carrying the file and line where numbl put
127-// it. Reported to the UI, which stops the render and shows nothing playable.
128-class CodeError extends Error {
129- file?: string;
130- line?: number;
131- constructor(message: string, file?: string, line?: number) {
132- super(message);
133- this.name = 'CodeError';
134- this.file = file;
135- this.line = line;
136- }
137-}
138-
139-const SKIP_WARNING = /skipping (\S+\.m) \((?:syntax error at line (\d+)|parse error)\)/;
140-
141-// numbl does not abort on a file it cannot parse: it drops that file with a
142-// warning and carries on with the rest, so a render could otherwise finish
143-// (or fail much later with a confusing 'undefined function') while the real
144-// fault is a typo in a file the stage never reached. Every numbl call runs
145-// with console.warn captured, and a dropped file becomes a CodeError. Only
146-// the files the app executes count: the CLI entry points (cymbal2_demo.m and
147-// the rest) are in the editor but play no part in a render, so a typo there
148-// must not stop one.
149-function skipWarningToError(warnings: string[]): CodeError | null {
150- for (const w of warnings) {
151- const hit = SKIP_WARNING.exec(w);
152- if (!hit) continue;
153- const file = hit[1];
154- if (!EXECUTED_BY_APP.has(file)) continue;
155- const line = hit[2] ? Number(hit[2]) : undefined;
156- return new CodeError(
157- line !== undefined
158- ? 'syntax error; the render was stopped'
159- : 'the file cannot be parsed; the render was stopped',
160- file,
161- line,
162- );
163- }
164- return null;
165-}
166-
167-function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
168- const warnings: string[] = [];
169- const prevWarn = console.warn;
170- console.warn = (...args: unknown[]) => {
171- warnings.push(args.map(String).join(' '));
172- prevWarn(...(args as []));
173- };
174- try {
175- const r = executeCode(
176- source,
177- { initialVariableValues: inputs, displayResults: false },
178- mfiles,
179- 'stage.m',
180- );
181- const skipped = skipWarningToError(warnings);
182- if (skipped) throw skipped;
183- return r;
184- } catch (err) {
185- // A file numbl dropped is the better explanation than whatever the run
186- // failed on afterwards (typically a call to the missing function).
187- const skipped = skipWarningToError(warnings);
188- if (skipped) throw skipped;
189- throw err;
190- } finally {
191- console.warn = prevWarn;
192- }
193-}
194-
195-// Parse every .m file before any expensive stage runs, so a syntax error is
196-// reported immediately instead of after a stage (or, if the modes table is
197-// cached, after the whole integration). Executing a trivial script with the
198-// files attached is enough: numbl pre-parses the whole workspace.
199-function checkMFiles(mfiles: MFile[]): void {
200- runStage(mfiles, '1;');
201-}
202-
203-function asScalar(v: RuntimeValue): number {
204- if (typeof v === 'number') return v;
205- const t = v as RuntimeTensor;
206- if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
207- if (typeof v === 'boolean') return v ? 1 : 0;
208- throw new Error('expected scalar result from numbl stage');
209-}
210-
211-function toF32(v: RuntimeValue): Float32Array {
212- const t = v as RuntimeTensor;
213- if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
214- const out = new Float32Array(t.data.length);
215- for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
216- return out;
217-}
218-
219-// Max-pool the contact force history to about 2000 points for the plot.
220-function poolForce(fhist: Float32Array, simRate: number): { tpool: Float32Array; fpool: Float32Array } {
221- const nsim = fhist.length;
222- const np0 = Math.min(2000, nsim);
223- const blk = Math.ceil(nsim / np0);
224- const np = Math.ceil(nsim / blk);
225- const fpool = new Float32Array(np);
226- const tpool = new Float32Array(np);
227- for (let b = 0; b < np; b++) {
228- let m = 0;
229- const end = Math.min(nsim, (b + 1) * blk);
230- for (let i = b * blk; i < end; i++) m = Math.max(m, fhist[i]);
231- fpool[b] = m;
232- tpool[b] = (b * blk + blk / 2) / simRate;
233- }
234- return { tpool, fpool };
235-}
236-
237-self.onmessage = (ev: MessageEvent<RenderRequest>) => {
238- const msg = ev.data;
239- if (msg.type !== 'render') return;
240-
241- const warnings: string[] = [];
242- const origWarn = console.warn;
243- console.warn = (...args: unknown[]) => {
244- warnings.push(args.map(String).join(' '));
245- };
246-
247- render(msg)
248- .catch((err) => {
249- // runStage turns a dropped file into a CodeError itself; this catches
250- // the same warning from anything that reached numbl another way.
251- const e = (skipWarningToError(warnings) ?? err) as {
252- message?: string;
253- file?: string;
254- line?: number;
255- };
256- const m: ErrorMessage = {
257- type: 'error',
258- id: msg.id,
259- message: e?.message ? String(e.message) : String(err),
260- file: typeof e?.file === 'string' ? e.file : undefined,
261- line: typeof e?.line === 'number' ? e.line : undefined,
262- };
263- (self as unknown as Worker).postMessage(m);
264- })
265- .finally(() => {
266- console.warn = origWarn;
267- });
268-};
269-
270-// Everything the response assembly needs, produced by either backend.
271-interface SimOutcome {
272- y32: Float32Array;
273- peakAccel: number;
274- nDropped: number;
275- tpool: Float32Array;
276- fpool: Float32Array;
277- backend: 'webgpu' | 'cpu';
278- d: {
279- strokesPlayed: number;
280- meanV: number;
281- contactVel: number;
282- peakForceN: number;
283- contactMs: number;
284- contactEvents: number;
285- peakXiMm: number;
286- peakDispM: number;
287- peakDispThickness: number;
288- nsteps: number;
289- simRate: number;
290- strikeRadiusM: number;
291- pickupRadiusM: number;
292- RsM: number;
293- };
294-}
295-
296-// The stroke-schedule block from cymbal2_sim.m, verbatim, so the GPU render
297-// integrates exactly the schedule the interpreter would.
298-const SCHEDULE_SRC = `
299-K = max(1, round(P.strike_count));
300-tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);
301-vstr = zeros(1, K);
302-state = P.seed * 2654435761;
303-for i = 1:K
304- state = mod(state * 1664525 + 1013904223, 4294967296);
305- x = state / 4294967296;
306- vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));
307-end
308-gap = P.strike_gap_m;
309-if gap <= 0
310- gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;
311-end
312-travel = gap ./ vstr + vstr / (2 * P.strike_accel);
313-tlaunch = tstr - travel;
314-off = -min(0, min(tlaunch));
315-tlaunch = tlaunch + off;
316-duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
317-nsteps = floor(duration * P.sample_rate * P.oversample);
318-`;
319-
320-async function renderOnGpu(
321- mfiles: MFile[],
322- params: Cymbal2Params,
323- onProgress: (frac: number) => void,
324-): Promise<SimOutcome> {
325- // A WebGPU device first: if there is none, fail before spending seconds on
326- // the setup stage.
327- const device = await requestGpuDevice();
328- try {
329- const r = runStage(mfiles, `
330-${pLiteral(params)}
331-S = cymbal2_system(P);
332-MB = cymbal2_mallet(P);
333-G = S.G;
334-o_D1 = G.D1; o_D2 = G.D2;
335-o_F = G.F; o_Gi = G.Gi;
336-o_Rd = G.Rd; o_R2d = G.R2d;
337-o_AiB = S.AiB; o_AiC = S.AiC; o_Ai = S.Ai; o_Fi = S.Fi;
338-o_ASiB = MB.ASiB; o_ASiC = MB.ASiC;
339-o_vM = MB.vM; o_vR = MB.vR; o_vK = MB.vK; o_tipcol = MB.tipcol;
340-${SCHEDULE_SRC}
341-o_tlaunch = tlaunch; o_vstr = vstr;
342-c_nrt = G.nrt; c_Nr = G.Nr; c_nth = G.nth; c_hr = G.hr; c_hth = G.hth;
343-c_ht = S.ht; c_Eh = S.Eh; c_Rs = S.Rs; c_rhoh = S.rhoh;
344-c_fsc = S.forceScale; c_nlr0 = S.eomRow0; c_nDropped = S.nDropped;
345-c_ls = S.ls; c_js = S.js; c_lp = S.lp; c_jp = S.jp;
346-c_strikeR = S.strikeRadiusM; c_pickupR = S.pickupRadiusM;
347-c_tipi = MB.tipi; c_nfree = MB.nfree; c_ps = MB.ps;
348-c_holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / S.ht));
349-c_missDepth = gap + 0.06;
350-c_gap = gap; c_K = K; c_nsteps = nsteps;
351-c_simRate = P.sample_rate * P.oversample;
352-`);
353- const v = r.variableValues;
354- const s = (n: string) => asScalar(v[n]);
355- const t = (n: string) => toF32(v[n]);
356- const setup: GpuSetup = {
357- ops: {
358- D1: t('o_D1'), D2: t('o_D2'),
359- F: t('o_F'), Gi: t('o_Gi'),
360- Rd: t('o_Rd'), R2d: t('o_R2d'),
361- AiB: t('o_AiB'), AiC: t('o_AiC'), Ai: t('o_Ai'), Fi: t('o_Fi'),
362- ASiB: t('o_ASiB'), ASiC: t('o_ASiC'),
363- vM: t('o_vM'), vR: t('o_vR'), vK: t('o_vK'), tipcol: t('o_tipcol'),
364- },
365- tlaunch: Array.from((v.o_tlaunch as RuntimeTensor).data),
366- vstr: Array.from((v.o_vstr as RuntimeTensor).data),
367- nrt: s('c_nrt'), Nr: s('c_Nr'), nth: s('c_nth'),
368- hr: s('c_hr'), hth: s('c_hth'),
369- ht: s('c_ht'), Eh: s('c_Eh'), Rs: s('c_Rs'), rhoh: s('c_rhoh'),
370- fsc: s('c_fsc'), nlr0: s('c_nlr0'),
371- ls: s('c_ls'), js: s('c_js'), lp: s('c_lp'), jp: s('c_jp'),
372- tipi: s('c_tipi'), nfree: s('c_nfree'), ps: s('c_ps'),
373- KM: params.head_stiffness, alM: params.head_alpha,
374- aH: params.strike_accel, aHr: params.retract_accel,
375- holdSteps: s('c_holdSteps'), missDepth: s('c_missDepth'),
376- gap: s('c_gap'), K: s('c_K'), nsteps: s('c_nsteps'),
377- };
378- const simRate = s('c_simRate');
379-
380- const sim = await runGpuSim(device, setup, (done, n) => onProgress(done / n));
381-
382- // The output stage through the unmodified cymbal2_output.m.
383- const accT = RTV.tensor(Float64Array.from(sim.acc), [1, sim.acc.length]);
384- const ro = runStage(mfiles, `
385-${pLiteral(params)}
386-[y, peakAccel] = cymbal2_output(P, acc);
387-`, { acc: accT });
388-
389- // Contact bookkeeping from the force history, as cymbal2_sim does it.
390- let peakF = 0, contactSteps = 0, events = 0, prev = 0;
391- for (let i = 0; i < sim.fhist.length; i++) {
392- const f = sim.fhist[i];
393- if (f > peakF) peakF = f;
394- if (f > 0) { contactSteps++; if (prev === 0) events++; }
395- prev = f;
396- }
397- const { tpool, fpool } = poolForce(sim.fhist, simRate);
398- const meanV = setup.vstr.reduce((a, b) => a + b, 0) / setup.vstr.length;
399-
400- return {
401- y32: toF32(ro.variableValues.y),
402- peakAccel: asScalar(ro.variableValues.peakAccel),
403- nDropped: s('c_nDropped'),
404- tpool,
405- fpool,
406- backend: 'webgpu',
407- d: {
408- strokesPlayed: sim.strokesPlayed,
409- meanV,
410- contactVel: sim.contactVel[0] ?? 0,
411- peakForceN: peakF,
412- contactMs: (1000 * contactSteps) / simRate,
413- contactEvents: events,
414- peakXiMm: sim.peakXi * 1000,
415- peakDispM: sim.peakDisp,
416- peakDispThickness: sim.peakDisp / params.thickness_m,
417- nsteps: setup.nsteps,
418- simRate,
419- strikeRadiusM: s('c_strikeR'),
420- pickupRadiusM: s('c_pickupR'),
421- RsM: setup.Rs,
422- },
423- };
424- } finally {
425- device.destroy();
426- }
427-}
428-
429-function renderOnCpu(mfiles: MFile[], params: Cymbal2Params): SimOutcome {
430- const r = runStage(mfiles, `
431-${pLiteral(params)}
432-S = cymbal2_system(P);
433-MB = cymbal2_mallet(P);
434-[acc, dinfo, fhist] = cymbal2_sim(P, S, MB);
435-[y, peakAccel] = cymbal2_output(P, acc);
436-% max-pool the contact force history to about 2000 points for the plot
437-nsim = length(fhist);
438-np = min(2000, nsim);
439-blk = ceil(nsim / np);
440-np = ceil(nsim / blk);
441-fpool = max(reshape([fhist, zeros(1, np * blk - nsim)], blk, np), [], 1);
442-tpool = ((0:np - 1) * blk + blk / 2) / dinfo.simRate;
443-d_strokes = dinfo.strokesPlayed;
444-d_meanV = dinfo.meanV;
445-d_cvel = dinfo.contactVel(1);
446-d_peakF = dinfo.peakForceN;
447-d_contactMs = dinfo.contactMs;
448-d_events = dinfo.contactEvents;
449-d_xi = dinfo.peakXiMm;
450-d_disp = dinfo.peakDispM;
451-d_dispH = dinfo.peakDispThickness;
452-d_nsteps = dinfo.nsteps;
453-d_simRate = dinfo.simRate;
454-d_dropped = S.nDropped;
455-d_strikeR = S.strikeRadiusM;
456-d_pickupR = S.pickupRadiusM;
457-d_Rs = S.Rs;
458-`);
459- const v = r.variableValues;
460- return {
461- y32: toF32(v.y),
462- peakAccel: asScalar(v.peakAccel),
463- nDropped: asScalar(v.d_dropped),
464- tpool: toF32(v.tpool),
465- fpool: toF32(v.fpool),
466- backend: 'cpu',
467- d: {
468- strokesPlayed: asScalar(v.d_strokes),
469- meanV: asScalar(v.d_meanV),
470- contactVel: asScalar(v.d_cvel),
471- peakForceN: asScalar(v.d_peakF),
472- contactMs: asScalar(v.d_contactMs),
473- contactEvents: asScalar(v.d_events),
474- peakXiMm: asScalar(v.d_xi),
475- peakDispM: asScalar(v.d_disp),
476- peakDispThickness: asScalar(v.d_dispH),
477- nsteps: asScalar(v.d_nsteps),
478- simRate: asScalar(v.d_simRate),
479- strikeRadiusM: asScalar(v.d_strikeR),
480- pickupRadiusM: asScalar(v.d_pickupR),
481- RsM: asScalar(v.d_Rs),
482- },
483- };
484-}
485-
486-async function render(msg: RenderRequest): Promise<void> {
487- const { id, mfiles, quality } = msg;
488- const params = { ...msg.params };
489- if (quality === 'preview') {
490- params.oversample = 1;
491- }
492-
493- const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
494-
495- // Cheap, and it must come first: a file that does not parse stops the
496- // render here rather than partway through a stage.
497- post({ type: 'progress', id, stage: 'parsing', frac: 0.01 });
498- checkMFiles(mfiles);
499-
500- const codeKey = hashMFiles(mfiles);
501-
502- // ---- linear modes, cheap, cached on geometry ----
503- const modesKey = codeKey + '|' + keyFor(params, MODES_KEYS);
504- let modesEntry = lruGet(modesCache, modesKey);
505- if (!modesEntry) {
506- post({ type: 'progress', id, stage: 'modes', frac: 0.05 });
507- const r = runStage(mfiles, `
508-${pLiteral(params)}
509-mF = [];
510-mM = [];
511-for mw = 0:min(10, floor(P.n_theta / 2))
512- f = cymbal2_modes(P, mw);
513- f = f(f > 0.5 & f < 4000);
514- nk = min(4, length(f));
515- mF = [mF; f(1:nk)];
516- mM = [mM; mw * ones(nk, 1)];
517-end
518-`);
519- modesEntry = {
520- modeF: toF32(r.variableValues.mF),
521- modeM: toF32(r.variableValues.mM),
522- };
523- lruSet(modesCache, modesKey, modesEntry, MODES_CACHE_MAX);
524- }
525-
526- // ---- the time integration ----
527- post({ type: 'progress', id, stage: 'integrate', frac: 0.1 });
528- const t0 = performance.now();
529- const warnings: string[] = [];
530-
531- // The WGSL mirrors the shipped cymbal2_sim.m, so an edited step falls back
532- // to the interpreter, where the edit is live. Edits to the setup and output
533- // stages flow through numbl on both paths.
534- const simSource = mfiles.find((f) => f.name === 'cymbal2_sim.m')?.source;
535- const simShipped = DEFAULT_MFILES.find((f) => f.name === 'cymbal2_sim.m')?.source;
536- const simEdited = simSource !== simShipped;
537-
538- let out: SimOutcome | null = null;
539- if (!simEdited) {
540- try {
541- out = await renderOnGpu(mfiles, params, (frac) =>
542- post({ type: 'progress', id, stage: 'integrate', frac: 0.1 + 0.85 * frac }),
543- );
544- } catch (err) {
545- const detail = err instanceof Error ? err.message : String(err);
546- // An error in the MATLAB source (it carries a file, or names the
547- // divergence guard) is a fault to report, not a reason to retry the
548- // same code on the CPU. Only a genuine GPU problem falls back.
549- const fromCode =
550- err instanceof CodeError ||
551- typeof (err as { file?: unknown }).file === 'string' ||
552- (err as { name?: string }).name === 'RuntimeError';
553- if (fromCode || /diverged/.test(detail)) throw err;
554- // console.warn is intercepted during a render; use info so the reason
555- // is visible in the devtools console.
556- console.info(`cymbal2: WebGPU unavailable (${detail}); rendering on the CPU`);
557- }
558- }
559- if (!out) {
560- out = renderOnCpu(mfiles, params);
561- if (simEdited) {
562- // Not a fault: just say why the render took the slow path.
563- warnings.push('cymbal2_sim.m differs from the shipped source, so the render used the interpreter (the GPU path mirrors the shipped step)');
564- }
565- }
566- const renderMs = performance.now() - t0;
567-
568- post({ type: 'progress', id, stage: 'output', frac: 0.97 });
569- const { y32, peakAccel, nDropped } = out;
570-
571- if (nDropped > 0) {
572- warnings.push(
573- `${nDropped} azimuthal columns were truncated: the radial grid cannot resolve their inner boundary layer; raise n_radial`,
574- );
575- }
576- if (!params.normalize) {
577- let peak = 0;
578- for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
579- if (peak > 1.0) {
580- warnings.push(
581- `clipping: peak acceleration ${peakAccel.toExponential(2)} m/s² exceeds full_scale_accel`,
582- );
583- }
584- }
585-
586- const resp: RenderResponse = {
587- type: 'result',
588- id,
589- quality,
590- y: y32,
591- sampleRate: params.sample_rate,
592- forceT: out.tpool,
593- forceN: out.fpool,
594- modeF: modesEntry.modeF.slice(),
595- modeM: modesEntry.modeM.slice(),
596- diag: {
597- strokesPlayed: out.d.strokesPlayed,
598- meanV: out.d.meanV,
599- contactVel: out.d.contactVel,
600- peakForceN: out.d.peakForceN,
601- contactMs: out.d.contactMs,
602- contactEvents: out.d.contactEvents,
603- peakXiMm: out.d.peakXiMm,
604- peakDispMm: out.d.peakDispM * 1000,
605- peakDispThickness: out.d.peakDispThickness,
606- peakAccel,
607- strikeRadiusM: out.d.strikeRadiusM,
608- pickupRadiusM: out.d.pickupRadiusM,
609- RsM: out.d.RsM,
610- nsteps: out.d.nsteps,
611- simRate: out.d.simRate,
612- nDropped,
613- durationS: y32.length / params.sample_rate,
614- renderMs,
615- oversample: params.oversample,
616- nTheta: params.n_theta,
617- backend: out.backend,
618- warnings,
619- },
620- };
621- (self as unknown as Worker).postMessage(resp, [
622- resp.y.buffer,
623- resp.forceT.buffer,
624- resp.forceN.buffer,
625- resp.modeF.buffer,
626- resp.modeM.buffer,
627- ]);
628-}
webapp/tsconfig.jsonmodified+1−0View file
@@ -6,6 +6,7 @@
66 "module": "ESNext",
77 "skipLibCheck": true,
88 "moduleResolution": "bundler",
9+ "allowImportingTsExtensions": true,
910 "resolveJsonModule": true,
1011 "isolatedModules": true,
1112 "noEmit": true,