10729f3WebGPU time loop for the cymbal2 shell model, with numbl setup and output stagesJeremy Magland 1// Setup and output stages of the cymbal2 render, run through numbl exactly as
2// the webapp runs them. The GPU path replaces only cymbal2_sim: everything the
3// time loop needs (grid operators, per-wavenumber inverses, mallet matrices,
4// stroke schedule) is computed here by the unmodified MATLAB code and pulled
5// out as tensors.
7import { readFileSync, readdirSync } from 'node:fs';
8import { fileURLToPath } from 'node:url';
9import path from 'node:path';
10import { executeCode, RTV } from 'numbl';
12const MATLAB_DIR = path.resolve(
13 path.dirname(fileURLToPath(import.meta.url)),
14 '../../../simulated_instruments/cymbal2/matlab',
15);
17export function loadMFiles() {
18 return readdirSync(MATLAB_DIR)
19 .filter((f) => f.endsWith('.m'))
20 .sort()
21 .map((name) => ({
22 name,
23 source: readFileSync(path.join(MATLAB_DIR, name), 'utf8'),
24 }));
25}
27// MATLAB lines overriding defaults from a plain params object.
28function pOverrides(params) {
29 const lines = [];
30 for (const [k, v] of Object.entries(params ?? {})) {
31 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
32 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
33 }
34 return lines.join('\n');
35}
37function stage(mfiles, source, inputs) {
38 return executeCode(
39 source,
40 { initialVariableValues: inputs, displayResults: false },
41 mfiles,
42 'stage.m',
43 );
44}
46function tensor(v, name) {
47 if (!v || v.kind !== 'tensor') throw new Error(`expected tensor for ${name}`);
48 return v; // {kind, shape, data: Float64Array} column-major
49}
51function scalar(v, name) {
52 if (typeof v === 'number') return v;
53 if (v && v.kind === 'tensor' && v.data.length === 1) return v.data[0];
54 throw new Error(`expected scalar for ${name}`);
55}
57function f32(v) {
58 const t = v.data;
59 const out = new Float32Array(t.length);
60 out.set(t);
61 return out;
62}
64// The stroke-schedule block is the one in cymbal2_sim.m, verbatim, so the GPU
65// render integrates exactly the schedule the CPU render would.
66const SCHEDULE_SRC = `
67K = max(1, round(P.strike_count));
68tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);
69vstr = zeros(1, K);
70state = P.seed * 2654435761;
71for i = 1:K
72 state = mod(state * 1664525 + 1013904223, 4294967296);
73 x = state / 4294967296;
74 vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));
75end
76gap = P.strike_gap_m;
77if gap <= 0
78 gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;
79end
80travel = gap ./ vstr + vstr / (2 * P.strike_accel);
81tlaunch = tstr - travel;
82off = -min(0, min(tlaunch));
83tlaunch = tlaunch + off;
84duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
85nsteps = floor(duration * P.sample_rate * P.oversample);
86`;
88export function runSetup(params) {
89 const mfiles = loadMFiles();
90 const src = `
91P = cymbal2_defaults();
92${pOverrides(params)}
93S = cymbal2_system(P);
94MB = cymbal2_mallet(P);
95G = S.G;
96o_D1 = G.D1; o_D2 = G.D2;
97o_F = G.F; o_Gi = G.Gi;
98o_Rd = G.Rd; o_R2d = G.R2d;
99o_AiB = S.AiB; o_AiC = S.AiC; o_Ai = S.Ai; o_Fi = S.Fi;
100o_ASiB = MB.ASiB; o_ASiC = MB.ASiC;
101o_vM = MB.vM; o_vR = MB.vR; o_vK = MB.vK; o_tipcol = MB.tipcol;
102${SCHEDULE_SRC}
103o_tlaunch = tlaunch; o_vstr = vstr;
104c_nrt = G.nrt; c_Nr = G.Nr; c_nth = G.nth; c_hr = G.hr; c_hth = G.hth;
105c_ht = S.ht; c_Eh = S.Eh; c_Rs = S.Rs; c_rhoh = S.rhoh;
106c_fsc = S.forceScale; c_nlr0 = S.eomRow0; c_nDropped = S.nDropped;
107c_ls = S.ls; c_js = S.js; c_lp = S.lp; c_jp = S.jp;
108c_strikeR = S.strikeRadiusM; c_pickupR = S.pickupRadiusM;
109c_tipi = MB.tipi; c_nfree = MB.nfree; c_ps = MB.ps;
110c_KM = P.head_stiffness; c_alM = P.head_alpha;
111c_aH = P.strike_accel; c_aHr = P.retract_accel;
112c_holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / S.ht));
113c_missDepth = gap + 0.06;
114c_gap = gap; c_K = K; c_nsteps = nsteps;
115c_simRate = P.sample_rate * P.oversample;
116c_thickness = P.thickness_m;
117`;
118 const r = stage(mfiles, src);
119 const v = r.variableValues;
120 const g = (n) => tensor(v[n], n);
121 const s = (n) => scalar(v[n], n);
123 const setup = {
124 mfiles,
125 ops: {
126 D1: f32(g('o_D1')), D2: f32(g('o_D2')),
127 F: f32(g('o_F')), Gi: f32(g('o_Gi')),
128 Rd: f32(g('o_Rd')), R2d: f32(g('o_R2d')),
129 AiB: f32(g('o_AiB')), AiC: f32(g('o_AiC')),
130 Ai: f32(g('o_Ai')), Fi: f32(g('o_Fi')),
131 ASiB: f32(g('o_ASiB')), ASiC: f32(g('o_ASiC')),
132 vM: f32(g('o_vM')), vR: f32(g('o_vR')),
133 vK: f32(g('o_vK')), tipcol: f32(g('o_tipcol')),
134 },
135 tlaunch: Array.from(g('o_tlaunch').data),
136 vstr: Array.from(g('o_vstr').data),
137 nrt: s('c_nrt'), Nr: s('c_Nr'), nth: s('c_nth'),
138 hr: s('c_hr'), hth: s('c_hth'),
139 ht: s('c_ht'), Eh: s('c_Eh'), Rs: s('c_Rs'), rhoh: s('c_rhoh'),
140 fsc: s('c_fsc'), nlr0: s('c_nlr0'), nDropped: s('c_nDropped'),
141 ls: s('c_ls'), js: s('c_js'), lp: s('c_lp'), jp: s('c_jp'),
142 strikeR: s('c_strikeR'), pickupR: s('c_pickupR'),
143 tipi: s('c_tipi'), nfree: s('c_nfree'), ps: s('c_ps'),
144 KM: s('c_KM'), alM: s('c_alM'), aH: s('c_aH'), aHr: s('c_aHr'),
145 holdSteps: s('c_holdSteps'), missDepth: s('c_missDepth'),
146 gap: s('c_gap'), K: s('c_K'), nsteps: s('c_nsteps'),
147 simRate: s('c_simRate'), thickness: s('c_thickness'),
148 };
149 if (setup.vstr.length !== setup.K) throw new Error('schedule length mismatch');
150 return setup;
151}
153// The output stage (DC blocker, decimation, fade, normalize), run through the
154// unmodified cymbal2_output.m with the GPU-computed acceleration handed in.
155export function runOutput(mfiles, params, acc) {
156 const accT = RTV.tensor(Float64Array.from(acc), [1, acc.length]);
157 const r = stage(
158 mfiles,
159 `
160P = cymbal2_defaults();
161${pOverrides(params)}
162[y, peakAccel] = cymbal2_output(P, acc);
163`,
164 { acc: accT },
165 );
166 const y = tensor(r.variableValues.y, 'y');
167 return { y: f32(y), peakAccel: scalar(r.variableValues.peakAccel, 'peakAccel') };
168}
170// CPU reference: the unmodified cymbal2_sim through numbl (slow; for
171// validation of the GPU path on reduced configurations).
172export function runCpuSim(mfiles, params) {
173 const t0 = performance.now();
174 const r = stage(
175 mfiles,
176 `
177P = cymbal2_defaults();
178${pOverrides(params)}
179S = cymbal2_system(P);
180MB = cymbal2_mallet(P);
181[acc, dinfo, fhist] = cymbal2_sim(P, S, MB);
182d_peakF = dinfo.peakForceN;
183d_cvel = dinfo.contactVel(1);
184d_contactMs = dinfo.contactMs;
185d_events = dinfo.contactEvents;
186d_xi = dinfo.peakXiMm;
187d_nsteps = dinfo.nsteps;
188`,
189 );
190 const ms = performance.now() - t0;
191 const v = r.variableValues;
192 return {
193 acc: f32(tensor(v.acc, 'acc')),
194 fhist: f32(tensor(v.fhist, 'fhist')),
195 peakForceN: scalar(v.d_peakF, 'd_peakF'),
196 contactVel: scalar(v.d_cvel, 'd_cvel'),
197 contactMs: scalar(v.d_contactMs, 'd_contactMs'),
198 contactEvents: scalar(v.d_events, 'd_events'),
199 peakXiMm: scalar(v.d_xi, 'd_xi'),
200 nsteps: scalar(v.d_nsteps, 'd_nsteps'),
201 ms,
202 };
203}