webapp/src/worker/renderWorker.ts at main - simulated

simulated_instruments / cymbal

502 lines · 17.0 KBBlame History Raw

1// Web Worker that renders the cymbal stroke 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 -> cymbal_basis (one shell eigensolve per nodal-diameter

6// order; the expensive stage, cached on geometry, material,

7// damping and discretization parameters)

8// radiate tier -> cymbal_strikes + cymbal_cascade + cymbal_radiate (cached

9// on the stick, the gesture, the cascade, the listener and

10// the synthesis rate)

11// output tier -> cymbal_output (always cheap, always re-run)

12// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).

14import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';

15import { hashMFiles, type MFile } from '../mfiles';

16import { SOLVE_KEYS, RADIATE_KEYS, type CymbalParams } from '../params';

18export interface RenderRequest {

19 type: 'render';

20 id: number;

21 params: CymbalParams;

22 quality: 'preview' | 'full';

23 mfiles: MFile[];

24}

26export interface RenderResponse {

27 type: 'result';

28 id: number;

29 quality: 'preview' | 'full';

30 y: Float32Array; // final signal, [-1, 1]

31 sampleRate: number;

32 strikeT: Float32Array; // stroke times, s

33 strikeV: Float32Array; // stroke velocities, m/s

34 strikeS: Float32Array; // stroke arclength fractions

35 modeF: Float32Array; // mode frequencies, Hz

36 modeZeta: Float32Array; // damping ratio of each mode

37 modeZetaRad: Float32Array; // radiation part of the damping

38 modeSigma: Float32Array; // radiation efficiency of each mode

39 modeN: Float32Array; // nodal diameters

40 nodeR: Float32Array; // meridian nodes, m

41 nodeZ: Float32Array;

42 shapeW: Float32Array; // normal displacement per mode, nn x n col-major

43 cascadeDb: Float32Array; // energy grid, nBands x nFrames col-major, dB

44 excDb: Float32Array; // first stroke's modal excitation per band, dB

45 cascadeBands: number;

46 cascadeFrames: number;

47 cascadeFLo: number;

48 cascadeFHi: number;

49 cascadeTEnd: number;

50 diag: {

51 modesKept: number;

52 maxN: number;

53 strikeCount: number;

54 meanV: number;

55 maxV: number;

56 meanTauUs: number;

57 peakForceN: number;

58 durationS: number;

59 peakPa: number;

60 peakSplDb: number;

61 directPeakPa: number;

62 listenerDistM: number;

63 strikeRadiusM: number;

64 strikeRimM: number;

65 roll2RadiusM: number;

66 peakWrmsUm: number;

67 peakStrength: number;

68 cascadeLostFrac: number;

69 plateMassKg: number;

70 warnings: string[];

71 solveMs: number; // 0 when the cached eigensolve was reused

72 radiateMs: number; // 0 when the cached radiation was reused

73 solveCached: boolean;

74 radiateCached: boolean;

75 oversample: number;

76 };

77}

79export interface ProgressMessage {

80 type: 'progress';

81 id: number;

82 stage: 'solve' | 'radiate' | 'output';

83 frac: number;

84}

86// Raised by a fault in the .m code (the user can edit it), so it is reported

87// rather than left to stall the request.

88export interface ErrorMessage {

89 type: 'error';

90 id: number;

91 message: string;

92 file?: string;

93 line?: number;

94}

96interface SolveEntry {

97 vars: Record<string, RuntimeValue>;

98 modeF: Float32Array;

99 modeZeta: Float32Array;

100 modeZetaRad: Float32Array;

101 modeSigma: Float32Array;

102 modeN: Float32Array;

103 nodeR: Float32Array;

104 nodeZ: Float32Array;

105 shapeW: Float32Array;

106 maxN: number;

107 plateMassKg: number;

108}

109

110interface RadiateEntry {

111 p: RuntimeValue;

112 strikeT: Float32Array;

113 strikeV: Float32Array;

114 strikeS: Float32Array;

115 cascadeDb: Float32Array;

116 excDb: Float32Array;

117 cascadeBands: number;

118 cascadeFrames: number;

119 cascadeTEnd: number;

120 strikeRadiusM: number;

121 strikeRimM: number;

122 roll2RadiusM: number;

123 strikeCount: number;

124 meanV: number;

125 maxV: number;

126 meanTauUs: number;

127 peakForceN: number;

128 durationS: number;

129 directPeakPa: number;

130 listenerDistM: number;

131 listenerInside: boolean;

132 peakWrmsUm: number;

133 peakStrength: number;

134 cascadeLostFrac: number;

135 oversample: number;

136}

137

138// Small LRU caches. The solve cache does not depend on the oversample rate,

139// so a preview render right after a full one reuses the same eigensolves.

140const solveCache = new Map<string, SolveEntry>();

141const radiateCache = new Map<string, RadiateEntry>();

142const SOLVE_CACHE_MAX = 2;

143const RADIATE_CACHE_MAX = 3;

144

145function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {

146 const v = cache.get(key);

147 if (v !== undefined) {

148 cache.delete(key);

149 cache.set(key, v);

150 }

151 return v;

152}

153

154function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {

155 cache.delete(key);

156 cache.set(key, value);

157 while (cache.size > max) {

158 const oldest = cache.keys().next().value as string;

159 cache.delete(oldest);

160 }

161}

162

163function keyFor(params: CymbalParams, keys: readonly (keyof CymbalParams)[]): string {

164 return keys.map((k) => String(params[k])).join('|');

165}

166

167// MATLAB source lines assigning every parameter onto a struct P.

168function pLiteral(params: CymbalParams): string {

169 const lines: string[] = [];

170 for (const [k, v] of Object.entries(params)) {

171 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);

172 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);

173 }

174 return lines.join('\n');

175}

176

177function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {

178 return executeCode(

179 source,

180 { initialVariableValues: inputs, displayResults: false },

181 mfiles,

182 'stage.m',

183 );

184}

185

186function asTensor(v: RuntimeValue): RuntimeTensor {

187 const t = v as RuntimeTensor;

188 if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');

189 return t;

190}

191

192function asScalar(v: RuntimeValue): number {

193 if (typeof v === 'number') return v;

194 const t = v as RuntimeTensor;

195 if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];

196 if (typeof v === 'boolean') return v ? 1 : 0;

197 throw new Error('expected scalar result from numbl stage');

198}

199

200function toF32(v: RuntimeValue): Float32Array {

201 const t = asTensor(v);

202 const out = new Float32Array(t.data.length);

203 for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];

204 return out;

205}

206

207// The basis fields handed from the solve stage into the radiate stage, in the

208// order the reassembly script expects.

209const BASIS_FIELDS = [

210 'f', 'omega', 'zeta', 'zetaRad', 'zetaMount', 'sigma', 'ndia', 'ldip', 'mix',

211 'NU', 'NW', 'NT',

212 'rnode', 'znode', 'snode', 'drdsnode', 'dzdsnode', 'thknode', 'csnode', 'ctnode',

213 'surfR', 'surfZ', 'surfDr', 'surfDz', 'surfLen', 'SW',

214 'S', 'area', 'mass', 'thkMean',

215];

216

217const CASCADE_BANDS = 48;

218const CASCADE_FRAMES = 240;

219

220const SKIP_WARNING = /skipping (\S+\.m) $syntax error at line (\d+)$/;

221

222self.onmessage = (ev: MessageEvent<RenderRequest>) => {

223 const msg = ev.data;

224 if (msg.type !== 'render') return;

225

226 // numbl reports a file it could not parse through console.warn and then

227 // skips it, so the failure surfaces at the call site as "Undefined

228 // function". Capture the warning to report the actual syntax error instead.

229 const warnings: string[] = [];

230 const origWarn = console.warn;

231 console.warn = (...args: unknown[]) => {

232 warnings.push(args.map(String).join(' '));

233 };

234

235 try {

236 render(msg);

237 } catch (err) {

238 const e = err as { message?: string; file?: string; line?: number };

239 let m: ErrorMessage = {

240 type: 'error',

241 id: msg.id,

242 message: e?.message ? String(e.message) : String(err),

243 file: typeof e?.file === 'string' ? e.file : undefined,

244 line: typeof e?.line === 'number' ? e.line : undefined,

245 };

246 for (const w of warnings) {

247 const hit = SKIP_WARNING.exec(w);

248 if (hit) {

249 m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };

250 break;

251 }

252 }

253 (self as unknown as Worker).postMessage(m);

254 } finally {

255 console.warn = origWarn;

256 }

257};

258

259function render(msg: RenderRequest): void {

260 const { id, params, mfiles } = msg;

261 let quality = msg.quality;

262

263 // The .m sources are part of the cache identity: editing the code must

264 // invalidate everything computed from the previous version.

265 const codeKey = hashMFiles(mfiles);

266 const solveKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);

267

268 const fullOversample = params.oversample;

269 const radiateParamsKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);

270 if (quality === 'preview' && radiateCache.has(radiateParamsKey + `|os=${fullOversample}`)) {

271 quality = 'full'; // full-rate radiation already cached

272 }

273 const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;

274 const fsHi = params.sample_rate * oversample;

275

276 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);

277

278 // ---- solve stage: one shell eigensolve per nodal-diameter order ----

279 let solveMs = 0;

280 let solveCached = true;

281 let solveEntry = lruGet(solveCache, solveKey);

282 if (!solveEntry) {

283 solveCached = false;

284 post({ type: 'progress', id, stage: 'solve', frac: 0.12 });

285 const t0 = performance.now();

286 const r = runStage(mfiles, `

287${pLiteral(params)}

288basis = cymbal_basis(P);

289${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}

290maxN = max(basis.ndia);

291`);

292 solveMs = performance.now() - t0;

293

294 const vars: Record<string, RuntimeValue> = {};

295 for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];

296 solveEntry = {

297 vars,

298 modeF: toF32(r.variableValues.b_f),

299 modeZeta: toF32(r.variableValues.b_zeta),

300 modeZetaRad: toF32(r.variableValues.b_zetaRad),

301 modeSigma: toF32(r.variableValues.b_sigma),

302 modeN: toF32(r.variableValues.b_ndia),

303 nodeR: toF32(r.variableValues.b_rnode),

304 nodeZ: toF32(r.variableValues.b_znode),

305 shapeW: toF32(r.variableValues.b_NW),

306 maxN: asScalar(r.variableValues.maxN),

307 plateMassKg: asScalar(r.variableValues.b_mass),

308 };

309 lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);

310 }

311

312 // ---- radiate stage: strikes + cascade + modal synthesis ----

313 const radiateKey = radiateParamsKey + `|os=${oversample}`;

314 let radiateMs = 0;

315 let radiateCached = true;

316 let radEntry = lruGet(radiateCache, radiateKey);

317 if (!radEntry) {

318 radiateCached = false;

319 post({ type: 'progress', id, stage: 'radiate', frac: 0.5 });

320 const t0 = performance.now();

321 const r = runStage(mfiles, `

322${pLiteral(params)}

323${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}

324basis.n = length(basis.f);

325[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = cymbal_strikes(P, basis);

326[env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, taustr, sstr, sinfo.durationS);

327[p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, ${fsHi});

328% collapse the modal energies onto a log-frequency by time grid for the plot

329nb = ${CASCADE_BANDS};

330nt = ${CASCADE_FRAMES};

331fLo = max(1, min(basis.f));

332fHi = max(fLo * 1.2, max(basis.f));

333bandOf = min(nb, max(1, 1 + floor(nb * (log(basis.f) - log(fLo)) / (log(fHi) - log(fLo) + 1e-12))));

334nf = size(env, 2);

335grid = zeros(nb, nt);

336for q = 1:nt

337 c0 = 1 + floor((q - 1) * nf / nt);

338 c1 = max(c0, floor(q * nf / nt));

339 e = sum(env(:, c0:c1).^2, 2) / (c1 - c0 + 1);

340 grid(:, q) = accumarray(bandOf, e, [nb, 1]);

341end

342cascadeDb = 10 * log10(grid + 1e-30);

343excDb = 10 * log10(accumarray(bandOf, 0.5 * DV(:, 1).^2, [nb, 1]) + 1e-30);

344cascadeTEnd = tframe(length(tframe));

345[strikeRadiusM, ~] = cymbal_point(basis, P.strike_r);

346strikeRimM = P.radius_m - strikeRadiusM;

347[roll2RadiusM, ~] = cymbal_point(basis, P.roll_r2);

348strikeCount = sinfo.count;

349meanV = sinfo.meanV;

350maxV = sinfo.maxV;

351meanTauUs = sinfo.meanTauUs;

352peakForceN = sinfo.peakForceN;

353durationS = sinfo.durationS;

354directPeakPa = radinfo.directPeakPa;

355listenerDistM = radinfo.listenerDistM;

356listenerInside = radinfo.listenerInside;

357peakWrmsUm = cinfo.peakWrmsUm;

358peakStrength = cinfo.peakStrength;

359cascadeLostFrac = cinfo.cascadeLostJ / max(1e-30, cinfo.injectedJ);

360`, solveEntry.vars);

361 radiateMs = performance.now() - t0;

362 radEntry = {

363 p: r.variableValues.p,

364 strikeT: toF32(r.variableValues.tstr),

365 strikeV: toF32(r.variableValues.vstr),

366 strikeS: toF32(r.variableValues.sstr),

367 cascadeDb: toF32(r.variableValues.cascadeDb),

368 excDb: toF32(r.variableValues.excDb),

369 cascadeBands: CASCADE_BANDS,

370 cascadeFrames: CASCADE_FRAMES,

371 cascadeTEnd: asScalar(r.variableValues.cascadeTEnd),

372 strikeRadiusM: asScalar(r.variableValues.strikeRadiusM),

373 strikeRimM: asScalar(r.variableValues.strikeRimM),

374 roll2RadiusM: asScalar(r.variableValues.roll2RadiusM),

375 strikeCount: asScalar(r.variableValues.strikeCount),

376 meanV: asScalar(r.variableValues.meanV),

377 maxV: asScalar(r.variableValues.maxV),

378 meanTauUs: asScalar(r.variableValues.meanTauUs),

379 peakForceN: asScalar(r.variableValues.peakForceN),

380 durationS: asScalar(r.variableValues.durationS),

381 directPeakPa: asScalar(r.variableValues.directPeakPa),

382 listenerDistM: asScalar(r.variableValues.listenerDistM),

383 listenerInside: asScalar(r.variableValues.listenerInside) !== 0,

384 peakWrmsUm: asScalar(r.variableValues.peakWrmsUm),

385 peakStrength: asScalar(r.variableValues.peakStrength),

386 cascadeLostFrac: asScalar(r.variableValues.cascadeLostFrac),

387 oversample,

388 };

389 lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);

390 }

391

392 // ---- output stage ----

393 post({ type: 'progress', id, stage: 'output', frac: 0.9 });

394 const rOut = runStage(mfiles, `

395${pLiteral(params)}

396[y, peakPa] = cymbal_output(P, p, ${radEntry.oversample});

397`, { p: radEntry.p });

398

399 const y32 = toF32(rOut.variableValues.y);

400 const peakPa = asScalar(rOut.variableValues.peakPa);

401

402 const warnings: string[] = [];

403 if (radEntry.listenerInside) {

404 warnings.push(

405 'the listener is within 2 cm of the plate; the dipole field is evaluated at a clamped minimum distance',

406 );

407 }

408 if (!params.normalize) {

409 let peak = 0;

410 for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));

411 if (peak > 1.0) {

412 warnings.push(

413 `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,

414 );

415 }

416 }

417 if (params.n_azimuth < 2.5 * solveEntry.maxN) {

418 warnings.push(

419 `only ${params.n_azimuth} azimuth samples for ${solveEntry.maxN} nodal diameters; the listener transfer of the high-order modes is aliased`,

420 );

421 }

422 if (radEntry.cascadeLostFrac > 0.25) {

423 warnings.push(

424 `${(100 * radEntry.cascadeLostFrac).toFixed(0)}% of the struck energy left the top of the retained band; raise max_mode_hz or lower cascade_leak`,

425 );

426 }

427 const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];

428 if (topMode < 0.6 * params.max_mode_hz) {

429 warnings.push(

430 `the highest retained mode is only ${topMode.toFixed(0)} Hz; raise max nodal diameters to fill the band`,

431 );

432 }

433

434 const resp: RenderResponse = {

435 type: 'result',

436 id,

437 quality,

438 y: y32,

439 sampleRate: params.sample_rate,

440 strikeT: radEntry.strikeT.slice(),

441 strikeV: radEntry.strikeV.slice(),

442 strikeS: radEntry.strikeS.slice(),

443 modeF: solveEntry.modeF.slice(),

444 modeZeta: solveEntry.modeZeta.slice(),

445 modeZetaRad: solveEntry.modeZetaRad.slice(),

446 modeSigma: solveEntry.modeSigma.slice(),

447 modeN: solveEntry.modeN.slice(),

448 nodeR: solveEntry.nodeR.slice(),

449 nodeZ: solveEntry.nodeZ.slice(),

450 shapeW: solveEntry.shapeW.slice(),

451 cascadeDb: radEntry.cascadeDb.slice(),

452 excDb: radEntry.excDb.slice(),

453 cascadeBands: radEntry.cascadeBands,

454 cascadeFrames: radEntry.cascadeFrames,

455 cascadeFLo: solveEntry.modeF[0],

456 cascadeFHi: solveEntry.modeF[solveEntry.modeF.length - 1],

457 cascadeTEnd: radEntry.cascadeTEnd,

458 diag: {

459 modesKept: solveEntry.modeF.length,

460 maxN: solveEntry.maxN,

461 strikeCount: radEntry.strikeCount,

462 meanV: radEntry.meanV,

463 maxV: radEntry.maxV,

464 meanTauUs: radEntry.meanTauUs,

465 peakForceN: radEntry.peakForceN,

466 durationS: radEntry.durationS,

467 peakPa,

468 peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),

469 directPeakPa: radEntry.directPeakPa,

470 listenerDistM: radEntry.listenerDistM,

471 strikeRadiusM: radEntry.strikeRadiusM,

472 strikeRimM: radEntry.strikeRimM,

473 roll2RadiusM: radEntry.roll2RadiusM,

474 peakWrmsUm: radEntry.peakWrmsUm,

475 peakStrength: radEntry.peakStrength,

476 cascadeLostFrac: radEntry.cascadeLostFrac,

477 plateMassKg: solveEntry.plateMassKg,

478 warnings,

479 solveMs,

480 radiateMs,

481 solveCached,

482 radiateCached,

483 oversample: radEntry.oversample,

484 },

485 };

486 (self as unknown as Worker).postMessage(resp, [

487 resp.y.buffer,

488 resp.strikeT.buffer,

489 resp.strikeV.buffer,

490 resp.strikeS.buffer,

491 resp.modeF.buffer,

492 resp.modeZeta.buffer,

493 resp.modeZetaRad.buffer,

494 resp.modeSigma.buffer,

495 resp.modeN.buffer,

496 resp.nodeR.buffer,

497 resp.nodeZ.buffer,

498 resp.shapeW.buffer,

499 resp.cascadeDb.buffer,

500 resp.excDb.buffer,

501 ]);

502}