Add instrument/: the triangle 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). One FEM eigensolve per pitch-scaled frame is shared
across velocities and knobs, and finished renders persist in IndexedDB
with LRU eviction and a 30-day age-out. Replaces the exported bundle;
export-src/ and exports/ are removed.
22 changed files+3622−714
.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 @@ There is no microphone model: no directivity, no capsule response, no gain stagi
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; `npm test` there runs the model in node) |
1213 | `webapp/` | Interactive browser app that executes the MATLAB code client-side |
1314
1415 The MATLAB code runs from the command line with `cd matlab && numbl run triangle_demo.m`, or `triangle_validate.m` to run the checks. The web app is a Vite project; see `webapp/README.md`.
export-src/triangle/export.jsondeleted+0−29View file
@@ -1,29 +0,0 @@
1-{
2- "id": "triangle",
3- "name": "Triangle",
4- "description": "Struck triangle; pitch scales the rod frame around MIDI 84 (72..96), velocity onto beater speed.",
5- "source": "simulated_instruments/triangle",
6- "knobs": [
7- {
8- "key": "strike_pos",
9- "label": "strike position",
10- "min": 0.02,
11- "max": 0.98,
12- "default": 0.17
13- },
14- {
15- "key": "strike_angle_deg",
16- "label": "strike angle (deg)",
17- "min": 0,
18- "max": 90,
19- "default": 20
20- },
21- {
22- "key": "tail_s",
23- "label": "ring-out (s)",
24- "min": 1,
25- "max": 8,
26- "default": 5
27- }
28- ]
29-}
export-src/triangle/render_note.mdeleted+0−28View file
@@ -1,28 +0,0 @@
1-function y = render_note(note)
2-% MIDI pitch scales the rod frame: side length, gap, and rod diameter by
3-% s = 2^-((pitch - 84)/12), which shifts the bending-mode spectrum by the
4-% semitone ratio exactly (rod modes go as diameter / length^2). The
5-% beater stays fixed. Pitch 84 is the default 6-inch triangle, playable
6-% 72..96; the spectrum is inharmonic, so the note number is a handle on
7-% register, not a tuned fundamental. Velocity drives the beater speed.
8-% Fixed: a single strike, a fixed seed, normalize off so velocity
9-% survives in pascals.
10-P = triangle_defaults();
11-P.sample_rate = note.sample_rate;
12-P.strike_count = 1;
13-P.seed = 1;
14-P.pre_delay_ms = 0;
15-P.normalize = 0;
16-if note.preview ~= 0
17- P.oversample = 1;
18-end
19-s = 2 ^ (-(max(72, min(96, note.pitch)) - 84) / 12);
20-P.side_length_m = P.side_length_m * s;
21-P.gap_m = P.gap_m * s;
22-P.rod_diameter_m = P.rod_diameter_m * s;
23-P.strike_velocity = 0.1 + (6 - 0.1) * note.velocity;
24-P.strike_pos = note.params.strike_pos;
25-P.strike_angle_deg = note.params.strike_angle_deg;
26-P.tail_s = note.params.tail_s;
27-y = triangle_render(P);
28-end
exports/triangle.jsondeleted+0−80View file
@@ -1,80 +0,0 @@
1-{
2- "format": "numbl-instrument-v1",
3- "id": "triangle",
4- "name": "Triangle",
5- "description": "Struck triangle; pitch scales the rod frame around MIDI 84 (72..96), velocity onto beater speed.",
6- "source": "simulated_instruments/triangle",
7- "knobs": [
8- {
9- "key": "strike_pos",
10- "label": "strike position",
11- "min": 0.02,
12- "max": 0.98,
13- "default": 0.17
14- },
15- {
16- "key": "strike_angle_deg",
17- "label": "strike angle (deg)",
18- "min": 0,
19- "max": 90,
20- "default": 20
21- },
22- {
23- "key": "tail_s",
24- "label": "ring-out (s)",
25- "min": 1,
26- "max": 8,
27- "default": 5
28- }
29- ],
30- "files": [
31- {
32- "name": "render_note.m",
33- "content": "function y = render_note(note)\n% MIDI pitch scales the rod frame: side length, gap, and rod diameter by\n% s = 2^-((pitch - 84)/12), which shifts the bending-mode spectrum by the\n% semitone ratio exactly (rod modes go as diameter / length^2). The\n% beater stays fixed. Pitch 84 is the default 6-inch triangle, playable\n% 72..96; the spectrum is inharmonic, so the note number is a handle on\n% register, not a tuned fundamental. Velocity drives the beater speed.\n% Fixed: a single strike, a fixed seed, normalize off so velocity\n% survives in pascals.\nP = triangle_defaults();\nP.sample_rate = note.sample_rate;\nP.strike_count = 1;\nP.seed = 1;\nP.pre_delay_ms = 0;\nP.normalize = 0;\nif note.preview ~= 0\n P.oversample = 1;\nend\ns = 2 ^ (-(max(72, min(96, note.pitch)) - 84) / 12);\nP.side_length_m = P.side_length_m * s;\nP.gap_m = P.gap_m * s;\nP.rod_diameter_m = P.rod_diameter_m * s;\nP.strike_velocity = 0.1 + (6 - 0.1) * note.velocity;\nP.strike_pos = note.params.strike_pos;\nP.strike_angle_deg = note.params.strike_angle_deg;\nP.tail_s = note.params.tail_s;\ny = triangle_render(P);\nend\n"
34- },
35- {
36- "name": "triangle_basis.m",
37- "content": "function basis = triangle_basis(P)\n% TRIANGLE_BASIS Mode set of the open triangle: FEM assembly and eigensolve.\n%\n% The planar geometry splits the frame into two exactly decoupled families\n% (see triangle_fem.m): in-plane modes (bending + stretching) and\n% out-of-plane modes (bending + torsion, converted into each other at the\n% two closed corners). Each family gives a generalized eigenproblem\n% K phi = om^2 M phi, reduced to symmetric standard form through the\n% Cholesky factor of the consistent mass matrix and diagonalized by svd\n% rather than eig: the reduced matrix is symmetric positive semidefinite,\n% so its singular vectors are its eigenvectors, and svd keeps them\n% orthonormal even for the near-degenerate pairs this geometry produces\n% (eig's nonsymmetric path returns complex, non-orthogonal vectors there).\n% The eigenvectors then come back mass-normalized: phi' M phi = I. The\n% three rigid-body modes of each family sit at zero frequency and are\n% dropped along with anything below min_mode_hz; modes above max_mode_hz\n% are truncated (the broadband part they carry is reinstated by the direct\n% term in triangle_radiate.m).\n%\n% For radiation and excitation the basis stores, per retained mode, the\n% displacement at every node (NUx, NUy, NUz; out-of-plane modes have only\n% NUz, in-plane modes only NUx, NUy) and the transverse displacement at\n% every element center (EUpx, EUpy, EUpz), with the axial component removed:\n% axial motion of a thin rod displaces almost no air. Element centers,\n% tangents, and lengths are kept alongside so the radiation stage can treat\n% each element as a point dipole.\n%\n% Damping: zeta = 0.5 (loss_factor + loss_slope om) plus radiation damping\n% from the incoherent sum of the element dipole powers,\n% zeta_rad = rho0 om^3 sum_e (Vdip_e^2 |u_perp,e|^2) / (24 pi c^3),\n% Vdip_e = 2 pi a^2 l_e (the factor 2 is the added-mass coefficient of a\n% transversally moving cylinder). The incoherent sum ignores the\n% cancellation between anti-phase sections, which matters only below the\n% coincidence frequency (about 1.9 kHz for an 8 mm steel rod); above it the\n% flexural waves are supersonic and the estimate is fair.\n\ngeo = triangle_geometry(P);\nnn = geo.nn;\nnel = nn - 1;\nE = P.youngs_modulus;\nrho = P.density;\nnu = P.poisson_ratio;\nd = P.rod_diameter_m;\nc = P.speed_of_sound;\nrho0 = P.air_density;\nnodes = geo.nodes;\n\n% ---- element geometry ----\nex = zeros(nel, 1);\ney = zeros(nel, 1);\ntx = zeros(nel, 1);\nty = zeros(nel, 1);\nel = zeros(nel, 1);\nfor e = 1:nel\n dx = nodes(1, e + 1) - nodes(1, e);\n dy = nodes(2, e + 1) - nodes(2, e);\n el(e) = sqrt(dx^2 + dy^2);\n tx(e) = dx / el(e);\n ty(e) = dy / el(e);\n ex(e) = (nodes(1, e) + nodes(1, e + 1)) / 2;\n ey(e) = (nodes(2, e) + nodes(2, e + 1)) / 2;\nend\n\n% ---- eigensolve per family, collect retained modes ----\nfs = [];\nfams = [];\nNUx = [];\nNUy = [];\nNUz = [];\nEUpx = [];\nEUpy = [];\nEUpz = [];\nomLo = 2 * pi * P.min_mode_hz;\nfor family = 0:1\n [K, M] = triangle_fem(nodes, family, E, rho, nu, d);\n Rc = chol(M);\n Am = (Rc') \\ (K / Rc);\n Am = (Am + Am') / 2;\n [Uv, Sv, ~] = svd(Am);\n om2 = flipud(diag(Sv)); % ascending\n Vv = fliplr(Uv);\n om = sqrt(max(0, om2)) * P.stiffness_trim;\n keep = find(om > omLo & om < 2 * pi * P.max_mode_hz);\n if length(keep) == 0\n continue;\n end\n Phi = Rc \\ Vv(:, keep); % mass-normalized shapes, 3 nn x nk\n nk = length(keep);\n if family == 0\n ux = Phi(1:3:3 * nn, :);\n uy = Phi(2:3:3 * nn, :);\n uz = zeros(nn, nk);\n eux = (ux(1:nn - 1, :) + ux(2:nn, :)) / 2;\n euy = (uy(1:nn - 1, :) + uy(2:nn, :)) / 2;\n ax = eux .* (tx * ones(1, nk)) + euy .* (ty * ones(1, nk));\n upx = eux - ax .* (tx * ones(1, nk));\n upy = euy - ax .* (ty * ones(1, nk));\n upz = zeros(nel, nk);\n else\n ux = zeros(nn, nk);\n uy = zeros(nn, nk);\n uz = Phi(1:3:3 * nn, :);\n upx = zeros(nel, nk);\n upy = zeros(nel, nk);\n upz = (uz(1:nn - 1, :) + uz(2:nn, :)) / 2;\n end\n fs = [fs; om(keep) / (2 * pi)];\n fams = [fams; family * ones(nk, 1)];\n NUx = [NUx, ux];\n NUy = [NUy, uy];\n NUz = [NUz, uz];\n EUpx = [EUpx, upx];\n EUpy = [EUpy, upy];\n EUpz = [EUpz, upz];\nend\nn = length(fs);\nif n == 0\n error('No rod modes between min_mode_hz and max_mode_hz -- check geometry/material.');\nend\n\n[fs, order] = sort(fs);\nfs = reshape(fs, n, 1);\nfams = reshape(fams(order), n, 1);\nNUx = NUx(:, order);\nNUy = NUy(:, order);\nNUz = NUz(:, order);\nEUpx = EUpx(:, order);\nEUpy = EUpy(:, order);\nEUpz = EUpz(:, order);\n\n% ---- damping: material + radiation ----\nomega = 2 * pi * fs;\na = d / 2;\nVdip = 2 * pi * a^2 * el; % nel x 1 dipole volume per element\nup2 = EUpx.^2 + EUpy.^2 + EUpz.^2; % nel x n\nS2 = (Vdip.^2).' * up2; % 1 x n\nzetaRad = rho0 * omega.^3 .* reshape(S2, n, 1) / (24 * pi * c^3);\nzeta = 0.5 * (P.loss_factor + P.loss_slope * omega) + P.radiation_loss * zetaRad;\n\nbasis.n = n;\nbasis.f = fs;\nbasis.omega = omega;\nbasis.zeta = zeta;\nbasis.zetaRad = zetaRad;\nbasis.family = fams;\nbasis.nodeX = nodes(1, :);\nbasis.nodeY = nodes(2, :);\nbasis.sfrac = geo.sfrac;\nbasis.totalLen = geo.totalLen;\nbasis.NUx = NUx;\nbasis.NUy = NUy;\nbasis.NUz = NUz;\nbasis.elemX = ex;\nbasis.elemY = ey;\nbasis.elemTanX = tx;\nbasis.elemTanY = ty;\nbasis.elemLen = el;\nbasis.EUpx = EUpx;\nbasis.EUpy = EUpy;\nbasis.EUpz = EUpz;\nend\n"
38- },
39- {
40- "name": "triangle_defaults.m",
41- "content": "function P = triangle_defaults()\n% TRIANGLE_DEFAULTS Default parameters, mirroring the web app's params.ts.\n% The output duration is not a parameter: it is\n% (strike_count - 1) / strike_rate_hz + tail_s.\n\n% Rod geometry (equilateral, open at the bottom-left corner)\nP.side_length_m = 0.152;\nP.gap_m = 0.015;\nP.rod_diameter_m = 0.008;\n\n% Rod material (steel)\nP.density = 7850.0;\nP.youngs_modulus = 2.0e11;\nP.poisson_ratio = 0.30;\nP.stiffness_trim = 1.00;\n\n% Discretization / mode retention\nP.n_elements = 60;\nP.min_mode_hz = 20.0;\nP.max_mode_hz = 22000.0;\n\n% Damping\nP.loss_factor = 6.0e-5;\nP.loss_slope = 4.0e-10;\nP.radiation_loss = 1.0;\n\n% Beater (steel rod tip)\nP.beater_radius_m = 0.0045;\nP.beater_mass_kg = 0.022;\nP.beater_youngs = 2.0e11;\nP.restitution = 0.60;\n\n% Strike / roll gesture\nP.strike_velocity = 1.5;\nP.strike_pos = 0.17; % arclength fraction: middle of the base\nP.strike_angle_deg = 20.0; % 0 = in the plane of the triangle, 90 = normal to it\nP.strike_count = 12;\nP.strike_rate_hz = 12.0;\nP.roll_pos2 = 0.40; % second contact point of a roll (the right leg)\nP.timing_jitter_ms = 6.0;\nP.velocity_spread = 0.15;\nP.seed = 1;\n\n% Radiation / listener\nP.direct_gain = 1.0;\nP.listener_x_m = 0.20;\nP.listener_y_m = 0.10;\nP.listener_z_m = 0.40;\nP.speed_of_sound = 343.0;\nP.air_density = 1.204;\n\n% Rendering / output\nP.sample_rate = 48000;\nP.oversample = 2;\nP.tail_s = 5.0;\nP.pre_delay_ms = 3.0;\nP.normalize = 1;\nP.peak_dbfs = -3.0;\nP.full_scale_pa = 2.0;\nend\n"
42- },
43- {
44- "name": "triangle_demo.m",
45- "content": "% Render a triangle roll with default parameters and print diagnostics.\n% Set P.strike_count = 1 for a single strike.\n% Run with: numbl run triangle_demo.m\n\nP = triangle_defaults();\ntic;\n[y, basis, dinfo] = triangle_render(P);\nelapsed = toc;\n\nfprintf('modes kept : %d (%d in-plane, %d out-of-plane)\\n', basis.n, dinfo.modesIn, dinfo.modesOut);\nfprintf('mode freq : ');\nfprintf('%.0f ', basis.f);\nfprintf('Hz\\n');\nfprintf('strikes : %d at %.1f Hz (mean v %.2f m/s, mean contact %.0f us, peak force %.0f N)\\n', ...\n dinfo.strikeCount, P.strike_rate_hz, dinfo.meanV, dinfo.meanTauUs, dinfo.peakForceN);\nfprintf('peak press.: %.4f Pa (%.1f dB SPL) at %.0f cm\\n', ...\n dinfo.peakPa, dinfo.peakSplDb, dinfo.listenerDistM * 100);\nfprintf('length : %.2f s @ %d Hz\\n', length(y) / P.sample_rate, P.sample_rate);\nfprintf('render time: %.2f s\\n', elapsed);\nif dinfo.listenerInside\n fprintf('WARNING : listener is at or inside the rod\\n');\nend\n"
46- },
47- {
48- "name": "triangle_fem.m",
49- "content": "function [K, M] = triangle_fem(nodes, family, E, rho, nu, d)\n% TRIANGLE_FEM Free-free frame system for one mode family of a bent rod.\n%\n% The rod centerline is the open polyline `nodes` (2 x nn, metres, in the\n% x-y plane), discretized into nn - 1 two-node elements. Because the geometry\n% is planar, the six degrees of freedom per node decouple exactly into two\n% independent families of three, and each family is assembled separately:\n% family 0 (in-plane): ux, uy, theta_z (in-plane bending + stretching)\n% family 1 (out-of-plane): uz, theta_x, theta_y (out-of-plane bending + torsion)\n% The bends are where the families earn their keep: within one straight leg\n% out-of-plane bending and torsion are independent, but a change of direction\n% rotates the (torsion, bending-rotation) pair, so the corner converts one\n% into the other. That conversion is what the per-element rotation T applies,\n% and it is what makes a triangle sound different from a straight rod.\n%\n% Elements are Euler-Bernoulli (cubic Hermite bending, linear axial and\n% torsion interpolation) with consistent mass including rotary inertia.\n% Shear deformation is neglected, which leaves the highest retained modes a\n% few percent sharp for a rod this thick; see the README.\n%\n% Sign conventions: in-plane bending uses theta_z = +dv/dx, giving the\n% textbook stiffness block; out-of-plane bending uses theta_y = -dw/dx, so\n% the odd coupling terms flip sign (the similarity transform Sf below).\n%\n% Returns K and M, each 3 nn x 3 nn, symmetric; M is positive definite, so\n% chol(M) reduces the generalized eigenproblem to a symmetric standard one.\n\nnn = size(nodes, 2);\nndof = 3 * nn;\nK = zeros(ndof);\nM = zeros(ndof);\n\nAa = pi * d^2 / 4; % cross-section area\nIz = pi * d^4 / 64; % bending second moment (circular)\nJp = pi * d^4 / 32; % polar moment = torsion constant (circular)\nGs = E / (2 * (1 + nu));\n\nfor e = 1:nn - 1\n dx = nodes(1, e + 1) - nodes(1, e);\n dy = nodes(2, e + 1) - nodes(2, e);\n l = sqrt(dx^2 + dy^2);\n c = dx / l;\n s = dy / l;\n\n % bending block on the (transverse, rotation) pairs of both nodes\n Kb = (E * Iz / l^3) * ...\n [12, 6 * l, -12, 6 * l; ...\n 6 * l, 4 * l^2, -6 * l, 2 * l^2; ...\n -12, -6 * l, 12, -6 * l; ...\n 6 * l, 2 * l^2, -6 * l, 4 * l^2];\n Mb = (rho * Aa * l / 420) * ...\n [156, 22 * l, 54, -13 * l; ...\n 22 * l, 4 * l^2, 13 * l, -3 * l^2; ...\n 54, 13 * l, 156, -22 * l; ...\n -13 * l, -3 * l^2, -22 * l, 4 * l^2] ...\n + (rho * Iz / (30 * l)) * ...\n [36, 3 * l, -36, 3 * l; ...\n 3 * l, 4 * l^2, -3 * l, -l^2; ...\n -36, -3 * l, 36, -3 * l; ...\n 3 * l, -l^2, -3 * l, 4 * l^2];\n\n Ke = zeros(6, 6);\n Me = zeros(6, 6);\n if family == 0\n % local order (u1, v1, t1, u2, v2, t2): axial on u, bending on (v, t)\n Ka2 = (E * Aa / l) * [1, -1; -1, 1];\n Ma2 = (rho * Aa * l / 6) * [2, 1; 1, 2];\n ia = [1 4];\n ib = [2 3 5 6];\n Ke(ia, ia) = Ka2;\n Me(ia, ia) = Ma2;\n Ke(ib, ib) = Kb;\n Me(ib, ib) = Mb;\n T = [c, s, 0; -s, c, 0; 0, 0, 1];\n else\n % local order (w1, tx1, ty1, w2, tx2, ty2): torsion on tx (about the leg\n % axis), bending on (w, ty) with theta_y = -dw/dx\n Kt2 = (Gs * Jp / l) * [1, -1; -1, 1];\n Mt2 = (rho * Jp * l / 6) * [2, 1; 1, 2];\n Sf = diag([1, -1, 1, -1]);\n ia = [2 5];\n ib = [1 3 4 6];\n Ke(ia, ia) = Kt2;\n Me(ia, ia) = Mt2;\n Ke(ib, ib) = Sf * Kb * Sf;\n Me(ib, ib) = Sf * Mb * Sf;\n T = [1, 0, 0; 0, c, s; 0, -s, c];\n end\n\n TT = blkdiag(T, T);\n Kg = TT' * Ke * TT;\n Mg = TT' * Me * TT;\n dofs = [3 * e - 2, 3 * e - 1, 3 * e, 3 * e + 1, 3 * e + 2, 3 * e + 3];\n K(dofs, dofs) = K(dofs, dofs) + Kg;\n M(dofs, dofs) = M(dofs, dofs) + Mg;\nend\nend\n"
50- },
51- {
52- "name": "triangle_geometry.m",
53- "content": "function G = triangle_geometry(P)\n% TRIANGLE_GEOMETRY Node polyline of the open triangle centerline.\n%\n% An equilateral triangle of side side_length_m in the x-y plane, apex up:\n% base from V2 (bottom left) to V3 (bottom right) along y = 0, apex V1 above\n% the origin. The rod is open at the bottom-left corner: gap_m of arclength\n% is removed there, half from the base and half from the left leg, so the\n% centerline runs A -> V3 -> V1 -> B counterclockwise, where A and B sit\n% gap_m / 2 from V2 along their legs. Total rod length is 3 L - gap.\n%\n% The polyline is discretized into about n_elements equal-length elements,\n% with the two closed corners landing exactly on nodes. Arclength fraction 0\n% is the free end on the base, 1 the free end on the left leg; strike_pos and\n% roll_pos2 are expressed in this coordinate.\n%\n% Returns G with:\n% nodes 2 x nn node coordinates (m)\n% sfrac 1 x nn arclength fraction of each node\n% totalLen rod arclength (m)\n% nn, ne node count and per-leg element counts (1 x 3)\n\nL = P.side_length_m;\ng = min(P.gap_m, 0.8 * L);\nV1 = [0; L * sqrt(3) / 2];\nV2 = [-L / 2; 0];\nV3 = [L / 2; 0];\ne1 = (V3 - V2) / L;\ne3 = (V2 - V1) / L;\nA = V2 + e1 * (g / 2);\nB = V2 - e3 * (g / 2);\n\nstarts = [A, V3, V1];\nstops = [V3, V1, B];\nlens = [L - g / 2, L, L - g / 2];\ntotal = lens(1) + lens(2) + lens(3);\n\nne = zeros(1, 3);\nfor k = 1:3\n ne(k) = max(2, round(P.n_elements * lens(k) / total));\nend\nnn = ne(1) + ne(2) + ne(3) + 1;\n\nnodes = zeros(2, nn);\nsarc = zeros(1, nn);\nnodes(:, 1) = A;\nidx = 1;\nacc = 0;\nfor k = 1:3\n for j = 1:ne(k)\n idx = idx + 1;\n frac = j / ne(k);\n nodes(:, idx) = starts(:, k) + (stops(:, k) - starts(:, k)) * frac;\n sarc(idx) = acc + lens(k) * frac;\n end\n acc = acc + lens(k);\nend\n\nG.nodes = nodes;\nG.sfrac = sarc / total;\nG.totalLen = total;\nG.nn = nn;\nG.ne = ne;\nend\n"
54- },
55- {
56- "name": "triangle_output.m",
57- "content": "function [y, peakPa] = triangle_output(P, p, oversample)\n% TRIANGLE_OUTPUT Anti-aliased decimation and output conditioning.\n% Kaiser-windowed-sinc FIR lowpass at 0.45 * fs_lo, decimate by oversample,\n% pre-roll, tail fade, then either normalize to peak_dbfs or divide by\n% full_scale_pa. Identical in structure to the woodblock and maraca output\n% stages.\n\npeakPa = max(abs(p));\nR = oversample;\nT = length(p);\n\nif R > 1\n taps = 12 * R + 1;\n M2 = (taps - 1) / 2;\n fcNorm = 0.45 / R; % cycles per input sample\n m = (0:taps - 1) - M2;\n h = sin(2 * pi * fcNorm * m) ./ (pi * m);\n h(M2 + 1) = 2 * fcNorm; % sinc value at m = 0\n r = m / M2;\n beta = 8;\n h = h .* (tr_i0(beta * sqrt(max(0, 1 - r.^2))) / tr_i0(beta));\n h = h / sum(h); % unity DC gain\n pf = conv(p, h);\n outLen = floor(T / R);\n y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen); % centered pick\nelse\n y = p;\nend\n\n% Pre-roll, then truncate to the nominal duration\npadN = floor(P.pre_delay_ms * 1e-3 * P.sample_rate);\nK = max(1, round(P.strike_count));\nduration = (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;\ntotal = floor(duration * P.sample_rate);\nif padN > 0\n y = [zeros(1, padN) y];\nend\ny = y(1:min(length(y), total));\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_pa;\nend\nend\n\nfunction v = tr_i0(x)\n% Modified Bessel function of the first kind, order zero (series expansion).\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"
58- },
59- {
60- "name": "triangle_radiate.m",
61- "content": "function [p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, fsHi)\n% TRIANGLE_RADIATE Modal synthesis and free-field pressure at the listener.\n%\n% Two contributions:\n%\n% 1. Modal ringing. A strike with impulse J at arclength position s gives\n% every mode a modal-velocity jump J S(omega tau) (phi(s) . dhat), where\n% S(x) = cos(x/2) / (1 - (x/pi)^2) is the spectrum of the half-sine\n% contact pulse and dhat the strike direction (strike_angle_deg between\n% the leg's in-plane normal and the plane normal, so one strike excites\n% both mode families). The rod radiates through its transverse motion:\n% each element is a point dipole of moment Vdip u_perp qdot (Vdip from\n% triangle_basis.m), and per mode the element dipoles are summed\n% coherently with their relative retardations,\n% G = sum_e (i om rho0 / 4 pi) Vdip_e (u_e . rhat_e)\n% (i k + 1/r_e) exp(-i k (r_e - rcen)) / r_e,\n% evaluated at the mode's centre frequency (modes are narrow-band). The\n% common carrier exp(-i k rcen) is stripped, so relative phase across the\n% triangle is kept while the overall propagation delay is not double\n% counted. This sum is where the interference across the frame lives:\n% moving the listener re-runs only this stage.\n%\n% 2. Direct strike clicks. Modes above max_mode_hz are truncated; their\n% collective off-resonance response is the immediate local motion of the\n% rod around the contact. It is reinstated from the classic driving-point\n% impedance of an infinite Euler beam, Z = 2 rho A c_B (1 + i) with\n% c_B = sqrt(om kappa cL), evaluated at the pulse's own centre frequency\n% om* = pi / tau: the contact force divided by |Z| is the local velocity,\n% which radiates as a dipole of length half a flexural wavelength along\n% the strike direction. Scaled by direct_gain. This is the least\n% principled term of the model, and direct_gain = 0 turns it off.\n%\n% The synthesis is exact in time: between strikes every mode evolves as\n% exp(lambda t), evaluated in chunks over each inter-strike segment (the\n% chunking only bounds memory; nothing is integrated).\n%\n% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.\n\nc = P.speed_of_sound;\nrho0 = P.air_density;\nK = max(1, round(P.strike_count));\nduration = (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;\nT = floor(duration * fsHi);\nN = basis.n;\nI = length(tstr);\nnel = length(basis.elemLen);\nnn = length(basis.sfrac);\n\nlx = P.listener_x_m;\nly = P.listener_y_m;\nlz = P.listener_z_m;\n\n% ---- element-to-listener geometry ----\nrex = lx - basis.elemX; % nel x 1\nrey = ly - basis.elemY;\nrez = lz * ones(nel, 1);\nre = sqrt(rex.^2 + rey.^2 + rez.^2);\nlistenerInside = min(re) < 0.02;\nre = max(re, 0.01);\nrhx = rex ./ re;\nrhy = rey ./ re;\nrhz = rez ./ re;\ncx = sum(basis.elemX .* basis.elemLen) / sum(basis.elemLen);\ncy = sum(basis.elemY .* basis.elemLen) / sum(basis.elemLen);\nrcen = sqrt((lx - cx)^2 + (ly - cy)^2 + lz^2);\n\n% ---- per-mode radiation transfer and complex pole ----\na = P.rod_diameter_m / 2;\nVdip = 2 * pi * a^2 * basis.elemLen; % nel x 1\nG = zeros(N, 1) + 0i;\nlam = zeros(N, 1) + 0i;\nphv = zeros(N, 1) + 0i; % modal-velocity phasor per unit jump\nfor m = 1:N\n om = basis.omega(m);\n zt = basis.zeta(m);\n k = om / c;\n proj = basis.EUpx(:, m) .* rhx + basis.EUpy(:, m) .* rhy + basis.EUpz(:, m) .* rhz;\n ph = exp(-1i * k * (re - rcen));\n G(m) = (1i * om * rho0 / (4 * pi)) * sum(Vdip .* proj .* (1i * k + 1 ./ re) .* ph ./ re);\n omd = om * sqrt(max(1e-12, 1 - zt^2));\n lam(m) = -zt * om + 1i * omd;\n phv(m) = 1 + 1i * zt * om / omd; % qdot(t) = Re[phv e^{lam t}] per unit jump\nend\n\n% ---- per-strike modal weights: A(n, i) modal-velocity-jump phasors ----\nA = zeros(N, I) + 0i;\ndirx = zeros(1, I);\ndiry = zeros(1, I);\ndirz = zeros(1, I);\nhitx = zeros(1, I);\nhity = zeros(1, I);\nang = P.strike_angle_deg * pi / 180;\nfor i = 1:I\n s = min(0.999, max(0.001, sstr(i)));\n idx = sum(basis.sfrac <= s);\n idx = min(max(idx, 1), nn - 1);\n w = (s - basis.sfrac(idx)) / (basis.sfrac(idx + 1) - basis.sfrac(idx));\n uxh = basis.NUx(idx, :) + w * (basis.NUx(idx + 1, :) - basis.NUx(idx, :));\n uyh = basis.NUy(idx, :) + w * (basis.NUy(idx + 1, :) - basis.NUy(idx, :));\n uzh = basis.NUz(idx, :) + w * (basis.NUz(idx + 1, :) - basis.NUz(idx, :));\n % strike direction: the leg's inward in-plane normal, tilted out of plane\n ex = basis.elemTanX(idx);\n ey = basis.elemTanY(idx);\n dirx(i) = cos(ang) * (-ey);\n diry(i) = cos(ang) * ex;\n dirz(i) = sin(ang);\n hitx(i) = basis.nodeX(idx) + w * (basis.nodeX(idx + 1) - basis.nodeX(idx));\n hity(i) = basis.nodeY(idx) + w * (basis.nodeY(idx + 1) - basis.nodeY(idx));\n proj = uxh * dirx(i) + uyh * diry(i) + uzh * dirz(i); % 1 x N\n x = basis.omega * taustr(i); % N x 1\n den = 1 - (x / pi).^2;\n S = cos(x / 2) ./ den;\n sing = abs(den) < 1e-6;\n S(sing) = pi / 4; % limit at x -> pi\n A(:, i) = phv .* (jstr(i) * S) .* reshape(proj, N, 1);\nend\n\n% ---- segment synthesis: exact damped-exponential evolution, chunked ----\nks = min(T, max(1, round(tstr * fsHi) + 1));\np = zeros(1, T);\ns = zeros(N, 1) + 0i;\nGt = G.';\ndt = 1 / fsHi;\nCH = 32768; % chunk length: bounds the exp matrix\nprevK = ks(1);\ni = 1;\nwhile i <= I\n s = s + A(:, i);\n while i < I && ks(i + 1) == ks(i)\n i = i + 1;\n s = s + A(:, i);\n end\n if i < I\n nextK = ks(i + 1);\n else\n nextK = T + 1;\n end\n pos = prevK;\n while pos < nextK\n len = min(CH, nextK - pos);\n seg = real(Gt * ((s * ones(1, len)) .* exp(lam * ((0:len - 1) * dt))));\n p(pos:pos + len - 1) = p(pos:pos + len - 1) + seg;\n s = s .* exp(lam * (len * dt));\n pos = pos + len;\n end\n prevK = nextK;\n i = i + 1;\nend\n\n% ---- direct strike clicks (driving-point velocity radiating as a dipole) ----\nAarea = pi * (2 * a)^2 / 4;\nkapcL = (a / 2) * sqrt(P.youngs_modulus / P.density); % kappa * cL\ndirectPeak = 0;\nfor i = 1:I\n rsx = lx - hitx(i);\n rsy = ly - hity(i);\n rs = max(0.01, sqrt(rsx^2 + rsy^2 + lz^2));\n cosg = (rsx * dirx(i) + rsy * diry(i) + lz * dirz(i)) / rs;\n nP = max(2, round(taustr(i) * fsHi));\n k0 = ks(i) + round((rs - rcen) / c * fsHi);\n if k0 < 1 || k0 + nP > T\n continue;\n end\n omStar = pi / taustr(i);\n cB = sqrt(omStar * kapcL);\n Zb = 2 * sqrt(2) * P.density * Aarea * cB; % |2 rho A cB (1 + i)|\n lamB = 2 * pi * cB / omStar; % flexural wavelength at om*\n mu = 2 * pi * a^2 * (lamB / 2) / Zb; % dipole moment per unit force\n amp = P.direct_gain * fstr(i) * mu * rho0 * omStar^2 * cosg / (4 * pi * rs * c);\n pulse = amp * sin(pi * (0:nP) / nP);\n p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;\n if abs(amp) > directPeak\n directPeak = abs(amp);\n end\nend\n\nradinfo.gainAbs = abs(G);\nradinfo.directPeakPa = directPeak;\nradinfo.listenerInside = listenerInside;\nradinfo.listenerDistM = rcen;\nend\n"
62- },
63- {
64- "name": "triangle_render.m",
65- "content": "function [y, basis, dinfo] = triangle_render(P)\n% TRIANGLE_RENDER Full render: basis + strikes -> radiate -> output.\n% Convenience composition of the stage functions; the web app calls the\n% stages separately so it can cache intermediates (the eigensolve above all).\n\nbasis = triangle_basis(P);\nfsHi = P.sample_rate * P.oversample;\n[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P);\n[p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, fsHi);\n[y, peakPa] = triangle_output(P, p, P.oversample);\n\ndinfo.strikeCount = sinfo.count;\ndinfo.meanV = sinfo.meanV;\ndinfo.maxV = sinfo.maxV;\ndinfo.meanTauUs = sinfo.meanTauUs;\ndinfo.peakForceN = sinfo.peakForceN;\ndinfo.durationS = sinfo.durationS;\ndinfo.modesIn = sum(basis.family == 0);\ndinfo.modesOut = sum(basis.family == 1);\ndinfo.f1 = basis.f(1);\ndinfo.peakPa = peakPa;\ndinfo.peakSplDb = 20 * log10(peakPa / 20e-6 + 1e-30);\ndinfo.directPeakPa = radinfo.directPeakPa;\ndinfo.listenerInside = radinfo.listenerInside;\ndinfo.listenerDistM = radinfo.listenerDistM;\ndinfo.sinfo = sinfo;\nend\n"
66- },
67- {
68- "name": "triangle_rng.m",
69- "content": "function [r, state] = triangle_rng(state, n)\n% TRIANGLE_RNG Portable seeded uniform generator (32-bit LCG).\n% Returns n uniforms in (0, 1) as a column vector, plus the advanced state.\n% Deliberately not MATLAB's rand: the same seed must produce the identical\n% roll pattern in the numbl CLI, Node, and the browser, and every double\n% here stays below 2^53 so the arithmetic is exact everywhere.\nr = zeros(n, 1);\nfor i = 1:n\n state = mod(state * 1664525 + 1013904223, 4294967296);\n r(i) = (state + 0.5) / 4294967296;\nend\nend\n"
70- },
71- {
72- "name": "triangle_strikes.m",
73- "content": "function [tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P)\n% TRIANGLE_STRIKES Beater strike events: a single strike, or a seeded roll.\n%\n% strike_count = 1 gives one deterministic strike of strike_velocity at\n% strike_pos at t = 0; no random numbers are drawn. strike_count > 1 gives a\n% roll: nominally periodic hits at strike_rate_hz, alternating between\n% strike_pos and roll_pos2 (the beater rattling in a corner strikes the two\n% legs that form it in turn), with Gaussian timing jitter (timing_jitter_ms)\n% and lognormal velocity spread (velocity_spread). The first hit stays at\n% t = 0 so the onset is fixed. All randomness comes from triangle_rng seeded\n% with P.seed, so a given seed reproduces the identical roll in the CLI,\n% Node, and the browser.\n%\n% Each strike is a Hertz contact of the beater tip (a steel sphere of radius\n% beater_radius_m) on the rod:\n% k_h = (4/3) E* sqrt(R_eff), R_eff = sqrt(r_b R''),\n% 1/R'' = 1/r_b + 2/d (sphere on cylinder: R_eff is the\n% geometric mean of the two relative\n% principal curvatures)\n% tau = 3.218 (m*^2 / (k_h^2 v))^{1/5} (Hertz contact time)\n% J = (1 + restitution) m* v (momentum transferred)\n% Fpk = pi J / (2 tau) (half-sine pulse peak)\n% The reduced mass m* pairs the beater with the rod mass engaged during the\n% contact. The engaged length is estimated from the flexural group velocity\n% at the pulse's own bandwidth, l = 2 sqrt(pi kappa cL tau), which depends\n% on tau, so a short fixed-point iteration closes the loop. This is an\n% estimate, not a contact simulation: what it captures is that a heavy\n% beater lengthens the contact and darkens the strike, while a light beater\n% approaches the free-sphere limit and tau ~ v^(-1/5).\n%\n% Returns row vectors sorted by time: tstr (s), vstr (m/s), jstr (N s),\n% taustr (s), fstr (peak N), sstr (arclength fraction of the contact), and\n% sinfo with summary statistics.\n\nK = max(1, round(P.strike_count));\nrate = max(1e-3, P.strike_rate_hz);\nduration = (K - 1) / rate + P.tail_s;\n\nE = P.youngs_modulus;\nrho = P.density;\nd = P.rod_diameter_m;\nnub = 0.30; % beater Poisson ratio, fixed\nEstar = 1 / ((1 - nub^2) / P.beater_youngs + (1 - P.poisson_ratio^2) / E);\nrb = P.beater_radius_m;\nRpp = 1 / (1 / rb + 2 / d);\nReff = sqrt(rb * Rpp);\nkh = (4 / 3) * Estar * sqrt(Reff);\nkapcL = (d / 4) * sqrt(E / rho); % kappa * cL: flexural dispersion constant\nAarea = pi * d^2 / 4;\n\nstate = mod(floor(abs(P.seed)), 4294967296);\nz1 = zeros(K, 1);\nz2 = zeros(K, 1);\nif K > 1\n [u, state] = triangle_rng(state, 4 * K);\n U = reshape(u, K, 4);\n z1 = sqrt(-2 * log(U(:, 1))) .* cos(2 * pi * U(:, 2));\n z2 = sqrt(-2 * log(U(:, 3))) .* cos(2 * pi * U(:, 4));\nend\n\ntstr = zeros(1, K);\nvstr = zeros(1, K);\njstr = zeros(1, K);\ntaustr = zeros(1, K);\nfstr = zeros(1, K);\nsstr = zeros(1, K);\nfor i = 1:K\n t = (i - 1) / rate;\n if i > 1\n t = max(0, t + P.timing_jitter_ms * 1e-3 * z1(i));\n end\n v = P.strike_velocity;\n if K > 1\n v = v * exp(P.velocity_spread * z2(i));\n end\n v = max(0.01, v);\n\n % engaged-mass fixed point for the Hertz contact time\n tau = 30e-6;\n mstar = P.beater_mass_kg;\n for it = 1:4\n leng = 2 * sqrt(pi * kapcL * tau);\n mloc = rho * Aarea * leng;\n mstar = P.beater_mass_kg * mloc / (P.beater_mass_kg + mloc);\n tau = 3.218 * (mstar^2 / (kh^2 * v))^(1 / 5);\n end\n\n tstr(i) = t;\n vstr(i) = v;\n jstr(i) = (1 + P.restitution) * mstar * v;\n taustr(i) = tau;\n fstr(i) = pi * jstr(i) / (2 * tau);\n if mod(i, 2) == 1\n sstr(i) = P.strike_pos;\n else\n sstr(i) = P.roll_pos2;\n end\nend\n\n[tstr, order] = sort(tstr);\ntstr = reshape(tstr, 1, K);\nvstr = reshape(vstr(order), 1, K);\njstr = reshape(jstr(order), 1, K);\ntaustr = reshape(taustr(order), 1, K);\nfstr = reshape(fstr(order), 1, K);\nsstr = reshape(sstr(order), 1, K);\n\nsinfo.count = K;\nsinfo.meanV = sum(vstr) / K;\nsinfo.maxV = max(vstr);\nsinfo.meanTauUs = sum(taustr) / K * 1e6;\nsinfo.peakForceN = max(fstr);\nsinfo.durationS = duration;\nend\n"
74- },
75- {
76- "name": "triangle_validate.m",
77- "content": "% Validate the triangle model. There is no external reference implementation,\n% so the checks are (a) closed-form identities the physics must satisfy,\n% (b) structural identities of the FEM eigensolve, (c) exact reproducibility\n% of the seeded roll process, and (d) pinned end-to-end diagnostics from the\n% default configuration, recorded when the model was built and treated as a\n% regression reference since.\n% Run with: numbl run triangle_validate.m\n\nfails = 0;\n\n% ---- 1. straight free-free rod against closed forms ----\n% The same element routines that build the triangle, assembled on a straight\n% polyline, must reproduce the textbook free-free bar: bending frequencies\n% (beta L)^2 kappa cL / (2 pi L^2), the first axial mode cL / 2L, the first\n% torsional mode cT / 2L, and exactly three rigid-body modes per family.\nE = 2e11; rho = 7850; nu = 0.30; d = 0.008; L = 0.44;\nnodes = [linspace(0, L, 61); zeros(1, 61)];\ncL = sqrt(E / rho);\nkap = d / 4;\nblTab = [4.73004074 7.85320462 10.99560784];\nfbClosed = blTab.^2 * kap * cL / (2 * pi * L^2);\nfaClosed = cL / (2 * L);\nftClosed = sqrt(E / (2 * (1 + nu)) / rho) / (2 * L);\nf0 = tv_eigfreqs(nodes, 0, E, rho, nu, d);\nf1 = tv_eigfreqs(nodes, 1, E, rho, nu, d);\n% the rigid-body zeros come back as numerical noise well below the first\n% elastic mode (185 Hz here); 5 Hz separates them cleanly\nnRigid0 = sum(f0 < 5);\nnRigid1 = sum(f1 < 5);\ne0 = f0(f0 > 1);\ne1 = f1(f1 > 1);\nerrB = max(abs(e0(1:3).' - fbClosed) ./ fbClosed);\n% rotary inertia lowers the FEM bending values slightly below the closed\n% form, which neglects it; 0.5% covers that plus discretization\nerrA = abs(min(abs(e0 - faClosed)) / faClosed);\nerrT = abs(min(abs(e1 - ftClosed)) / ftClosed);\nok = nRigid0 == 3 && nRigid1 == 3 && errB < 5e-3 && errA < 2e-3 && errT < 2e-3;\nfails = fails + (~ok);\nfprintf('straight rod: bend err %.1e, axial %.1e, torsion %.1e, rigid %d/%d %s\\n', ...\n errB, errA, errT, nRigid0, nRigid1, tv_verdict(ok));\n\n% ---- 2. eigensolve identities on the triangle itself ----\n% The Cholesky reduction must hand back mass-orthonormal shapes that\n% diagonalize the stiffness: phi' M phi = I and phi' K phi = diag(om^2).\nP = triangle_defaults();\ngeo = triangle_geometry(P);\n[K, M] = triangle_fem(geo.nodes, 0, P.youngs_modulus, P.density, P.poisson_ratio, P.rod_diameter_m);\nRc = chol(M);\nAm = (Rc') \\ (K / Rc);\nAm = (Am + Am') / 2;\n[Uv, Sv, ~] = svd(Am);\nom2 = flipud(diag(Sv));\nVv = fliplr(Uv);\nPhi = Rc \\ Vv;\neOrth = max(max(abs(Phi' * M * Phi - eye(length(om2)))));\nKP = Phi' * K * Phi;\neDiag = max(max(abs(KP - diag(om2)))) / max(om2);\nok = eOrth < 1e-8 && eDiag < 1e-10;\nfails = fails + (~ok);\nfprintf('eigensolve: orthonormality %.1e, diagonalization %.1e %s\\n', eOrth, eDiag, tv_verdict(ok));\n\n% ---- 3. Hertz contact-time scaling tau ~ v^(-1/5) ----\n% With a beater much lighter than the engaged rod mass the reduced mass is a\n% constant and the free-sphere scaling must emerge. A single strike, so the\n% check does not depend on the default gesture or on the velocity spread.\nQ = P;\nQ.strike_count = 1;\nQ.beater_mass_kg = 5e-4;\n[~, ~, ~, tau1, ~, ~, ~] = triangle_strikes(Q);\nQ.strike_velocity = 2 * P.strike_velocity;\n[~, ~, ~, tau2, ~, ~, ~] = triangle_strikes(Q);\nratio = tau2(1) / tau1(1);\nok = abs(ratio - 2^(-1 / 5)) < 0.03;\nfails = fails + (~ok);\nfprintf('hertz scaling: tau ratio %.4f (target %.4f) %s\\n', ratio, 2^(-1 / 5), tv_verdict(ok));\n\n% ---- 4. seeded roll is exactly reproducible ----\nQ = P;\nQ.strike_count = 12;\n[t1r, v1r, j1r, ~, ~, s1r, g1] = triangle_strikes(Q);\n[t2r, v2r, j2r, ~, ~, s2r, g2] = triangle_strikes(Q);\nsame = max(abs(t1r - t2r)) == 0 && max(abs(v1r - v2r)) == 0 && max(abs(s1r - s2r)) == 0;\nok = same ...\n && g1.count == 12 ...\n && abs(g1.meanV - 1.501211020) < 1e-6 ...\n && abs(sum(j1r) - 0.383069278968) < 1e-9 ...\n && abs(t1r(2) - 0.087732781124) < 1e-9;\nfails = fails + (~ok);\nfprintf('roll (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\\n', ...\n Q.seed, g1.count, g1.meanV, sum(j1r), tv_verdict(ok));\n\n% ---- 5. far-field falloff: -6 dB per doubling of distance ----\n% RMS over the render, with the two listeners on the same ray from the\n% origin so every element dipole keeps its directivity factor: off the ray\n% the in-plane dipoles fall faster than 1/r, which is physics rather than\n% error. Modes are capped at 8 kHz so 2 m is beyond the Rayleigh distance\n% of the frame; at 20 kHz it is not, and Fresnel-zone interference shifts\n% the ratio by a few tenths of a dB, again physics rather than error.\nQ = P;\nQ.normalize = 0;\nQ.strike_count = 1;\nQ.tail_s = 1.0;\nQ.max_mode_hz = 8000;\nQ.listener_x_m = 0.6; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;\n[y1, ~, ~] = triangle_render(Q);\nQ.listener_x_m = 1.2; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;\n[y2, ~, ~] = triangle_render(Q);\ndrop = 10 * log10(sum(y2.^2) / sum(y1.^2));\nok = abs(drop + 6.02) < 0.2;\nfails = fails + (~ok);\nfprintf('falloff 2 m -> 4 m: %.2f dB rms (target -6.02) %s\\n', drop, tv_verdict(ok));\n\n% ---- 6. pinned end-to-end diagnostics at the defaults ----\n[y, basis, dinfo] = triangle_render(P);\nok = abs(dinfo.peakPa - 13.769445) < 0.02 ...\n && abs(basis.f(1) - 163.9200) < 0.1 ...\n && basis.n == 40;\nfails = fails + (~ok);\nfprintf('pinned render: peak %.4f Pa (target %.4f), mode1 %.2f Hz (target %.2f), %d modes %s\\n', ...\n dinfo.peakPa, 13.769445, basis.f(1), 163.9200, basis.n, tv_verdict(ok));\n\n% ---- 7. the ring decays ----\n% A triangle rings for seconds, so the tail is not silent; but the energy in\n% the last half second must sit well below that of the opening half second,\n% which at the defaults carries the first six hits of the roll.\nfs = P.sample_rate;\nn = length(y);\nhead = y(1:round(0.5 * fs));\ntail = y(n - round(0.5 * fs):n);\ndropE = 10 * log10(sum(tail.^2) / sum(head.^2));\nok = dropE < -25;\nfails = fails + (~ok);\nfprintf('ring decay: last/first half second %.1f dB %s\\n', dropE, tv_verdict(ok));\n\nif fails > 0\n fprintf('%d CHECK(S) FAILED\\n', fails);\nelse\n fprintf('all checks passed\\n');\nend\n\nfunction f = tv_eigfreqs(nodes, family, E, rho, nu, d)\n[K, M] = triangle_fem(nodes, family, E, rho, nu, d);\nRc = chol(M);\nAm = (Rc') \\ (K / Rc);\nAm = (Am + Am') / 2;\n[~, Sv, ~] = svd(Am);\nf = sort(sqrt(diag(Sv))) / (2 * pi);\nend\n\nfunction s = tv_verdict(ok)\nif ok\n s = 'ok';\nelse\n s = 'MISMATCH';\nend\nend\n"
78- }
79- ]
80-}
instrument/.gitignoreadded+1−0View file
@@ -0,0 +1 @@
1+node_modules/
instrument/core.tsadded+573−0View file
@@ -0,0 +1,573 @@
1+// The triangle as a track-render instrument: pure rendering logic with no
2+// DOM or Worker dependency, shared by the browser worker (worker.ts) and
3+// the node tests (test.ts). Two request shapes arrive here:
4+//
5+// - the track path (midiboot): every note of a MIDI track. Notes are
6+// grouped by cache key (pitch and bucketed velocity; a strike ignores
7+// note duration, ring-out is the tail_s knob), each distinct sound is
8+// rendered once through the staged pipeline, and copies are summed at
9+// the onsets. The FEM eigensolve depends only on the pitch-scaled
10+// geometry, so it is cached per pitch and shared across velocities,
11+// knob values, and repeats.
12+//
13+// - the explorer path (the triangle web app, wantDiagnostics set): the
14+// full parameter surface and live-edited .m sources drive one staged
15+// render with the diagnostics payload the app's plots consume, exactly
16+// the arrangement of the app's former render worker.
17+//
18+// Stage intermediates cross executeCode calls as numbl RuntimeValues
19+// (column-major Float64).
20+
21+import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
22+import type { MFile, TrackRenderRequest } from './types.ts';
23+import { TRIANGLE_KNOBS, type TriangleVariant } from './knobs.ts';
24+import { NoteCache } from './noteCache.ts';
25+
26+export { TRIANGLE_KNOBS, type TriangleVariant };
27+
28+/** The environment the host provides: sources, progress, cancellation, and
29+ * a macrotask yield so a browser worker can receive a cancel message. */
30+export interface RenderEnv {
31+ mfiles: MFile[];
32+ post: (frac: number, stage: string) => void;
33+ /** Throws to abandon the render; called between stages and notes. */
34+ checkCancel: () => void;
35+ yieldToHost: () => Promise<void>;
36+}
37+
38+// --- diagnostics payload for the explorer path (the web app's plots) ---
39+
40+export interface TriangleExploreDiag {
41+ quality: 'preview' | 'full';
42+ strikeT: Float32Array;
43+ strikeV: Float32Array;
44+ strikeS: Float32Array;
45+ modeF: Float32Array;
46+ modeZeta: Float32Array;
47+ modeZetaRad: Float32Array;
48+ modeFamily: Float32Array;
49+ nodeX: Float32Array;
50+ nodeY: Float32Array;
51+ nodeS: Float32Array;
52+ shapeUx: Float32Array;
53+ shapeUy: Float32Array;
54+ shapeUz: Float32Array;
55+ diag: {
56+ modesKept: number;
57+ modesIn: number;
58+ modesOut: number;
59+ strikeCount: number;
60+ meanV: number;
61+ maxV: number;
62+ meanTauUs: number;
63+ peakForceN: number;
64+ durationS: number;
65+ peakPa: number;
66+ peakSplDb: number;
67+ directPeakPa: number;
68+ listenerDistM: number;
69+ warnings: string[];
70+ solveMs: number;
71+ radiateMs: number;
72+ solveCached: boolean;
73+ radiateCached: boolean;
74+ oversample: number;
75+ };
76+}
77+
78+// --- shared helpers ---
79+
80+export function hashMFiles(files: MFile[]): string {
81+ let h = 0x811c9dc5;
82+ for (const f of files) {
83+ const s = f.name + ' ' + f.source;
84+ for (let i = 0; i < s.length; i++) {
85+ h ^= s.charCodeAt(i);
86+ h = Math.imul(h, 0x01000193);
87+ }
88+ }
89+ return (h >>> 0).toString(36);
90+}
91+
92+function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
93+ const v = cache.get(key);
94+ if (v !== undefined) {
95+ cache.delete(key);
96+ cache.set(key, v);
97+ }
98+ return v;
99+}
100+
101+function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
102+ cache.delete(key);
103+ cache.set(key, value);
104+ while (cache.size > max) {
105+ const oldest = cache.keys().next().value as string;
106+ cache.delete(oldest);
107+ }
108+}
109+
110+type Params = Record<string, number | boolean>;
111+
112+function pLiteral(params: Params): string {
113+ const lines: string[] = [];
114+ for (const [k, v] of Object.entries(params)) {
115+ if (typeof v === 'number') {
116+ if (!Number.isFinite(v)) throw new Error(`parameter ${k} is not finite`);
117+ lines.push(`P.${k} = ${v};`);
118+ } else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
119+ }
120+ return lines.join('\n');
121+}
122+
123+function keyFor(params: Params, keys: readonly string[]): string {
124+ return keys.map((k) => String(params[k] ?? 'd')).join('|');
125+}
126+
127+function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
128+ return executeCode(source, { initialVariableValues: inputs, displayResults: false }, mfiles, 'stage.m');
129+}
130+
131+function asTensor(v: RuntimeValue): RuntimeTensor {
132+ const t = v as RuntimeTensor;
133+ if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
134+ return t;
135+}
136+
137+function asScalar(v: RuntimeValue): number {
138+ if (typeof v === 'number') return v;
139+ const t = v as RuntimeTensor;
140+ if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
141+ if (typeof v === 'boolean') return v ? 1 : 0;
142+ throw new Error('expected scalar result from numbl stage');
143+}
144+
145+function toF32(v: RuntimeValue): Float32Array {
146+ const t = asTensor(v);
147+ const out = new Float32Array(t.data.length);
148+ for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
149+ return out;
150+}
151+
152+function knobValue(settings: Record<string, number>, key: string): number {
153+ const spec = TRIANGLE_KNOBS.find((k) => k.key === key)!;
154+ const raw = settings[key];
155+ return raw === undefined ? spec.default : Math.max(spec.min, Math.min(spec.max, raw));
156+}
157+
158+/** Velocity quantized to a few levels so identical strikes share a cached
159+ * render; 8 levels is inaudible next to the model's own dynamics. */
160+export function bucketVelocity(velocity01: number, levels = 8): number {
161+ return Math.ceil(Math.max(0.001, Math.min(1, velocity01)) * levels) / levels;
162+}
163+
164+// The basis fields handed from the solve stage into the radiate stage, in
165+// the order the reassembly script expects.
166+const BASIS_FIELDS = [
167+ 'f', 'omega', 'zeta', 'zetaRad', 'family',
168+ 'nodeX', 'nodeY', 'sfrac', 'totalLen',
169+ 'NUx', 'NUy', 'NUz',
170+ 'elemX', 'elemY', 'elemTanX', 'elemTanY', 'elemLen',
171+ 'EUpx', 'EUpy', 'EUpz',
172+];
173+
174+const SOLVE_KEYS = [
175+ 'side_length_m', 'gap_m', 'rod_diameter_m',
176+ 'density', 'youngs_modulus', 'poisson_ratio', 'stiffness_trim',
177+ 'n_elements', 'min_mode_hz', 'max_mode_hz',
178+ 'loss_factor', 'loss_slope', 'radiation_loss',
179+ 'speed_of_sound', 'air_density',
180+];
181+
182+const RADIATE_KEYS = [
183+ 'beater_radius_m', 'beater_mass_kg', 'beater_youngs', 'restitution',
184+ 'strike_velocity', 'strike_pos', 'strike_angle_deg',
185+ 'strike_count', 'strike_rate_hz', 'roll_pos2',
186+ 'timing_jitter_ms', 'velocity_spread', 'seed',
187+ 'direct_gain', 'listener_x_m', 'listener_y_m', 'listener_z_m',
188+ 'sample_rate', 'tail_s',
189+];
190+
191+interface SolveEntry {
192+ vars: Record<string, RuntimeValue>;
193+ modeF: Float32Array;
194+ modeZeta: Float32Array;
195+ modeZetaRad: Float32Array;
196+ modeFamily: Float32Array;
197+ nodeX: Float32Array;
198+ nodeY: Float32Array;
199+ nodeS: Float32Array;
200+ shapeUx: Float32Array;
201+ shapeUy: Float32Array;
202+ shapeUz: Float32Array;
203+ modesIn: number;
204+ modesOut: number;
205+}
206+
207+interface RadiateEntry {
208+ p: RuntimeValue;
209+ strikeT: Float32Array;
210+ strikeV: Float32Array;
211+ strikeS: Float32Array;
212+ strikeCount: number;
213+ meanV: number;
214+ maxV: number;
215+ meanTauUs: number;
216+ peakForceN: number;
217+ durationS: number;
218+ directPeakPa: number;
219+ listenerDistM: number;
220+ listenerInside: boolean;
221+ fsHi: number;
222+ oversample: number;
223+}
224+
225+// Caches, all keyed under the code hash so editing the .m sources
226+// invalidates everything computed from the previous version. The explorer
227+// caches are small; the note-solve cache holds one FEM eigensolve per
228+// pitch-scaled frame, which is what lets a track reuse the expensive stage
229+// across velocities and repeats.
230+const solveCache = new Map<string, SolveEntry>();
231+const radiateCache = new Map<string, RadiateEntry>();
232+const SOLVE_CACHE_MAX = 2;
233+const RADIATE_CACHE_MAX = 3;
234+interface NoteSolveEntry { vars: Record<string, RuntimeValue> }
235+const noteSolveCache = new Map<string, NoteSolveEntry>();
236+const NOTE_SOLVE_CACHE_MAX = 16;
237+
238+// Finished per-note audio: an in-memory LRU in front of an IndexedDB
239+// store, so renders survive page reloads (budgets and eviction live in
240+// noteCache.ts).
241+const noteCache = new NoteCache('triangle-instrument-notes');
242+
243+// --- the explorer path (full params, staged, with diagnostics) ---
244+
245+async function renderExplore(
246+ req: TrackRenderRequest,
247+ env: RenderEnv,
248+): Promise<{ audio: Float32Array; diag: TriangleExploreDiag }> {
249+ const params = req.params ?? {};
250+ const mfiles = env.mfiles;
251+ let quality = req.quality;
252+
253+ const codeKey = hashMFiles(mfiles);
254+ const solveKey = 'x|' + codeKey + '|' + keyFor(params, SOLVE_KEYS);
255+ const fullOversample = Number(params.oversample ?? 2);
256+ const radiateParamsKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);
257+ if (quality === 'preview' && radiateCache.has(radiateParamsKey + `|os=${fullOversample}`)) {
258+ quality = 'full'; // full-rate radiation already cached
259+ }
260+ const oversample = quality === 'preview' ? Math.min(1, fullOversample) : fullOversample;
261+ const sampleRate = Number(params.sample_rate ?? req.sampleRate);
262+ const fsHi = sampleRate * oversample;
263+
264+ const base = `P = triangle_defaults();\n${pLiteral(params)}`;
265+
266+ // ---- solve stage: FEM assembly and eigensolve ----
267+ let solveMs = 0;
268+ let solveCached = true;
269+ let solveEntry = lruGet(solveCache, solveKey);
270+ if (!solveEntry) {
271+ solveCached = false;
272+ env.post(0.15, 'solve');
273+ await env.yieldToHost();
274+ env.checkCancel();
275+ const t0 = performance.now();
276+ const r = runStage(mfiles, `
277+${base}
278+basis = triangle_basis(P);
279+${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
280+modesIn = sum(basis.family == 0);
281+modesOut = sum(basis.family == 1);
282+`);
283+ solveMs = performance.now() - t0;
284+
285+ const vars: Record<string, RuntimeValue> = {};
286+ for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];
287+ solveEntry = {
288+ vars,
289+ modeF: toF32(r.variableValues.b_f),
290+ modeZeta: toF32(r.variableValues.b_zeta),
291+ modeZetaRad: toF32(r.variableValues.b_zetaRad),
292+ modeFamily: toF32(r.variableValues.b_family),
293+ nodeX: toF32(r.variableValues.b_nodeX),
294+ nodeY: toF32(r.variableValues.b_nodeY),
295+ nodeS: toF32(r.variableValues.b_sfrac),
296+ shapeUx: toF32(r.variableValues.b_NUx),
297+ shapeUy: toF32(r.variableValues.b_NUy),
298+ shapeUz: toF32(r.variableValues.b_NUz),
299+ modesIn: asScalar(r.variableValues.modesIn),
300+ modesOut: asScalar(r.variableValues.modesOut),
301+ };
302+ lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
303+ }
304+
305+ // ---- radiate stage: strike process + modal synthesis ----
306+ const radiateKey = radiateParamsKey + `|os=${oversample}`;
307+ let radiateMs = 0;
308+ let radiateCached = true;
309+ let radEntry = lruGet(radiateCache, radiateKey);
310+ if (!radEntry) {
311+ radiateCached = false;
312+ env.post(0.55, 'radiate');
313+ await env.yieldToHost();
314+ env.checkCancel();
315+ const t0 = performance.now();
316+ const r = runStage(mfiles, `
317+${base}
318+${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}
319+basis.n = length(basis.f);
320+[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P);
321+[p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, ${fsHi});
322+strikeCount = sinfo.count;
323+meanV = sinfo.meanV;
324+maxV = sinfo.maxV;
325+meanTauUs = sinfo.meanTauUs;
326+peakForceN = sinfo.peakForceN;
327+durationS = sinfo.durationS;
328+directPeakPa = radinfo.directPeakPa;
329+listenerDistM = radinfo.listenerDistM;
330+listenerInside = radinfo.listenerInside;
331+`, solveEntry.vars);
332+ radiateMs = performance.now() - t0;
333+ radEntry = {
334+ p: r.variableValues.p,
335+ strikeT: toF32(r.variableValues.tstr),
336+ strikeV: toF32(r.variableValues.vstr),
337+ strikeS: toF32(r.variableValues.sstr),
338+ strikeCount: asScalar(r.variableValues.strikeCount),
339+ meanV: asScalar(r.variableValues.meanV),
340+ maxV: asScalar(r.variableValues.maxV),
341+ meanTauUs: asScalar(r.variableValues.meanTauUs),
342+ peakForceN: asScalar(r.variableValues.peakForceN),
343+ durationS: asScalar(r.variableValues.durationS),
344+ directPeakPa: asScalar(r.variableValues.directPeakPa),
345+ listenerDistM: asScalar(r.variableValues.listenerDistM),
346+ listenerInside: asScalar(r.variableValues.listenerInside) !== 0,
347+ fsHi,
348+ oversample,
349+ };
350+ lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
351+ }
352+
353+ // ---- output stage ----
354+ env.post(0.9, 'output');
355+ await env.yieldToHost();
356+ env.checkCancel();
357+ const rOut = runStage(mfiles, `
358+${base}
359+[y, peakPa] = triangle_output(P, p, ${radEntry.oversample});
360+`, { p: radEntry.p });
361+
362+ const y32 = toF32(rOut.variableValues.y);
363+ const peakPa = asScalar(rOut.variableValues.peakPa);
364+
365+ const warnings: string[] = [];
366+ if (radEntry.listenerInside) {
367+ warnings.push('the listener is within 2 cm of the rod; the dipole field is evaluated at a clamped minimum distance');
368+ }
369+ if (!params.normalize) {
370+ let peak = 0;
371+ for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
372+ if (peak > 1.0) {
373+ warnings.push(`clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${Number(params.full_scale_pa ?? 0.2).toFixed(2)}`);
374+ }
375+ }
376+ const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];
377+ const maxModeHz = Number(params.max_mode_hz ?? 0);
378+ if (maxModeHz > 0 && topMode < 0.4 * maxModeHz) {
379+ warnings.push(`the highest retained mode is only ${topMode.toFixed(0)} Hz; the shimmer above it falls entirely on the direct click`);
380+ }
381+
382+ const diag: TriangleExploreDiag = {
383+ quality,
384+ strikeT: radEntry.strikeT.slice(),
385+ strikeV: radEntry.strikeV.slice(),
386+ strikeS: radEntry.strikeS.slice(),
387+ modeF: solveEntry.modeF.slice(),
388+ modeZeta: solveEntry.modeZeta.slice(),
389+ modeZetaRad: solveEntry.modeZetaRad.slice(),
390+ modeFamily: solveEntry.modeFamily.slice(),
391+ nodeX: solveEntry.nodeX.slice(),
392+ nodeY: solveEntry.nodeY.slice(),
393+ nodeS: solveEntry.nodeS.slice(),
394+ shapeUx: solveEntry.shapeUx.slice(),
395+ shapeUy: solveEntry.shapeUy.slice(),
396+ shapeUz: solveEntry.shapeUz.slice(),
397+ diag: {
398+ modesKept: solveEntry.modeF.length,
399+ modesIn: solveEntry.modesIn,
400+ modesOut: solveEntry.modesOut,
401+ strikeCount: radEntry.strikeCount,
402+ meanV: radEntry.meanV,
403+ maxV: radEntry.maxV,
404+ meanTauUs: radEntry.meanTauUs,
405+ peakForceN: radEntry.peakForceN,
406+ durationS: radEntry.durationS,
407+ peakPa,
408+ peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
409+ directPeakPa: radEntry.directPeakPa,
410+ listenerDistM: radEntry.listenerDistM,
411+ warnings,
412+ solveMs,
413+ radiateMs,
414+ solveCached,
415+ radiateCached,
416+ oversample: radEntry.oversample,
417+ },
418+ };
419+ return { audio: y32, diag };
420+}
421+
422+// --- the track path (midiboot): per-strike rendering with caching ---
423+
424+function notePLines(note: { pitch: number; velocity: number }, req: TrackRenderRequest): string {
425+ // The note-to-parameter mapping: MIDI pitch scales the rod frame, side
426+ // length, gap, and rod diameter by s = 2^-((pitch - 84)/12), which
427+ // shifts the bending-mode spectrum by the semitone ratio exactly (rod
428+ // modes go as diameter / length^2). The beater stays fixed. Pitch 84 is
429+ // the default 6-inch triangle, playable 72..96; the spectrum is
430+ // inharmonic, so the note number is a handle on register, not a tuned
431+ // fundamental. Velocity drives the beater speed. Fixed: a single
432+ // strike, a fixed seed, normalize off so velocity survives in pascals.
433+ const settings = req.settings ?? {};
434+ const s = Math.pow(2, -(Math.max(72, Math.min(96, Math.round(note.pitch))) - 84) / 12);
435+ const lines = [
436+ 'P = triangle_defaults();',
437+ `P.sample_rate = ${req.sampleRate};`,
438+ 'P.strike_count = 1;',
439+ 'P.seed = 1;',
440+ 'P.pre_delay_ms = 0;',
441+ 'P.normalize = 0;',
442+ ];
443+ if (req.quality === 'preview') lines.push('P.oversample = 1;');
444+ lines.push(
445+ `P.side_length_m = P.side_length_m * ${s};`,
446+ `P.gap_m = P.gap_m * ${s};`,
447+ `P.rod_diameter_m = P.rod_diameter_m * ${s};`,
448+ `P.strike_velocity = ${0.1 + (6 - 0.1) * note.velocity};`,
449+ `P.strike_pos = ${knobValue(settings, 'strike_pos')};`,
450+ `P.strike_angle_deg = ${knobValue(settings, 'strike_angle_deg')};`,
451+ `P.tail_s = ${knobValue(settings, 'tail_s')};`,
452+ );
453+ return lines.join('\n');
454+}
455+
456+async function renderNoteAudio(
457+ note: { pitch: number; velocity: number },
458+ req: TrackRenderRequest,
459+ env: RenderEnv,
460+ codeKey: string,
461+ progress: (frac: number, stage: string) => void,
462+): Promise<Float32Array> {
463+ const pLines = notePLines(note, req);
464+
465+ // solve: one FEM eigensolve per pitch-scaled frame, shared across
466+ // velocities and knob values (the knobs are all radiate-tier
467+ // parameters).
468+ const pitch = Math.max(72, Math.min(96, Math.round(note.pitch)));
469+ const solveKey = `n|${codeKey}|${pitch}`;
470+ let entry = lruGet(noteSolveCache, solveKey);
471+ if (!entry) {
472+ progress(0.1, 'solve');
473+ await env.yieldToHost();
474+ env.checkCancel();
475+ const r = runStage(env.mfiles, `
476+${pLines}
477+basis = triangle_basis(P);
478+${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
479+`);
480+ const vars: Record<string, RuntimeValue> = {};
481+ for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];
482+ entry = { vars };
483+ lruSet(noteSolveCache, solveKey, entry, NOTE_SOLVE_CACHE_MAX);
484+ }
485+
486+ progress(0.5, 'radiate');
487+ await env.yieldToHost();
488+ env.checkCancel();
489+ const r = runStage(env.mfiles, `
490+${pLines}
491+${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}
492+basis.n = length(basis.f);
493+[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P);
494+fsHi = P.sample_rate * P.oversample;
495+[p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, fsHi);
496+[y, peakPa] = triangle_output(P, p, P.oversample);
497+`, entry.vars);
498+ return toF32(r.variableValues.y);
499+}
500+
501+async function renderTrackNotes(
502+ req: TrackRenderRequest,
503+ env: RenderEnv,
504+): Promise<{ audio: Float32Array }> {
505+ if (req.notes.length === 0) return { audio: new Float32Array(req.lengthSamples) };
506+ const codeKey = hashMFiles(env.mfiles);
507+ const settings = req.settings ?? {};
508+ const knobsKey = TRIANGLE_KNOBS.map((k) => `${k.key}=${knobValue(settings, k.key)}`).join(',');
509+
510+ // Group notes by everything that affects their sound. A strike ignores
511+ // note duration (ring-out is the tail_s knob), so a repeated pattern
512+ // costs one render per distinct pitch and velocity level.
513+ const groups = new Map<string, { pitch: number; velocity: number; onsets: number[] }>();
514+ for (const n of req.notes) {
515+ const pitch = Math.round(n.pitch);
516+ const vel = bucketVelocity(n.velocity);
517+ const key = `${codeKey}|${req.quality}|${req.sampleRate}|${knobsKey}|${pitch}|${vel}`;
518+ let g = groups.get(key);
519+ if (!g) groups.set(key, (g = { pitch, velocity: vel, onsets: [] }));
520+ g.onsets.push(n.timeSec);
521+ }
522+
523+ const placed: { buf: Float32Array; onsets: number[] }[] = [];
524+ // Resolve cache hits first (memory, then IndexedDB), so the progress
525+ // denominator counts only the sounds that actually need rendering.
526+ const resolved = new Map<string, Float32Array>();
527+ for (const key of groups.keys()) {
528+ const hit = await noteCache.get(key);
529+ if (hit) resolved.set(key, hit);
530+ }
531+ const misses = groups.size - resolved.size;
532+ let done = 0;
533+ for (const [key, g] of groups) {
534+ let buf = resolved.get(key);
535+ if (!buf) {
536+ const i = done;
537+ buf = await renderNoteAudio(g, req, env, codeKey, (frac, stage) =>
538+ env.post((i + frac) / Math.max(1, misses), `note ${i + 1}/${misses}: ${stage}`));
539+ noteCache.put(key, buf);
540+ done++;
541+ }
542+ placed.push({ buf, onsets: g.onsets });
543+ }
544+
545+ // lengthSamples is a minimum: size the output to the longest tail so a
546+ // ring that outlives the last note is returned in full, and place every
547+ // copy.
548+ let length = req.lengthSamples;
549+ for (const p of placed) {
550+ for (const t of p.onsets) length = Math.max(length, Math.round(t * req.sampleRate) + p.buf.length);
551+ }
552+ const out = new Float32Array(length);
553+ for (const p of placed) {
554+ for (const t of p.onsets) {
555+ const off = Math.round(t * req.sampleRate);
556+ const end = Math.min(out.length, off + p.buf.length);
557+ for (let j = Math.max(0, off); j < end; j++) out[j] += p.buf[j - off];
558+ }
559+ }
560+ return { audio: out };
561+}
562+
563+/** The single entry point: explorer path when diagnostics are requested,
564+ * the per-strike track path otherwise. */
565+export async function renderRequest(
566+ variant: TriangleVariant,
567+ req: TrackRenderRequest,
568+ env: RenderEnv,
569+): Promise<{ audio: Float32Array; diag?: TriangleExploreDiag }> {
570+ void variant;
571+ if (req.wantDiagnostics) return renderExplore(req, env);
572+ return renderTrackNotes(req, env);
573+}
instrument/index.tsadded+112−0View file
@@ -0,0 +1,112 @@
1+// The triangle behind the generic track-render instrument interface, for
2+// direct import by midiboot and by the triangle web app. The MATLAB model
3+// in ../matlab is executed by numbl in a lazily spawned worker that this
4+// module owns; the worker is kept alive across renders so its caches (the
5+// per-pitch eigensolves and finished strike renders) persist. Aborting a
6+// request is cooperative, checked between render stages; if the worker
7+// does not respond within a grace period it is torn down (and the caches
8+// with it) and respawned.
9+import type { Instrument, TrackRenderRequest, TrackRenderResult } from './types.ts';
10+import { TRIANGLE_KNOBS, type TriangleVariant } from './knobs.ts';
11+import type { WorkerCancelMsg, WorkerRenderMsg, WorkerResponse } from './worker.ts';
12+
13+const CANCEL_GRACE_MS = 5000;
14+
15+interface Pending {
16+ resolve: (r: TrackRenderResult) => void;
17+ reject: (e: Error) => void;
18+ onProgress?: (p: { frac: number; stage: string }) => void;
19+}
20+
21+let worker: Worker | null = null;
22+let nextId = 1;
23+const pending = new Map<number, Pending>();
24+
25+function ensureWorker(): Worker {
26+ if (worker) return worker;
27+ worker = new Worker(new URL('./worker.ts', import.meta.url), { type: 'module' });
28+ worker.onmessage = (ev: MessageEvent<WorkerResponse>) => {
29+ const msg = ev.data;
30+ const p = pending.get(msg.id);
31+ if (!p) return;
32+ if (msg.kind === 'progress') {
33+ p.onProgress?.({ frac: msg.frac, stage: msg.stage });
34+ return;
35+ }
36+ pending.delete(msg.id);
37+ if (msg.kind === 'result') p.resolve({ audio: msg.audio, diag: msg.diag });
38+ else if (msg.kind === 'cancelled') p.reject(new DOMException('render cancelled', 'AbortError'));
39+ else p.reject(Object.assign(new Error(msg.message), { file: msg.file, line: msg.line }));
40+ };
41+ worker.onerror = (ev) => {
42+ const err = new Error(`triangle instrument worker failed: ${ev.message}`);
43+ for (const p of pending.values()) p.reject(err);
44+ pending.clear();
45+ };
46+ return worker;
47+}
48+
49+function teardownWorker(reason: Error): void {
50+ if (worker) worker.terminate();
51+ worker = null;
52+ for (const p of pending.values()) p.reject(reason);
53+ pending.clear();
54+}
55+
56+function renderOnWorker(variant: TriangleVariant, req: TrackRenderRequest): Promise<TrackRenderResult> {
57+ if (req.signal?.aborted) return Promise.reject(new DOMException('render cancelled', 'AbortError'));
58+ const w = ensureWorker();
59+ const id = nextId++;
60+ const msg: WorkerRenderMsg = {
61+ kind: 'render',
62+ id,
63+ variant,
64+ notes: req.notes.map((n) => ({ timeSec: n.timeSec, durationSec: n.durationSec, pitch: n.pitch, velocity: n.velocity })),
65+ sampleRate: req.sampleRate,
66+ lengthSamples: req.lengthSamples,
67+ quality: req.quality,
68+ settings: req.settings ?? {},
69+ params: req.params,
70+ mfiles: req.mfiles,
71+ wantDiagnostics: req.wantDiagnostics,
72+ };
73+ return new Promise<TrackRenderResult>((resolve, reject) => {
74+ let graceTimer: ReturnType<typeof setTimeout> | null = null;
75+ const onAbort = () => {
76+ const cancel: WorkerCancelMsg = { kind: 'cancel', id };
77+ worker?.postMessage(cancel);
78+ // Escalation: a stage that never yields (or a hung worker) is torn
79+ // down after the grace period, losing the caches.
80+ graceTimer = setTimeout(() => {
81+ if (pending.has(id)) teardownWorker(new DOMException('render cancelled', 'AbortError'));
82+ }, CANCEL_GRACE_MS);
83+ };
84+ req.signal?.addEventListener('abort', onAbort, { once: true });
85+ pending.set(id, {
86+ resolve: (r) => {
87+ req.signal?.removeEventListener('abort', onAbort);
88+ if (graceTimer !== null) clearTimeout(graceTimer);
89+ resolve(r);
90+ },
91+ reject: (e) => {
92+ req.signal?.removeEventListener('abort', onAbort);
93+ if (graceTimer !== null) clearTimeout(graceTimer);
94+ reject(e);
95+ },
96+ onProgress: req.onProgress,
97+ });
98+ w.postMessage(msg);
99+ });
100+}
101+
102+class TriangleInstrument implements Instrument {
103+ readonly id = 'triangle';
104+ readonly name = 'Triangle';
105+ readonly knobs = TRIANGLE_KNOBS;
106+ renderTrack(req: TrackRenderRequest): Promise<TrackRenderResult> {
107+ return renderOnWorker('triangle', req);
108+ }
109+}
110+
111+export const triangleInstrument: Instrument = new TriangleInstrument();
112+export const triangleInstruments: Instrument[] = [triangleInstrument];
instrument/knobs.tsadded+11−0View file
@@ -0,0 +1,11 @@
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 TRIANGLE_KNOBS: KnobSpec[] = [
6+ { key: 'strike_pos', label: 'strike position', min: 0.02, max: 0.98, default: 0.17 },
7+ { key: 'strike_angle_deg', label: 'strike angle (deg)', min: 0, max: 90, default: 20 },
8+ { key: 'tail_s', label: 'ring-out (s)', min: 1, max: 8, default: 5 },
9+];
10+
11+export type TriangleVariant = 'triangle';
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+2301−0View file
This diff is 2,306 lines long and is not shown.
instrument/package.jsonadded+16−0View file
@@ -0,0 +1,16 @@
1+{
2+ "name": "triangle-instrument",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "description": "The triangle model 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+ "tsx": "^4.19.0"
15+ }
16+}
instrument/test.tsadded+122−0View file
@@ -0,0 +1,122 @@
1+// Node-side test of the triangle track-render instrument: executes the
2+// actual MATLAB model through numbl (no browser, no worker) and checks the
3+// audio. Run with: npx tsx test.ts
4+import { readdirSync, readFileSync } from 'node:fs';
5+import { join, dirname } from 'node:path';
6+import { fileURLToPath } from 'node:url';
7+import { renderRequest, type RenderEnv } from './core.ts';
8+import type { NoteSpec, TrackRenderRequest } from './types.ts';
9+
10+const matlabDir = join(dirname(fileURLToPath(import.meta.url)), '..', 'matlab');
11+const mfiles = readdirSync(matlabDir)
12+ .filter((f) => f.endsWith('.m'))
13+ .sort()
14+ .map((f) => ({ name: f, source: readFileSync(join(matlabDir, f), 'utf8') }));
15+
16+const env: RenderEnv = {
17+ mfiles,
18+ post: () => {},
19+ checkCancel: () => {},
20+ yieldToHost: async () => {},
21+};
22+
23+let failures = 0;
24+function check(label: string, ok: boolean, detail = ''): void {
25+ console.log(`${ok ? 'ok ' : 'FAIL'} ${label}${detail ? ` (${detail})` : ''}`);
26+ if (!ok) failures++;
27+}
28+
29+function stats(y: Float32Array) {
30+ let peak = 0;
31+ let finite = true;
32+ for (let i = 0; i < y.length; i++) {
33+ if (!Number.isFinite(y[i])) finite = false;
34+ peak = Math.max(peak, Math.abs(y[i]));
35+ }
36+ return { peak, finite };
37+}
38+
39+const SR = 22050;
40+
41+function req(notes: NoteSpec[], overrides: Partial<TrackRenderRequest> = {}): TrackRenderRequest {
42+ const end = notes.reduce((m, n) => Math.max(m, n.timeSec + n.durationSec), 0);
43+ return {
44+ notes,
45+ sampleRate: SR,
46+ lengthSamples: Math.ceil((end + 1.5) * SR),
47+ quality: 'preview',
48+ settings: { tail_s: 1 },
49+ ...overrides,
50+ };
51+}
52+
53+const note = (timeSec: number, pitch: number, velocity: number): NoteSpec =>
54+ ({ timeSec, durationSec: 0.1, pitch, velocity });
55+
56+const t0 = performance.now();
57+
58+// One strike: finite, audible, decaying.
59+const a = await renderRequest('triangle', req([note(0, 84, 0.7)]), env);
60+{
61+ const s = stats(a.audio);
62+ check('single strike is finite', s.finite);
63+ check('single strike is audible', s.peak > 1e-5, `peak ${s.peak.toExponential(2)}`);
64+ const n = a.audio.length;
65+ const head = stats(a.audio.slice(0, Math.floor(n / 3))).peak;
66+ const tail = stats(a.audio.slice(Math.floor((2 * n) / 3))).peak;
67+ check('single strike decays', tail < head * 0.8, `head ${head.toExponential(2)} tail ${tail.toExponential(2)}`);
68+}
69+
70+// The note cache.
71+{
72+ const t = performance.now();
73+ await renderRequest('triangle', req([note(0, 84, 0.7)]), env);
74+ const ms = performance.now() - t;
75+ check('repeat render hits the note cache', ms < 250, `${ms.toFixed(0)} ms`);
76+}
77+
78+// Velocity changes the sound (normalize is off, so amplitude scales); the
79+// second velocity reuses the pitch's cached eigensolve, so it only re-runs
80+// the radiate stage.
81+const soft = await (async () => {
82+ const t = performance.now();
83+ const r = await renderRequest('triangle', req([note(0, 84, 0.2)]), env);
84+ return { r, ms: performance.now() - t };
85+})();
86+{
87+ const pa = stats(a.audio).peak;
88+ const ps = stats(soft.r.audio).peak;
89+ check('velocity changes the sound', ps < pa, `soft ${ps.toExponential(2)} < hard ${pa.toExponential(2)}`);
90+}
91+
92+// Pitch and knobs change the sound.
93+const high = await renderRequest('triangle', req([note(0, 89, 0.7)]), env);
94+{
95+ let diff = high.audio.length !== a.audio.length;
96+ if (!diff) for (let i = 0; i < a.audio.length; i++) if (Math.abs(a.audio[i] - high.audio[i]) > 1e-9) { diff = true; break; }
97+ check('pitch changes the sound', diff);
98+}
99+const angled = await renderRequest('triangle', req([note(0, 84, 0.7)], { settings: { tail_s: 1, strike_angle_deg: 80 } }), env);
100+{
101+ let diff = false;
102+ const n = Math.min(a.audio.length, angled.audio.length);
103+ for (let i = 0; i < n; i++) if (Math.abs(a.audio[i] - angled.audio[i]) > 1e-9) { diff = true; break; }
104+ check('strike_angle_deg knob changes the sound', diff);
105+}
106+
107+// A pattern: repeats from the cache, distinct sounds at their onsets, and
108+// the buffer covers at least the requested length.
109+const tr = await renderRequest('triangle', req([note(0, 84, 0.7), note(0.6, 84, 0.7), note(1.2, 89, 0.7)]), env);
110+{
111+ const s = stats(tr.audio);
112+ check('track render is finite', s.finite);
113+ const at = (sec: number) => stats(tr.audio.slice(Math.floor(sec * SR), Math.floor((sec + 0.3) * SR))).peak;
114+ check('energy at first onset', at(0.0) > 1e-5);
115+ check('energy at second onset', at(0.6) > 1e-5);
116+ check('energy at third onset', at(1.2) > 1e-5);
117+ const asked = req([note(0, 84, 0.7), note(0.6, 84, 0.7), note(1.2, 89, 0.7)]).lengthSamples;
118+ check('buffer covers the requested length', tr.audio.length >= asked, `${tr.audio.length} >= ${asked}`);
119+}
120+
121+console.log(`total time: ${((performance.now() - t0) / 1000).toFixed(1)} s`);
122+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+120−0View file
@@ -0,0 +1,120 @@
1+// The render worker: numbl runs here so a seconds-long eigensolve never
2+// blocks the UI. Requests are processed one at a time; the core yields to
3+// the macrotask queue between stages, which is when queued cancel messages
4+// are delivered, so cancellation is cooperative and the caches survive it.
5+import { renderRequest, type TriangleVariant } from './core.ts';
6+import { DEFAULT_MFILES } from './mfiles.ts';
7+import type { MFile, NoteSpec } from './types.ts';
8+
9+export interface WorkerRenderMsg {
10+ kind: 'render';
11+ id: number;
12+ variant: TriangleVariant;
13+ notes: NoteSpec[];
14+ sampleRate: number;
15+ lengthSamples: number;
16+ quality: 'preview' | 'full';
17+ settings?: Record<string, number>;
18+ params?: Record<string, number | boolean>;
19+ mfiles?: MFile[];
20+ wantDiagnostics?: boolean;
21+}
22+
23+export interface WorkerCancelMsg {
24+ kind: 'cancel';
25+ id: number;
26+}
27+
28+export type WorkerResponse =
29+ | { kind: 'progress'; id: number; frac: number; stage: string }
30+ | { kind: 'result'; id: number; audio: Float32Array; diag?: unknown }
31+ | { kind: 'cancelled'; id: number }
32+ | { kind: 'error'; id: number; message: string; file?: string; line?: number };
33+
34+class Cancelled extends Error {}
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+const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
44+
45+async function run(msg: WorkerRenderMsg): Promise<void> {
46+ // numbl reports a file it could not parse through console.warn and then
47+ // skips it, so the failure surfaces at the call site as "Undefined
48+ // function". Capture the warning to report the actual syntax error.
49+ const warnings: string[] = [];
50+ const origWarn = console.warn;
51+ console.warn = (...args: unknown[]) => { warnings.push(args.map(String).join(' ')); };
52+ try {
53+ const result = await renderRequest(msg.variant, {
54+ notes: msg.notes,
55+ sampleRate: msg.sampleRate,
56+ lengthSamples: msg.lengthSamples,
57+ quality: msg.quality,
58+ settings: msg.settings,
59+ params: msg.params,
60+ wantDiagnostics: msg.wantDiagnostics,
61+ }, {
62+ mfiles: msg.mfiles ?? DEFAULT_MFILES,
63+ post: (frac, stage) => post({ kind: 'progress', id: msg.id, frac, stage }),
64+ checkCancel: () => { if (cancelled.has(msg.id)) throw new Cancelled(); },
65+ yieldToHost: () => new Promise((res) => setTimeout(res, 0)),
66+ });
67+ const transfer: Transferable[] = [result.audio.buffer];
68+ if (result.diag) {
69+ for (const v of Object.values(result.diag as unknown as Record<string, unknown>)) {
70+ if (v instanceof Float32Array && v.buffer !== result.audio.buffer) transfer.push(v.buffer);
71+ }
72+ }
73+ post({ kind: 'result', id: msg.id, audio: result.audio, diag: result.diag }, transfer);
74+ } catch (err) {
75+ if (err instanceof Cancelled) {
76+ post({ kind: 'cancelled', id: msg.id });
77+ } else {
78+ const e = err as { message?: string; file?: string; line?: number };
79+ let m: WorkerResponse = {
80+ kind: 'error',
81+ id: msg.id,
82+ message: e?.message ? String(e.message) : String(err),
83+ file: typeof e?.file === 'string' ? e.file : undefined,
84+ line: typeof e?.line === 'number' ? e.line : undefined,
85+ };
86+ for (const w of warnings) {
87+ const hit = SKIP_WARNING.exec(w);
88+ if (hit) {
89+ m = { kind: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
90+ break;
91+ }
92+ }
93+ post(m);
94+ }
95+ } finally {
96+ console.warn = origWarn;
97+ cancelled.delete(msg.id);
98+ }
99+}
100+
101+async function pump(): Promise<void> {
102+ if (running) return;
103+ running = true;
104+ while (queue.length > 0) {
105+ const msg = queue.shift()!;
106+ if (cancelled.has(msg.id)) { post({ kind: 'cancelled', id: msg.id }); cancelled.delete(msg.id); continue; }
107+ await run(msg);
108+ }
109+ running = false;
110+}
111+
112+self.onmessage = (ev: MessageEvent<WorkerRenderMsg | WorkerCancelMsg>) => {
113+ const msg = ev.data;
114+ if (msg.kind === 'render') {
115+ queue.push(msg);
116+ void pump();
117+ } else if (msg.kind === 'cancel') {
118+ cancelled.add(msg.id);
119+ }
120+};
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
@@ -36,7 +36,7 @@ The New seed button redraws the roll timing and velocities without touching any
3636 - `../matlab/` — the algorithm: `triangle_geometry.m`, `triangle_fem.m`, `triangle_basis.m`, `triangle_strikes.m`, `triangle_radiate.m`, `triangle_output.m`, `triangle_rng.m`, plus `triangle_render.m` / `triangle_demo.m` / `triangle_validate.m` for CLI use
3737 - `src/params.ts` — parameter schema, defaults, tier assignment, persistence
3838 - `src/mfiles.ts` — the `.m` sources, bundled at build time and hashed for cache keying
39-- `src/worker/renderWorker.ts` — off-thread numbl execution with tiered caching
39+- `../instrument/` — off-thread numbl execution with tiered caching, behind the shared track-render instrument interface (also used by midiboot)
4040 - `src/lib/outline.ts` — the open-triangle centerline, mirrored from `triangle_geometry.m` for the scene and the mode-shape plot
4141 - `src/components/` — 3D scene (three.js), parameter panel, code panel and editor, plots (strikes, waveform, spectrum, spectrogram, mode map, mode shapes)
4242 - `src/selftest.ts` — Node harness for the MATLAB validation script
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 const modes = Array.from(result.modeF);
66 const shown = modes.slice(0, 10);
webapp/src/hooks/useRenderer.tsmodified+72−65View file
@@ -1,12 +1,24 @@
1-// Hook owning the render worker. Renders are explicit: nothing starts until
2-// Render preview or Render full is pressed, so changing a parameter or a .m
3-// source never kicks off work on its own. A request in flight can be
4-// cancelled, 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.
59
610 import { useCallback, useEffect, useRef, useState } from 'react';
711 import type { TriangleParams } from '../params';
812 import type { MFile } from '../mfiles';
9-import type { ErrorMessage, ProgressMessage, RenderResponse } from '../worker/renderWorker';
13+import { triangleInstrument } from '../../../instrument/index.ts';
14+import type { TriangleExploreDiag } 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 = TriangleExploreDiag & {
19+ y: Float32Array;
20+ sampleRate: number;
21+};
1022
1123 export interface RenderError {
1224 message: string;
@@ -15,7 +27,7 @@ export interface RenderError {
1527 }
1628
1729 export interface RenderState {
18- result: RenderResponse | null;
30+ result: RenderResult | null;
1931 error: RenderError | null;
2032 busy: boolean;
2133 stage: string;
@@ -24,9 +36,7 @@ export interface RenderState {
2436 }
2537
2638 export function useRenderer(params: TriangleParams, mfiles: MFile[]) {
27- const workerRef = useRef<Worker | null>(null);
28- const nextId = useRef(1);
29- const inFlight = useRef<number | null>(null);
39+ const abortRef = useRef<AbortController | null>(null);
3040 const startedAt = useRef(0);
3141 const [state, setState] = useState<RenderState>({
3242 result: null,
@@ -37,49 +47,7 @@ export function useRenderer(params: TriangleParams, mfiles: MFile[]) {
3747 elapsedS: 0,
3848 });
3949
40- const attachHandlers = useCallback((worker: Worker) => {
41- worker.onmessage = (ev: MessageEvent<RenderResponse | ProgressMessage | ErrorMessage>) => {
42- const msg = ev.data;
43- if (msg.type === 'progress') {
44- setState((s) => ({ ...s, stage: msg.stage, progress: msg.frac }));
45- return;
46- }
47- inFlight.current = null;
48- if (msg.type === 'error') {
49- // A failed render leaves nothing playable: the previous result is
50- // dropped, so Play and Download cannot hand back audio that no
51- // longer corresponds to the code or parameters on screen.
52- setState((s) => ({
53- ...s,
54- result: null,
55- error: { message: msg.message, file: msg.file, line: msg.line },
56- busy: false,
57- stage: '',
58- progress: 0,
59- }));
60- return;
61- }
62- setState((s) => ({ ...s, result: msg, error: null, busy: false, stage: '', progress: 0 }));
63- };
64- }, []);
65-
66- const makeWorker = useCallback(() => {
67- const worker = new Worker(new URL('../worker/renderWorker.ts', import.meta.url), {
68- type: 'module',
69- });
70- attachHandlers(worker);
71- workerRef.current = worker;
72- return worker;
73- }, [attachHandlers]);
74-
75- useEffect(() => {
76- const worker = makeWorker();
77- return () => {
78- worker.terminate();
79- workerRef.current = null;
80- inFlight.current = null;
81- };
82- }, [makeWorker]);
50+ useEffect(() => () => abortRef.current?.abort(), []);
8351
8452 // elapsed-time ticker while a render is in flight
8553 useEffect(() => {
@@ -93,27 +61,66 @@ export function useRenderer(params: TriangleParams, mfiles: MFile[]) {
9361
9462 const render = useCallback(
9563 (quality: 'preview' | 'full') => {
96- const worker = workerRef.current;
97- if (!worker || inFlight.current !== null) return;
98- const id = nextId.current++;
99- inFlight.current = id;
64+ if (abortRef.current) return; // a request is already in flight
65+ const ctrl = new AbortController();
66+ abortRef.current = ctrl;
10067 startedAt.current = performance.now();
10168 setState((s) => ({ ...s, busy: true, stage: 'starting', progress: 0, elapsedS: 0 }));
102- worker.postMessage({ type: 'render', id, params, mfiles, quality });
69+ const sampleRate = Number(params.sample_rate);
70+ void triangleInstrument.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 TriangleExploreDiag;
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+ });
103115 },
104116 [params, mfiles],
105117 );
106118
107- // Cancelling terminates the worker outright: the .m interpreter runs a tight
108- // loop that never yields, so there is nothing to poll a flag from. The
109- // replacement worker starts with empty stage caches.
110119 const cancel = useCallback(() => {
111- const worker = workerRef.current;
112- if (worker) worker.terminate();
113- inFlight.current = null;
114- makeWorker();
120+ abortRef.current?.abort();
121+ abortRef.current = null;
115122 setState((s) => ({ ...s, busy: false, stage: '', progress: 0, elapsedS: 0 }));
116- }, [makeWorker]);
123+ }, []);
117124
118125 const renderPreview = useCallback(() => render('preview'), [render]);
119126 const renderFull = useCallback(() => render('full'), [render]);
webapp/src/worker/renderWorker.tsdeleted+0−430View file
@@ -1,430 +0,0 @@
1-// Web Worker that renders the triangle strike by executing the MATLAB
2-// implementation in ../matlab through numbl. The stages run as separate
3-// executeCode calls so intermediates can be cached and cheap parameter
4-// changes skip the expensive ones:
5-// solve tier -> triangle_basis (FEM assembly + eigensolve; the expensive
6-// stage, cached on geometry/material/damping params)
7-// radiate tier -> triangle_strikes + triangle_radiate (cached on the
8-// strike process, the listener, and the synthesis rate)
9-// output tier -> triangle_output (always cheap, always re-run)
10-// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).
11-
12-import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
13-import { hashMFiles, type MFile } from '../mfiles';
14-import { SOLVE_KEYS, RADIATE_KEYS, type TriangleParams } from '../params';
15-
16-export interface RenderRequest {
17- type: 'render';
18- id: number;
19- params: TriangleParams;
20- quality: 'preview' | 'full';
21- mfiles: MFile[];
22-}
23-
24-export interface RenderResponse {
25- type: 'result';
26- id: number;
27- quality: 'preview' | 'full';
28- y: Float32Array; // final signal, [-1, 1]
29- sampleRate: number;
30- strikeT: Float32Array; // strike times, s
31- strikeV: Float32Array; // strike velocities, m/s
32- strikeS: Float32Array; // strike arclength fractions
33- modeF: Float32Array; // mode frequencies, Hz
34- modeZeta: Float32Array; // damping ratio of each mode
35- modeZetaRad: Float32Array; // radiation part of the damping
36- modeFamily: Float32Array; // 0 = in-plane, 1 = out-of-plane
37- nodeX: Float32Array; // rod centerline nodes, m
38- nodeY: Float32Array;
39- nodeS: Float32Array; // arclength fraction of each node
40- shapeUx: Float32Array; // node displacements per mode, nn x n col-major
41- shapeUy: Float32Array;
42- shapeUz: Float32Array;
43- diag: {
44- modesKept: number;
45- modesIn: number;
46- modesOut: number;
47- strikeCount: number;
48- meanV: number;
49- maxV: number;
50- meanTauUs: number;
51- peakForceN: number;
52- durationS: number;
53- peakPa: number;
54- peakSplDb: number;
55- directPeakPa: number;
56- listenerDistM: number;
57- warnings: string[];
58- solveMs: number; // 0 when the cached eigensolve was reused
59- radiateMs: number; // 0 when the cached radiation was reused
60- solveCached: boolean;
61- radiateCached: boolean;
62- oversample: number;
63- };
64-}
65-
66-export interface ProgressMessage {
67- type: 'progress';
68- id: number;
69- stage: 'solve' | 'radiate' | 'output';
70- frac: number;
71-}
72-
73-// Raised by a fault in the .m code (the user can edit it), so it is reported
74-// rather than left to stall the request.
75-export interface ErrorMessage {
76- type: 'error';
77- id: number;
78- message: string;
79- file?: string;
80- line?: number;
81-}
82-
83-// Basis fields that cross from the solve stage into the radiate stage as
84-// numbl RuntimeValues, plus their Float32 copies for the plots.
85-interface SolveEntry {
86- vars: Record<string, RuntimeValue>;
87- modeF: Float32Array;
88- modeZeta: Float32Array;
89- modeZetaRad: Float32Array;
90- modeFamily: Float32Array;
91- nodeX: Float32Array;
92- nodeY: Float32Array;
93- nodeS: Float32Array;
94- shapeUx: Float32Array;
95- shapeUy: Float32Array;
96- shapeUz: Float32Array;
97- modesIn: number;
98- modesOut: number;
99-}
100-
101-interface RadiateEntry {
102- p: RuntimeValue;
103- strikeT: Float32Array;
104- strikeV: Float32Array;
105- strikeS: Float32Array;
106- strikeCount: number;
107- meanV: number;
108- maxV: number;
109- meanTauUs: number;
110- peakForceN: number;
111- durationS: number;
112- directPeakPa: number;
113- listenerDistM: number;
114- listenerInside: boolean;
115- fsHi: number;
116- oversample: number;
117-}
118-
119-// Small LRU caches. The solve cache does not depend on the oversample rate,
120-// so a preview render right after a full one reuses the same eigensolve.
121-const solveCache = new Map<string, SolveEntry>();
122-const radiateCache = new Map<string, RadiateEntry>();
123-const SOLVE_CACHE_MAX = 2;
124-const RADIATE_CACHE_MAX = 3;
125-
126-function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
127- const v = cache.get(key);
128- if (v !== undefined) {
129- cache.delete(key);
130- cache.set(key, v);
131- }
132- return v;
133-}
134-
135-function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
136- cache.delete(key);
137- cache.set(key, value);
138- while (cache.size > max) {
139- const oldest = cache.keys().next().value as string;
140- cache.delete(oldest);
141- }
142-}
143-
144-function keyFor(params: TriangleParams, keys: readonly (keyof TriangleParams)[]): string {
145- return keys.map((k) => String(params[k])).join('|');
146-}
147-
148-// MATLAB source lines assigning every parameter onto a struct P.
149-function pLiteral(params: TriangleParams): string {
150- const lines: string[] = [];
151- for (const [k, v] of Object.entries(params)) {
152- if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
153- else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
154- }
155- return lines.join('\n');
156-}
157-
158-function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
159- return executeCode(
160- source,
161- { initialVariableValues: inputs, displayResults: false },
162- mfiles,
163- 'stage.m',
164- );
165-}
166-
167-function asTensor(v: RuntimeValue): RuntimeTensor {
168- const t = v as RuntimeTensor;
169- if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
170- return t;
171-}
172-
173-function asScalar(v: RuntimeValue): number {
174- if (typeof v === 'number') return v;
175- const t = v as RuntimeTensor;
176- if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
177- if (typeof v === 'boolean') return v ? 1 : 0;
178- throw new Error('expected scalar result from numbl stage');
179-}
180-
181-function toF32(v: RuntimeValue): Float32Array {
182- const t = asTensor(v);
183- const out = new Float32Array(t.data.length);
184- for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
185- return out;
186-}
187-
188-// The basis fields handed from the solve stage into the radiate stage, in the
189-// order the reassembly script expects.
190-const BASIS_FIELDS = [
191- 'f', 'omega', 'zeta', 'zetaRad', 'family',
192- 'nodeX', 'nodeY', 'sfrac', 'totalLen',
193- 'NUx', 'NUy', 'NUz',
194- 'elemX', 'elemY', 'elemTanX', 'elemTanY', 'elemLen',
195- 'EUpx', 'EUpy', 'EUpz',
196-];
197-
198-const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
199-
200-self.onmessage = (ev: MessageEvent<RenderRequest>) => {
201- const msg = ev.data;
202- if (msg.type !== 'render') return;
203-
204- // numbl reports a file it could not parse through console.warn and then
205- // skips it, so the failure surfaces at the call site as "Undefined
206- // function". Capture the warning to report the actual syntax error instead.
207- const warnings: string[] = [];
208- const origWarn = console.warn;
209- console.warn = (...args: unknown[]) => {
210- warnings.push(args.map(String).join(' '));
211- };
212-
213- try {
214- render(msg);
215- } catch (err) {
216- // Most likely a fault in the (user-editable) .m code. Report it so the
217- // request completes and the message reaches the code panel.
218- const e = err as { message?: string; file?: string; line?: number };
219- let m: ErrorMessage = {
220- type: 'error',
221- id: msg.id,
222- message: e?.message ? String(e.message) : String(err),
223- file: typeof e?.file === 'string' ? e.file : undefined,
224- line: typeof e?.line === 'number' ? e.line : undefined,
225- };
226- for (const w of warnings) {
227- const hit = SKIP_WARNING.exec(w);
228- if (hit) {
229- m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
230- break;
231- }
232- }
233- (self as unknown as Worker).postMessage(m);
234- } finally {
235- console.warn = origWarn;
236- }
237-};
238-
239-function render(msg: RenderRequest): void {
240- const { id, params, mfiles } = msg;
241- let quality = msg.quality;
242-
243- // The .m sources are part of the cache identity: editing the code must
244- // invalidate everything computed from the previous version.
245- const codeKey = hashMFiles(mfiles);
246- const solveKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);
247-
248- const fullOversample = params.oversample;
249- const radiateParamsKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);
250- if (quality === 'preview' && radiateCache.has(radiateParamsKey + `|os=${fullOversample}`)) {
251- quality = 'full'; // full-rate radiation already cached
252- }
253- const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;
254- const fsHi = params.sample_rate * oversample;
255-
256- const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
257-
258- // ---- solve stage: FEM assembly and eigensolve ----
259- let solveMs = 0;
260- let solveCached = true;
261- let solveEntry = lruGet(solveCache, solveKey);
262- if (!solveEntry) {
263- solveCached = false;
264- post({ type: 'progress', id, stage: 'solve', frac: 0.15 });
265- const t0 = performance.now();
266- const r = runStage(mfiles, `
267-${pLiteral(params)}
268-basis = triangle_basis(P);
269-${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
270-modesIn = sum(basis.family == 0);
271-modesOut = sum(basis.family == 1);
272-`);
273- solveMs = performance.now() - t0;
274-
275- const vars: Record<string, RuntimeValue> = {};
276- for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];
277- solveEntry = {
278- vars,
279- modeF: toF32(r.variableValues.b_f),
280- modeZeta: toF32(r.variableValues.b_zeta),
281- modeZetaRad: toF32(r.variableValues.b_zetaRad),
282- modeFamily: toF32(r.variableValues.b_family),
283- nodeX: toF32(r.variableValues.b_nodeX),
284- nodeY: toF32(r.variableValues.b_nodeY),
285- nodeS: toF32(r.variableValues.b_sfrac),
286- shapeUx: toF32(r.variableValues.b_NUx),
287- shapeUy: toF32(r.variableValues.b_NUy),
288- shapeUz: toF32(r.variableValues.b_NUz),
289- modesIn: asScalar(r.variableValues.modesIn),
290- modesOut: asScalar(r.variableValues.modesOut),
291- };
292- lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
293- }
294-
295- // ---- radiate stage: strike process + modal synthesis ----
296- const radiateKey = radiateParamsKey + `|os=${oversample}`;
297- let radiateMs = 0;
298- let radiateCached = true;
299- let radEntry = lruGet(radiateCache, radiateKey);
300- if (!radEntry) {
301- radiateCached = false;
302- post({ type: 'progress', id, stage: 'radiate', frac: 0.55 });
303- const t0 = performance.now();
304- const r = runStage(mfiles, `
305-${pLiteral(params)}
306-${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}
307-basis.n = length(basis.f);
308-[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P);
309-[p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, ${fsHi});
310-strikeCount = sinfo.count;
311-meanV = sinfo.meanV;
312-maxV = sinfo.maxV;
313-meanTauUs = sinfo.meanTauUs;
314-peakForceN = sinfo.peakForceN;
315-durationS = sinfo.durationS;
316-directPeakPa = radinfo.directPeakPa;
317-listenerDistM = radinfo.listenerDistM;
318-listenerInside = radinfo.listenerInside;
319-`, solveEntry.vars);
320- radiateMs = performance.now() - t0;
321- radEntry = {
322- p: r.variableValues.p,
323- strikeT: toF32(r.variableValues.tstr),
324- strikeV: toF32(r.variableValues.vstr),
325- strikeS: toF32(r.variableValues.sstr),
326- strikeCount: asScalar(r.variableValues.strikeCount),
327- meanV: asScalar(r.variableValues.meanV),
328- maxV: asScalar(r.variableValues.maxV),
329- meanTauUs: asScalar(r.variableValues.meanTauUs),
330- peakForceN: asScalar(r.variableValues.peakForceN),
331- durationS: asScalar(r.variableValues.durationS),
332- directPeakPa: asScalar(r.variableValues.directPeakPa),
333- listenerDistM: asScalar(r.variableValues.listenerDistM),
334- listenerInside: asScalar(r.variableValues.listenerInside) !== 0,
335- fsHi,
336- oversample,
337- };
338- lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
339- }
340-
341- // ---- output stage ----
342- post({ type: 'progress', id, stage: 'output', frac: 0.9 });
343- const rOut = runStage(mfiles, `
344-${pLiteral(params)}
345-[y, peakPa] = triangle_output(P, p, ${radEntry.oversample});
346-`, { p: radEntry.p });
347-
348- const y32 = toF32(rOut.variableValues.y);
349- const peakPa = asScalar(rOut.variableValues.peakPa);
350-
351- const warnings: string[] = [];
352- if (radEntry.listenerInside) {
353- warnings.push(
354- 'the listener is within 2 cm of the rod; the dipole field is evaluated at a clamped minimum distance',
355- );
356- }
357- if (!params.normalize) {
358- let peak = 0;
359- for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
360- if (peak > 1.0) {
361- warnings.push(
362- `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,
363- );
364- }
365- }
366- const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];
367- if (topMode < 0.4 * params.max_mode_hz) {
368- warnings.push(
369- `the highest retained mode is only ${topMode.toFixed(0)} Hz; the shimmer above it falls entirely on the direct click`,
370- );
371- }
372-
373- const resp: RenderResponse = {
374- type: 'result',
375- id,
376- quality,
377- y: y32,
378- sampleRate: params.sample_rate,
379- strikeT: radEntry.strikeT.slice(),
380- strikeV: radEntry.strikeV.slice(),
381- strikeS: radEntry.strikeS.slice(),
382- modeF: solveEntry.modeF.slice(),
383- modeZeta: solveEntry.modeZeta.slice(),
384- modeZetaRad: solveEntry.modeZetaRad.slice(),
385- modeFamily: solveEntry.modeFamily.slice(),
386- nodeX: solveEntry.nodeX.slice(),
387- nodeY: solveEntry.nodeY.slice(),
388- nodeS: solveEntry.nodeS.slice(),
389- shapeUx: solveEntry.shapeUx.slice(),
390- shapeUy: solveEntry.shapeUy.slice(),
391- shapeUz: solveEntry.shapeUz.slice(),
392- diag: {
393- modesKept: solveEntry.modeF.length,
394- modesIn: solveEntry.modesIn,
395- modesOut: solveEntry.modesOut,
396- strikeCount: radEntry.strikeCount,
397- meanV: radEntry.meanV,
398- maxV: radEntry.maxV,
399- meanTauUs: radEntry.meanTauUs,
400- peakForceN: radEntry.peakForceN,
401- durationS: radEntry.durationS,
402- peakPa,
403- peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
404- directPeakPa: radEntry.directPeakPa,
405- listenerDistM: radEntry.listenerDistM,
406- warnings,
407- solveMs,
408- radiateMs,
409- solveCached,
410- radiateCached,
411- oversample: radEntry.oversample,
412- },
413- };
414- (self as unknown as Worker).postMessage(resp, [
415- resp.y.buffer,
416- resp.strikeT.buffer,
417- resp.strikeV.buffer,
418- resp.strikeS.buffer,
419- resp.modeF.buffer,
420- resp.modeZeta.buffer,
421- resp.modeZetaRad.buffer,
422- resp.modeFamily.buffer,
423- resp.nodeX.buffer,
424- resp.nodeY.buffer,
425- resp.nodeS.buffer,
426- resp.shapeUx.buffer,
427- resp.shapeUy.buffer,
428- resp.shapeUz.buffer,
429- ]);
430-}
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,