5623b11Recorder simulator: MATLAB jet-drive model and interactive web appJeremy Magland 1// Web Worker that renders the recorder 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// basis tier -> recorder_basis per note (transmission-line solve +
6// modal fit; cached on bore/hole/window params and the
7// note's fingering)
8// jet tier -> recorder_jet per note (the time-domain integration;
9// the expensive stage, cached on the breath, the jet,
10// and the synthesis rate)
11// radiate/output -> recorder_radiate per note + assembly + recorder_output
12// (cheap, always re-run: moving the listener costs a
13// matrix multiply, not a re-integration)
14// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).
16import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
17import { hashMFiles, type MFile } from '../mfilesCore';
18import { BASIS_KEYS, JET_KEYS, type RecorderParams } from '../params';
20export interface RenderRequest {
21 type: 'render';
22 id: number;
23 params: RecorderParams;
24 quality: 'preview' | 'full';
25 mfiles: MFile[];
26}
28export interface RenderResponse {
29 type: 'result';
30 id: number;
31 quality: 'preview' | 'full';
32 y: Float32Array; // final signal, [-1, 1]
33 sampleRate: number;
34 noteIdx: Float32Array; // fingering chart row per note
35 noteOn: Float32Array; // onset time of each note, s
36 fPlay: Float32Array; // played frequency per note, Hz
37 fTarget: Float32Array; // equal-tempered target per note, Hz
38 cents: Float32Array; // played vs target, cents
39 qPlot: Float32Array; // jet inflow envelope, decimated, whole timeline
40 qPlotRate: number; // samples per second of qPlot
41 modeF: Float32Array; // first note: fitted resonances, Hz
42 modeQ: Float32Array; // first note: fitted Q per mode
43 gridF: Float32Array; // first note: admittance grid, Hz
44 gridY: Float32Array; // first note: |Y| over the grid
45 profX: Float32Array; // first note: bore stations, m
46 profP: Float32Array; // first note: |p(x)| per mode, nsta x n col-major
47 diag: {
48 noteCount: number;
49 modesKept: number;
50 durationS: number;
51 peakPa: number;
52 peakSplDb: number;
53 listenerDistM: number;
54 warnings: string[];
55 basisMs: number; // 0 when every note's basis was cached
56 jetMs: number; // 0 when every note's jet run was cached
57 radiateMs: number;
58 basisCached: boolean;
59 jetCached: boolean;
60 oversample: number;
61 };
62}
64export interface ProgressMessage {
65 type: 'progress';
66 id: number;
67 stage: 'basis' | 'jet' | 'radiate';
68 frac: number;
69}
71// Raised by a fault in the .m code (the user can edit it), so it is reported
72// rather than left to stall the request.
73export interface ErrorMessage {
74 type: 'error';
75 id: number;
76 message: string;
77 file?: string;
78 line?: number;
79}
81// Basis fields that cross from the basis stage into the jet and radiate
82// stages as numbl RuntimeValues.
83const BASIS_FIELDS = ['f', 'omega', 'Q', 'g', 'B', 'apX', 'apW', 'Sm', 'fTarget'];
85interface BasisEntry {
86 vars: Record<string, RuntimeValue>;
87 modeF: Float32Array;
88 modeQ: Float32Array;
89 gridF: Float32Array;
90 gridY: Float32Array;
91 profX: Float32Array;
92 profP: Float32Array;
93}
95interface JetEntry {
96 Zr: RuntimeValue;
97 Zi: RuntimeValue;
98 Qin: RuntimeValue;
99 fPlay: number;
100 fTarget: number;
101 cents: number;
102 qPlot: Float32Array;
103 qPlotRate: number;
104 fsHi: number;
105}
107const basisCache = new Map<string, BasisEntry>();
108const jetCache = new Map<string, JetEntry>();
109const BASIS_CACHE_MAX = 6;
110const JET_CACHE_MAX = 4;
112function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
113 const v = cache.get(key);
114 if (v !== undefined) {
115 cache.delete(key);
116 cache.set(key, v);
117 }
118 return v;
119}
121function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
122 cache.delete(key);
123 cache.set(key, value);
124 while (cache.size > max) {
125 const oldest = cache.keys().next().value as string;
126 cache.delete(oldest);
127 }
128}
130function keyFor(params: RecorderParams, keys: readonly (keyof RecorderParams)[]): string {
131 return keys.map((k) => String(params[k])).join('|');
132}
134// MATLAB source lines assigning every parameter onto a struct P.
135function pLiteral(params: RecorderParams): string {
136 const lines: string[] = [];
137 for (const [k, v] of Object.entries(params)) {
138 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
139 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
140 }
141 return lines.join('\n');
142}
144function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
145 return executeCode(
146 source,
147 { initialVariableValues: inputs, displayResults: false },
148 mfiles,
149 'stage.m',
150 );
151}
153function asTensor(v: RuntimeValue): RuntimeTensor {
154 const t = v as RuntimeTensor;
155 if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
156 return t;
157}
159function asScalar(v: RuntimeValue): number {
160 if (typeof v === 'number') return v;
161 const t = v as RuntimeTensor;
162 if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
163 if (typeof v === 'boolean') return v ? 1 : 0;
164 throw new Error('expected scalar result from numbl stage');
165}
167function toF32(v: RuntimeValue): Float32Array {
168 const t = asTensor(v);
169 const out = new Float32Array(t.data.length);
170 for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
171 return out;
172}
174const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
176self.onmessage = (ev: MessageEvent<RenderRequest>) => {
177 const msg = ev.data;
178 if (msg.type !== 'render') return;
180 // numbl reports a file it could not parse through console.warn and then
181 // skips it, so the failure surfaces at the call site as "Undefined
182 // function". Capture the warning to report the actual syntax error instead.
183 const warnings: string[] = [];
184 const origWarn = console.warn;
185 console.warn = (...args: unknown[]) => {
186 warnings.push(args.map(String).join(' '));
187 };
189 try {
190 render(msg);
191 } catch (err) {
192 const e = err as { message?: string; file?: string; line?: number };
193 let m: ErrorMessage = {
194 type: 'error',
195 id: msg.id,
196 message: e?.message ? String(e.message) : String(err),
197 file: typeof e?.file === 'string' ? e.file : undefined,
198 line: typeof e?.line === 'number' ? e.line : undefined,
199 };
200 for (const w of warnings) {
201 const hit = SKIP_WARNING.exec(w);
202 if (hit) {
203 m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
204 break;
205 }
206 }
207 (self as unknown as Worker).postMessage(m);
208 } finally {
209 console.warn = origWarn;
210 }
211};
213function render(msg: RenderRequest): void {
214 const { id, params, mfiles } = msg;
215 let quality = msg.quality;
217 const codeKey = hashMFiles(mfiles);
218 const basisParamsKey = codeKey + '|' + keyFor(params, BASIS_KEYS);
219 const jetParamsKey = keyFor(params, JET_KEYS);
221 // The gesture, mirrored from recorder_notes.m.
222 const K = Math.min(3, Math.max(1, Math.round(params.note_count)));
223 const notes = [params.note1, params.note2, params.note3]
224 .slice(0, K)
225 .map((n) => Math.min(14, Math.max(1, Math.round(n))));
226 const gap = Math.max(0, params.gap_ms) * 1e-3;
227 const tOn = notes.map((_, i) => i * (params.note_dur_s + gap));
228 const spans = notes.map((_, i) =>
229 params.note_dur_s + (i < K - 1 ? Math.min(params.tail_s, 0.3) : params.tail_s),
230 );
231 const durationS = tOn[K - 1] + params.note_dur_s + params.tail_s;
233 const jetKeyAt = (i: number, os: number) =>
234 basisParamsKey + `|n=${notes[i]}` + '||' + jetParamsKey + `|span=${spans[i]}|ord=${i + 1}|os=${os}`;
236 // A preview right after a full render reuses the full-rate jet runs.
237 if (quality === 'preview' && notes.every((_, i) => jetCache.has(jetKeyAt(i, params.oversample)))) {
238 quality = 'full';
239 }
240 const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;
241 const fsHi = params.sample_rate * oversample;
243 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
245 let basisMs = 0;
246 let jetMs = 0;
247 let basisCached = true;
248 let jetCached = true;
249 const basisEntries: BasisEntry[] = [];
250 const jetEntries: JetEntry[] = [];
252 for (let i = 0; i < K; i++) {
253 // ---- basis stage: transmission-line solve + modal fit ----
254 const bKey = basisParamsKey + `|n=${notes[i]}`;
255 let bEntry = lruGet(basisCache, bKey);
256 if (!bEntry) {
257 basisCached = false;
258 post({ type: 'progress', id, stage: 'basis', frac: (i + 0.5) / K * 0.2 });
259 const t0 = performance.now();
260 const r = runStage(mfiles, `
261${pLiteral(params)}
262basis = recorder_basis(P, ${notes[i]});
263${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
264b_profP = basis.profP;
265b_profX = basis.profX;
266b_gridF = basis.gridF;
267b_gridY = basis.gridY;
268`);
269 basisMs += performance.now() - t0;
270 const vars: Record<string, RuntimeValue> = {};
271 for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];
272 bEntry = {
273 vars,
274 modeF: toF32(r.variableValues.b_f),
275 modeQ: toF32(r.variableValues.b_Q),
276 gridF: toF32(r.variableValues.b_gridF),
277 gridY: toF32(r.variableValues.b_gridY),
278 profX: toF32(r.variableValues.b_profX),
279 profP: toF32(r.variableValues.b_profP),
280 };
281 lruSet(basisCache, bKey, bEntry, BASIS_CACHE_MAX);
282 }
283 basisEntries.push(bEntry);
285 // ---- jet stage: the time-domain integration ----
286 const jKey = jetKeyAt(i, oversample);
287 let jEntry = lruGet(jetCache, jKey);
288 if (!jEntry) {
289 jetCached = false;
290 post({ type: 'progress', id, stage: 'jet', frac: 0.2 + ((i + 0.5) / K) * 0.6 });
291 const t0 = performance.now();
292 const ds = Math.max(1, Math.round(fsHi / 2000));
293 const r = runStage(mfiles, `
294${pLiteral(params)}
295${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}
296basis.n = length(basis.f);
297[Zr, Zi, Qin, jinfo] = recorder_jet(P, basis, ${spans[i]}, ${i + 1}, ${fsHi});
298fPlay = jinfo.fPlay;
299fTarget = jinfo.fTarget;
300cents = jinfo.cents;
301qPlot = Qin(1:${ds}:length(Qin));
302`, bEntry.vars);
303 jetMs += performance.now() - t0;
304 jEntry = {
305 Zr: r.variableValues.Zr,
306 Zi: r.variableValues.Zi,
307 Qin: r.variableValues.Qin,
308 fPlay: asScalar(r.variableValues.fPlay),
309 fTarget: asScalar(r.variableValues.fTarget),
310 cents: asScalar(r.variableValues.cents),
311 qPlot: toF32(r.variableValues.qPlot),
312 qPlotRate: fsHi / ds,
313 fsHi,
314 };
315 lruSet(jetCache, jKey, jEntry, JET_CACHE_MAX);
316 }
317 jetEntries.push(jEntry);
318 }
320 // ---- radiate + assemble + output ----
321 post({ type: 'progress', id, stage: 'radiate', frac: 0.9 });
322 const t0 = performance.now();
323 const inputs: Record<string, RuntimeValue> = {};
324 const noteBlocks: string[] = [];
325 for (let i = 0; i < K; i++) {
326 for (const f of BASIS_FIELDS) inputs[`b${i + 1}_${f}`] = basisEntries[i].vars[`b_${f}`];
327 inputs[`n${i + 1}_Zr`] = jetEntries[i].Zr;
328 inputs[`n${i + 1}_Zi`] = jetEntries[i].Zi;
329 inputs[`n${i + 1}_Qin`] = jetEntries[i].Qin;
330 noteBlocks.push(`
331${BASIS_FIELDS.map((f) => `basis.${f} = b${i + 1}_${f};`).join('\n')}
332basis.n = length(basis.f);
333[pn, radinfo] = recorder_radiate(P, basis, n${i + 1}_Zr, n${i + 1}_Zi, n${i + 1}_Qin, ${fsHi});
334k0 = min(T, round(${tOn[i]} * ${fsHi}) + 1);
335len = min(length(pn), T - k0 + 1);
336p(k0:k0 + len - 1) = p(k0:k0 + len - 1) + pn(1:len);
337inside = max(inside, radinfo.listenerInside);
338distM = radinfo.listenerDistM;
339`);
340 }
341 const rOut = runStage(mfiles, `
342${pLiteral(params)}
343T = floor(${durationS} * ${fsHi});
344p = zeros(1, T);
345inside = 0;
346distM = 0;
347${noteBlocks.join('\n')}
348[y, peakPa] = recorder_output(P, p, ${oversample});
349`, inputs);
350 const radiateMs = performance.now() - t0;
352 const y32 = toF32(rOut.variableValues.y);
353 const peakPa = asScalar(rOut.variableValues.peakPa);
354 const listenerInside = asScalar(rOut.variableValues.inside) !== 0;
355 const listenerDistM = asScalar(rOut.variableValues.distM);
357 // assemble the jet-flow envelope over the whole timeline
358 const qRate = jetEntries[0].qPlotRate;
359 const qPlot = new Float32Array(Math.ceil(durationS * qRate));
360 for (let i = 0; i < K; i++) {
361 const off = Math.round(tOn[i] * qRate);
362 const q = jetEntries[i].qPlot;
363 for (let k = 0; k < q.length && off + k < qPlot.length; k++) qPlot[off + k] += q[k];
364 }
366 const warns: string[] = [];
367 if (listenerInside) {
368 warns.push('the listener is within 2 cm of the instrument; the monopole field is evaluated at a clamped minimum distance');
369 }
370 for (let i = 0; i < K; i++) {
371 const j = jetEntries[i];
372 if (j.fPlay <= 0) {
373 warns.push(`note ${i + 1} did not speak (no oscillation locked); more breath or jet receptivity may help`);
374 } else if (Math.abs(j.cents) > 60) {
375 warns.push(`note ${i + 1} locked to ${j.fPlay.toFixed(0)} Hz, not its fingering's register (${j.fTarget.toFixed(0)} Hz)`);
376 }
377 }
378 if (!params.normalize) {
379 let peak = 0;
380 for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
381 if (peak > 1.0) {
382 warns.push(`clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`);
383 }
384 }
386 const resp: RenderResponse = {
387 type: 'result',
388 id,
389 quality,
390 y: y32,
391 sampleRate: params.sample_rate,
392 noteIdx: new Float32Array(notes),
393 noteOn: new Float32Array(tOn),
394 fPlay: new Float32Array(jetEntries.map((j) => j.fPlay)),
395 fTarget: new Float32Array(jetEntries.map((j) => j.fTarget)),
396 cents: new Float32Array(jetEntries.map((j) => j.cents)),
397 qPlot,
398 qPlotRate: qRate,
399 modeF: basisEntries[0].modeF.slice(),
400 modeQ: basisEntries[0].modeQ.slice(),
401 gridF: basisEntries[0].gridF.slice(),
402 gridY: basisEntries[0].gridY.slice(),
403 profX: basisEntries[0].profX.slice(),
404 profP: basisEntries[0].profP.slice(),
405 diag: {
406 noteCount: K,
407 modesKept: basisEntries[0].modeF.length,
408 durationS,
409 peakPa,
410 peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
411 listenerDistM,
412 warnings: warns,
413 basisMs,
414 jetMs,
415 radiateMs,
416 basisCached,
417 jetCached,
418 oversample,
419 },
420 };
421 (self as unknown as Worker).postMessage(resp, [
422 resp.y.buffer,
423 resp.noteIdx.buffer,
424 resp.noteOn.buffer,
425 resp.fPlay.buffer,
426 resp.fTarget.buffer,
427 resp.cents.buffer,
428 resp.qPlot.buffer,
429 resp.modeF.buffer,
430 resp.modeQ.buffer,
431 resp.gridF.buffer,
432 resp.gridY.buffer,
433 resp.profX.buffer,
434 resp.profP.buffer,
435 ]);
436}