Recorder simulator: MATLAB jet-drive model and interactive web app
46 changed files+8450−0
.github/workflows/deploy.ymladded+70−0View file
@@ -0,0 +1,70 @@
1+name: Deploy to Pages
2+
3+# Builds webapp/ and publishes it as this repository's site.
4+#
5+# Two notes on the vault:
6+# - there is no cache service, so setup-node runs without `cache: npm`;
7+# - a mochi site is served at /<collection>/<repo>/site/ or at the root of an
8+# origin of its own. The Vite build uses a relative base (`base: './'`), so
9+# it works at either without being told the deploy path, and configure-pages
10+# is therefore only here for parity with the deploy step.
11+
12+on:
13+ push:
14+ branches: [main]
15+ workflow_dispatch:
16+
17+permissions:
18+ contents: read
19+ pages: write
20+ id-token: write
21+
22+concurrency:
23+ group: pages
24+ cancel-in-progress: false
25+
26+jobs:
27+ build:
28+ runs-on: ubuntu-latest
29+ steps:
30+ - uses: actions/checkout@v4
31+
32+ - uses: actions/setup-node@v4
33+ with:
34+ node-version: 22
35+
36+ - name: Install
37+ working-directory: webapp
38+ run: npm ci
39+
40+ # Runs matlab/recorder_validate.m through numbl: the closed-form
41+ # open-open cylinder, the modal-fit identity, emergent register
42+ # selection and overblowing, reproducibility of the seeded turbulence,
43+ # and the pinned end-to-end diagnostics. Fails the build on any
44+ # mismatch.
45+ - name: Validate the model
46+ working-directory: webapp
47+ run: npm run selftest
48+
49+ - name: Build
50+ working-directory: webapp
51+ run: npm run build
52+
53+ - name: Setup Pages
54+ uses: actions/configure-pages@v5
55+
56+ - name: Upload artifact
57+ uses: actions/upload-pages-artifact@v3
58+ with:
59+ path: webapp/dist
60+
61+ deploy:
62+ needs: build
63+ runs-on: ubuntu-latest
64+ environment:
65+ name: github-pages
66+ url: ${{ steps.deployment.outputs.page_url }}
67+ steps:
68+ - name: Deploy to Pages
69+ id: deployment
70+ uses: actions/deploy-pages@v4
README.mdadded+48−0View file
@@ -0,0 +1,48 @@
1+# Recorder
2+
3+A physically-based simulation of a soprano recorder playing one to three tongued notes, sampled as acoustic pressure at a point in space. The instrument is a tapered bore with eight tone holes under a baroque fingering chart; an air jet from the windway is displaced by the acoustic flow through the window, amplified along its travel to the labium, and the flow it injects drives the bore resonances that displace it, so the tone is a self-sustained oscillation of the loop rather than anything prescribed. The bore is solved as a lossy transmission line, and every open aperture radiates to a listener point.
4+
5+There is no microphone model: no directivity, no capsule response, no gain staging. The output is p(r, t) in pascals.
6+
7+## Contents
8+
9+| Path | What it is |
10+|---|---|
11+| `matlab/` | The algorithm in MATLAB syntax, run by [numbl](https://numbl.org) |
12+| `webapp/` | Interactive browser app that executes the MATLAB code client-side |
13+
14+The MATLAB code runs from the command line with `cd matlab && numbl run recorder_demo.m`, or `recorder_validate.m` to run the checks. The web app is a Vite project; see `webapp/README.md`.
15+
16+## The model
17+
18+The signal chain, with each stage in its own `.m` file:
19+
20+```
21+fingering -> tone hole lattice on the tapered bore recorder_geometry.m
22+ -> lossy transmission line -> input admittance recorder_bore.m
23+ -> modal fit of the resonances recorder_basis.m
24+notes of the gesture recorder_notes.m
25+ -> breath envelope -> jet / bore loop, integrated
26+ sample by sample (the sound happens here) recorder_jet.m
27+ -> coherent sum of aperture monopoles recorder_radiate.m
28+ -> free-field pressure at the listener
29+ -> anti-aliased decimation -> PCM recorder_output.m
30+```
31+
32+Three points are worth stating plainly, because they are where the model commits itself.
33+
34+**The passive instrument is a transmission line with real tone holes.** The bore is a chain of short lossy cylinders (viscothermal wall losses through the standard wide-tube propagation constant, a linear taper standing in for the maker's reaming), each open hole a shunt air mass with end corrections and radiation resistance, each closed hole a shunt compliance, the foot an unflanged open end, and the window an air mass in series. What the jet sees is the admittance of that loop, refitted as a handful of second-order resonators whenever the fingering changes (`recorder_basis.m`). The hole positions and sizes were tuned numerically so that the ten first-register fingerings of the baroque chart land on an equal-tempered C major scale within a cent of passive resonance; the pinched thumb hole is modeled as a small open hole with the viscous resistance of a narrow crack, which is what damps the fundamental and lets the second register speak. Cross-fingerings (the F5 fork) and the pinched fingerings inherit every simplification here, and the forked F6 plays some 40 cents flat, which is honestly the sort of note it is on the instrument.
35+
36+**The oscillation is emergent, not prescribed.** `recorder_jet.m` integrates the jet-drive loop of Verge and Fabre sample by sample at the oversampled rate: breath pressure sets the jet velocity through Bernoulli, the acoustic flow through the window deflects the jet, the deflection convects to the labium (a delay that varies with breath), grows exponentially, and is low-passed at the jet's own instability cutoff; the flow split at the labium saturates as a tanh, and its injection drives the modal resonators through the jet-drive pressure source, closed by a nonlinear vortex loss. Nothing in the loop says which note to play. The note grows out of seeded turbulence noise where the loop gain crosses one, locks near the bore resonance whose quarter-period matches the jet transit time, sharpens as the breath rises, and overblows to the octave when the transit time leaves the fundamental's range; the validation script checks that register selection and overblowing actually emerge rather than assuming they do. The trade-off is that this stage cannot be closed-form, so a two-second note takes a few seconds to render where the sibling instruments render faster than real time. Two stated limitations: notes are rendered independently and summed, which is right for tongued articulation and wrong for legato, and there is no player, so nothing listens to the pitch and corrects the breath (the `breath_track` parameter stands in for the player's habit of blowing higher notes harder).
37+
38+**Radiation is a coherent sum of aperture monopoles.** The window, every open tone hole, and the foot each carry a complex flow per unit loop flow, taken from the transmission-line solve at each mode's centre frequency, and each radiates as a point monopole with its own retardation to the listener, so the interference between apertures, and its dependence on where the listener stands, is kept. The broadband remainder, breath noise, the attack transient, and everything the modal truncation dropped, radiates from the window as the direct term; that is the least principled term of the model, and `direct_gain` turns it off.
39+
40+## Validation
41+
42+`matlab/recorder_validate.m` (also `npm run selftest` in `webapp/`) checks the transmission line against the closed-form open-open cylinder (with the viscothermal dispersion the closed form needs), the modal fit against the admittance it was fitted to, emergent register selection on three fingerings including the all-closed C5 (the easiest fingering to squeak), emergent overblowing at quadrupled breath pressure, exact reproducibility of the seeded turbulence, the far-field -6 dB per doubling of distance, and pinned end-to-end diagnostics.
43+
44+As in the sibling `maraca` and `triangle` repositories, there is no separate Python reference implementation: the MATLAB code is the reference, and the pinned numbers were recorded from it when the model was built.
45+
46+## Defaults
47+
48+A soprano recorder in C: a 30.2 cm bore tapering from 6.7 to 5.2 mm radius, eight tone holes, a 4.5 mm cutup and a 1.0 mm flue. The default gesture is a single G5 at 668 Pa breath (350 Pa at C5, tracked up with pitch), which locks at 781.9 Hz, five cents flat of equal temperament, over twelve fitted bore resonances; the peak is 1.80 Pa (99 dB SPL) at the 44 cm listener before normalization. Setting `note_count` to 3 plays a tongued G5 A5 B5 figure. The full chart runs C5 to A6, the top four notes on pinched second-register fingerings that want more breath, as they do on the instrument.
matlab/recorder_basis.madded+96−0View file
@@ -0,0 +1,96 @@
1+function basis = recorder_basis(P, note)
2+% RECORDER_BASIS Modal fit of the passive resonator for one fingering.
3+%
4+% The series-loop admittance Y(omega) = 1 / (Z_window + Z_bore) from
5+% recorder_bore.m relates the acoustic flow U through the mouth to the
6+% pressure the jet imposes across it; its peaks are the bore resonances the
7+% jet can lock to. Each peak is fitted with one second-order section,
8+% Y_n(omega) = i omega B_n / (omega_n^2 - omega^2 + i omega g_n),
9+% so that in the time domain U_n'' + g_n U_n' + omega_n^2 U_n = B_n dp/dt.
10+% omega_n comes from parabolic interpolation of log|Y| on the grid, g_n from
11+% the half-power bandwidth, and B_n = g_n Re(Y(omega_n)) so the fitted peak
12+% height matches. The fit discards nothing a recorder player can reach:
13+% every local maximum below max_mode_hz is kept, up to max_modes, however
14+% weak (a cross-fingered lattice makes the upper resonances weak but they
15+% still color the attack).
16+%
17+% For radiation the basis stores, per mode, the complex flow out of each
18+% radiating aperture (the window, each open hole, the foot) per unit loop
19+% flow, evaluated at the mode's centre frequency, plus the pressure
20+% standing-wave profile along the bore for the mode-shape plot.
21+
22+geo = recorder_geometry(P);
23+note = min(max(1, round(note)), geo.nNotes);
24+fing = geo.fingerings(note, :);
25+
26+fgrid = 60:P.grid_hz:P.max_mode_hz;
27+out = recorder_bore(P, geo, fing, fgrid);
28+mag = abs(out.Y);
29+
30+% ---- peak picking with parabolic refinement ----
31+fpk = [];
32+qpk = [];
33+ypk = [];
34+nf = length(fgrid);
35+for i = 2:nf - 1
36+ if ~(mag(i) > mag(i - 1) && mag(i) >= mag(i + 1))
37+ continue;
38+ end
39+ y1 = log(mag(i - 1)); y2 = log(mag(i)); y3 = log(mag(i + 1));
40+ den = y1 - 2 * y2 + y3;
41+ sh = 0;
42+ if abs(den) > 1e-12
43+ sh = min(0.5, max(-0.5, 0.5 * (y1 - y3) / den));
44+ end
45+ fc = fgrid(i) + sh * P.grid_hz;
46+ yc = exp(y2 - 0.25 * (y1 - y3) * sh);
47+ % local bandwidth from the two neighbouring grid points: for one section,
48+ % |Y(f_n + df)|^2 / |Y(f_n)|^2 = g^2 / ((4 pi df)^2 + g^2), so each
49+ % neighbour gives g = 4 pi |df| sqrt(r / (1 - r)). Staying local keeps the
50+ % estimate off the skirts of adjacent resonances, which a half-power scan
51+ % walks onto.
52+ gs = [];
53+ for s = [-1, 1]
54+ df = fgrid(i + s) - fc;
55+ r = (mag(i + s) / yc)^2;
56+ if r < 0.999
57+ gs = [gs, 4 * pi * abs(df) * sqrt(r / (1 - r))];
58+ end
59+ end
60+ if length(gs) == 0
61+ g = 2 * pi * P.grid_hz;
62+ else
63+ g = sum(gs) / length(gs);
64+ end
65+ fpk = [fpk, fc];
66+ qpk = [qpk, 2 * pi * fc / g];
67+ ypk = [ypk, yc];
68+end
69+
70+n = min(length(fpk), round(P.max_modes));
71+if n == 0
72+ error('No bore resonances below max_mode_hz -- check the geometry.');
73+end
74+fpk = fpk(1:n);
75+qpk = qpk(1:n);
76+ypk = ypk(1:n);
77+
78+% ---- aperture weights and standing-wave profiles at the mode centres ----
79+mo = recorder_bore(P, geo, fing, fpk);
80+
81+basis.n = n;
82+basis.note = note;
83+basis.f = reshape(fpk, n, 1);
84+basis.omega = 2 * pi * basis.f;
85+basis.Q = reshape(qpk, n, 1);
86+basis.g = basis.omega ./ basis.Q;
87+basis.B = basis.g .* reshape(ypk, n, 1);
88+basis.apX = mo.apX; % nap x 1 axial positions
89+basis.apW = mo.W; % nap x n complex flow weights
90+basis.profX = mo.profX; % 1 x nsta
91+basis.profP = abs(mo.profP); % nsta x n standing-wave magnitudes
92+basis.fTarget = geo.fTarget(note);
93+basis.Sm = geo.Sm;
94+basis.gridF = fgrid;
95+basis.gridY = mag;
96+end
matlab/recorder_bore.madded+178−0View file
@@ -0,0 +1,178 @@
1+function out = recorder_bore(P, geo, fing, fvec)
2+% RECORDER_BORE Transmission-line acoustics of the bore at given frequencies.
3+%
4+% The bore is a chain of short lossy cylindrical segments (the linear taper
5+% is discretized at 6 mm pitch) with the eight tone holes as lumped shunts:
6+% an open hole is the air mass of its chimney plus two end corrections and a
7+% baffled radiation resistance, a closed hole is the compliance of its
8+% chimney volume, and a pinched thumb hole is an open hole at a reduced
9+% effective radius. Viscothermal wall losses enter through the standard
10+% wide-tube propagation constant,
11+% Gamma = i k (1 + (1 - i) lossFac), lossFac = (dv / 2a)(1 + (gam-1)/sqrt(Pr)),
12+% with dv = sqrt(2 nu / omega) the viscous boundary layer. The foot is an
13+% unflanged open end; the window is an air mass (mouth end correction) with
14+% a baffled monopole resistance.
15+%
16+% Everything is vectorized over the frequency vector fvec (row, Hz).
17+%
18+% Returns, per frequency:
19+% out.Y series-loop admittance 1 / (Z_window + Z_bore): the passive
20+% resonator the jet drives, whose peaks are the playable notes
21+% out.Zbore bore input impedance at the window
22+% out.W (nap x nf) complex flow out of each radiating aperture per
23+% unit loop flow into the bore (aperture 1 is the window with
24+% weight -1, then each open hole, then the foot)
25+% out.apX (1 x nap) axial position of each radiating aperture, m
26+% out.profX (1 x nsta) axial stations along the bore
27+% out.profP (nsta x nf) complex pressure at each station per unit inflow
28+% (the standing-wave profile of each resonance when evaluated
29+% at its frequency)
30+
31+c = P.speed_of_sound;
32+rho0 = P.air_density;
33+nf = length(fvec);
34+om = 2 * pi * reshape(fvec, 1, nf);
35+k = om / c;
36+nu = 1.5e-5; % kinematic viscosity of air
37+gam = 1.4;
38+sqPr = sqrt(0.71);
39+dv = sqrt(2 * nu ./ max(om, 1)); % viscous boundary layer thickness
40+
41+L = geo.boreLen;
42+rH = geo.rHead;
43+rF = geo.rFoot;
44+
45+% ---- segment list: window (x = 0) to foot (x = L), split at the holes ----
46+bounds = [0, geo.holePos, L];
47+segX0 = [];
48+segX1 = [];
49+for j = 1:length(bounds) - 1
50+ span = bounds(j + 1) - bounds(j);
51+ if span <= 1e-9
52+ continue;
53+ end
54+ m = max(1, ceil(span / 0.006));
55+ e = bounds(j) + span * (0:m) / m;
56+ segX0 = [segX0, e(1:m)];
57+ segX1 = [segX1, e(2:m + 1)];
58+end
59+nseg = length(segX0);
60+
61+% ---- hole shunt impedances ----
62+% holeZ(j, :) for holes 1..8; holeOpen(j) marks a radiating (open) hole.
63+holeZ = zeros(8, nf) + 0i;
64+holeOpen = zeros(1, 8);
65+for j = 1:8
66+ b = geo.holeRad(j);
67+ t = geo.holeChim(j);
68+ s = fing(j);
69+ if s >= 0.25
70+ pinched = s < 0.75;
71+ if pinched
72+ b = b * P.thumb_pinch; % pinched: reduced effective radius
73+ end
74+ Sh = pi * b^2;
75+ teff = t + 1.5 * b; % inner + outer end corrections
76+ holeZ(j, :) = 1i * om * rho0 * teff / Sh + rho0 * c * (k * b).^2 / (2 * Sh);
77+ if pinched
78+ % The thumb crack is a narrow slit, not a round hole: its viscous
79+ % resistance is what damps the fundamental (whose pressure maximum
80+ % the vent sits near) while barely touching the second register,
81+ % which is the whole point of pinching.
82+ holeZ(j, :) = holeZ(j, :) + P.pinch_loss * rho0 * c / Sh;
83+ end
84+ holeOpen(j) = 1;
85+ else
86+ V = pi * b^2 * t;
87+ holeZ(j, :) = -1i * rho0 * c^2 ./ (max(om, 1) * V);
88+ end
89+end
90+
91+% ---- upward march: input impedance seen from the window ----
92+aF = rF;
93+SF = pi * aF^2;
94+ZcF = rho0 * c / SF;
95+Z = ZcF * ((k * aF).^2 / 4 + 1i * 0.61 * (k * aF)); % unflanged foot
96+
97+% precompute per-segment exp(Gamma l) and Zc, reused by both passes
98+segE = zeros(nseg, nf) + 0i;
99+segZc = zeros(nseg, 1);
100+for e = 1:nseg
101+ xm = (segX0(e) + segX1(e)) / 2;
102+ a = rH + (rF - rH) * xm / L;
103+ S = pi * a^2;
104+ segZc(e) = rho0 * c / S;
105+ lossFac = (dv / (2 * a)) * (1 + (gam - 1) / sqPr);
106+ G = 1i * k .* (1 + lossFac) + k .* lossFac;
107+ segE(e, :) = exp(G * (segX1(e) - segX0(e)));
108+end
109+
110+% which hole (if any) sits at the top boundary of each segment
111+holeAt = zeros(1, nseg);
112+for j = 1:8
113+ idx = find(abs(segX0 - geo.holePos(j)) < 1e-9);
114+ if length(idx) == 1
115+ holeAt(idx) = j;
116+ end
117+end
118+
119+for e = nseg:-1:1
120+ Zc = segZc(e);
121+ E = segE(e, :);
122+ E2 = E .* E;
123+ th = (E2 - 1) ./ (E2 + 1);
124+ Z = Zc * (Z + Zc * th) ./ (Zc + Z .* th);
125+ j = holeAt(e);
126+ if j > 0
127+ Zh = holeZ(j, :);
128+ Z = Z .* Zh ./ (Z + Zh);
129+ end
130+end
131+Zbore = Z;
132+
133+% ---- window impedance and the series-loop admittance ----
134+Zw = 1i * om * rho0 * geo.deltaM / geo.Sm + rho0 * c * k.^2 / (2 * pi);
135+Y = 1 ./ (Zw + Zbore);
136+
137+% ---- downward pass: pressure profile and aperture flows per unit inflow ----
138+nsta = nseg + 1;
139+profX = zeros(1, nsta);
140+profP = zeros(nsta, nf) + 0i;
141+p = Zbore; % pressure at the bore entry, U = 1
142+U = ones(1, nf) + 0i;
143+profX(1) = 0;
144+profP(1, :) = p;
145+apX = 0;
146+W = -(ones(1, nf) + 0i); % window: flow out = -inflow
147+for e = 1:nseg
148+ % at the top boundary of this segment, shunt flow leaves through a hole
149+ j = holeAt(e);
150+ if j > 0
151+ q = p ./ holeZ(j, :);
152+ U = U - q;
153+ if holeOpen(j)
154+ apX = [apX; segX0(e)];
155+ W = [W; q];
156+ end
157+ end
158+ Zc = segZc(e);
159+ E = segE(e, :);
160+ ch = (E + 1 ./ E) / 2;
161+ sh = (E - 1 ./ E) / 2;
162+ p2 = p .* ch - Zc * U .* sh;
163+ U = U .* ch - (p / Zc) .* sh;
164+ p = p2;
165+ profX(e + 1) = segX1(e);
166+ profP(e + 1, :) = p;
167+end
168+apX = [apX; L];
169+W = [W; U]; % what remains flows out the foot
170+
171+out.f = reshape(fvec, 1, nf);
172+out.Y = Y;
173+out.Zbore = Zbore;
174+out.W = W;
175+out.apX = reshape(apX, length(apX), 1);
176+out.profX = profX;
177+out.profP = profP;
178+end
matlab/recorder_defaults.madded+64−0View file
@@ -0,0 +1,64 @@
1+function P = recorder_defaults()
2+% RECORDER_DEFAULTS Default parameters, mirroring the web app's params.ts.
3+% A soprano (descant) recorder in C. The output duration is not a parameter:
4+% it is note_count * note_dur_s + (note_count - 1) * gap_ms + tail_s.
5+
6+% Bore (window at x = 0, open foot at x = bore_length_m)
7+P.bore_length_m = 0.3015;
8+P.bore_r_head_m = 0.0067;
9+P.bore_r_foot_m = 0.0052;
10+P.wall_m = 0.003;
11+P.hole_scale = 1.0;
12+P.thumb_pinch = 0.35;
13+P.pinch_loss = 0.5;
14+
15+% Window, windway, and jet
16+P.cutup_m = 0.0045; % flue exit to labium (window length along the bore)
17+P.window_width_m = 0.010;
18+P.flue_h_m = 0.0010; % flue exit height
19+P.flue_w_m = 0.010; % flue exit width
20+P.labium_offset_m = 0.0004;
21+P.jet_gain = 1.0; % receptivity trim
22+P.jet_growth = 0.4; % instability growth per flue height of travel
23+P.jet_conv = 0.55; % jet wave speed as a fraction of Uj
24+P.jet_stcut = 0.45; % jet instability cutoff, Strouhal-like: f_c = jet_stcut Uj / W
25+P.vortex_loss = 1.0; % nonlinear mouth loss coefficient
26+P.noise_level = 0.03; % turbulence displacement noise on the jet
27+P.direct_gain = 1.0;
28+
29+% Breath and gesture
30+P.breath_pa = 350.0; % blowing pressure for C5; higher notes track up
31+P.breath_track = 1.6; % pressure grows as (f_note / 523.25)^breath_track
32+P.attack_ms = 12.0;
33+P.release_ms = 25.0;
34+P.vibrato_hz = 5.0;
35+P.vibrato_depth = 0.0;
36+P.note_count = 1; % 1..3
37+P.note1 = 5; % index into the fingering chart: 5 = G5
38+P.note2 = 6; % A5
39+P.note3 = 7; % B5
40+P.note_dur_s = 1.6;
41+P.gap_ms = 60.0;
42+P.seed = 1;
43+
44+% Modal fit
45+P.max_mode_hz = 10000.0;
46+P.max_modes = 12;
47+P.grid_hz = 2.5;
48+
49+% Radiation / listener (window at the origin, foot along -y, holes toward +z)
50+P.listener_x_m = 0.25;
51+P.listener_y_m = -0.10;
52+P.listener_z_m = 0.35;
53+P.speed_of_sound = 343.0;
54+P.air_density = 1.204;
55+
56+% Rendering / output
57+P.sample_rate = 48000;
58+P.oversample = 2;
59+P.tail_s = 0.4;
60+P.pre_delay_ms = 3.0;
61+P.normalize = 1;
62+P.peak_dbfs = -3.0;
63+P.full_scale_pa = 0.5;
64+end
matlab/recorder_demo.madded+28−0View file
@@ -0,0 +1,28 @@
1+% Render a recorder note with default parameters and print diagnostics.
2+% Set P.note_count = 3 for a three-note figure; P.note1..note3 pick the
3+% fingerings (1 = C5 ... 10 = D6, 11..14 = pinched E6, F6, G6, A6).
4+% Run with: numbl run recorder_demo.m
5+
6+P = recorder_defaults();
7+tic;
8+[y, basis, dinfo] = recorder_render(P);
9+elapsed = toc;
10+
11+fprintf('bore modes : %d resonances, first at %.1f Hz (Q %.0f)\n', basis.n, basis.f(1), basis.Q(1));
12+fprintf('notes : %d\n', dinfo.noteCount);
13+for i = 1:dinfo.noteCount
14+ fprintf(' note %d : fingering %d, played %.1f Hz (target %.1f, %+.1f cents)\n', ...
15+ i, dinfo.noteIdx(i), dinfo.fPlay(i), 523.25 * 2^((dz_semi(dinfo.noteIdx(i))) / 12), dinfo.cents(i));
16+end
17+fprintf('peak press.: %.4f Pa (%.1f dB SPL) at %.0f cm\n', ...
18+ dinfo.peakPa, dinfo.peakSplDb, dinfo.listenerDistM * 100);
19+fprintf('length : %.2f s @ %d Hz\n', length(y) / P.sample_rate, P.sample_rate);
20+fprintf('render time: %.2f s\n', elapsed);
21+if dinfo.listenerInside
22+ fprintf('WARNING : listener is at or inside the instrument\n');
23+end
24+
25+function s = dz_semi(note)
26+semis = [0 2 4 5 7 9 11 12 13 14 16 17 19 21];
27+s = semis(note);
28+end
matlab/recorder_geometry.madded+63−0View file
@@ -0,0 +1,63 @@
1+function geo = recorder_geometry(P)
2+% RECORDER_GEOMETRY Bore, tone hole lattice, window, and fingering chart.
3+%
4+% The bore runs from the window (x = 0, where the jet drives it) to the open
5+% foot (x = bore_length_m), with a linear taper from r_head to r_foot; real
6+% recorder bores are stepped and locally reworked by the maker, so the taper
7+% here is a first-order stand-in. Eight tone holes: index 1 is the thumb
8+% hole on the back, 2..8 are finger holes 1..7 from the top down. Hole
9+% positions are stored as fractions of the design bore length, so scaling
10+% bore_length_m transposes the whole instrument; radii scale with
11+% P.hole_scale. The positions and radii below were tuned numerically (see
12+% the repository README) so that the first-register fingerings land on an
13+% equal-tempered C major scale at A4 = 440 with the default geometry.
14+%
15+% The fingering chart is baroque (English). Each row is one note, columns
16+% are the eight holes: 0 closed, 1 open, 0.5 pinched (the thumb nail crack
17+% that vents the second register; modeled as a reduced effective radius).
18+% Rows 1..10 are the first register C5..D6; rows 11..14 are pinched
19+% second-register fingerings E6..A6, whose intonation inherits every
20+% simplification of the model and is honestly worse.
21+
22+L = P.bore_length_m;
23+
24+% Tuned hole lattice (fractions of bore length from the window, and radii in
25+% m at the default geometry). Hole 1 is the thumb; 2..8 are finger holes.
26+geo.holeFrac = [0.3699 0.4271 0.4640 0.5292 0.6204 0.6519 0.7424 0.8278];
27+geo.holeRad = [0.0027 0.0032 0.0030 0.0029 0.0033 0.0026 0.0032 0.0025] * P.hole_scale;
28+geo.holePos = geo.holeFrac * L;
29+geo.holeChim = P.wall_m * ones(1, 8);
30+
31+% Fingering chart: thumb, then finger holes 1..7.
32+% T 1 2 3 4 5 6 7
33+geo.fingerings = [ ...
34+ 0 0 0 0 0 0 0 0; % 1 C5
35+ 0 0 0 0 0 0 0 1; % 2 D5
36+ 0 0 0 0 0 0 1 1; % 3 E5
37+ 0 0 0 0 0 1 0 0; % 4 F5 (fork)
38+ 0 0 0 0 1 1 1 1; % 5 G5
39+ 0 0 0 1 1 1 1 1; % 6 A5
40+ 0 0 1 1 1 1 1 1; % 7 B5
41+ 0 1 0 1 1 1 1 1; % 8 C6
42+ 1 0 0 1 1 1 1 1; % 9 C#6
43+ 1 1 0 1 1 1 1 1; % 10 D6
44+ 0.5 0 0 0 0 0 1 1; % 11 E6 (pinched E5)
45+ 0.5 0 0 0 0 1 0 1; % 12 F6 (pinched)
46+ 0.5 0 0 0 1 1 1 1; % 13 G6 (pinched G5)
47+ 0.5 0 0 1 1 1 1 1; % 14 A6 (pinched A5)
48+];
49+geo.nNotes = 14;
50+
51+% Equal-tempered targets, A4 = 440: semitones of each row above C5.
52+semis = [0 2 4 5 7 9 11 12 13 14 16 17 19 21];
53+geo.fTarget = 440 * 2^(3 / 12) * 2.^(reshape(semis, 1, 14) / 12);
54+
55+% Bore radius at position x (column or row vector), linear taper.
56+geo.rHead = P.bore_r_head_m;
57+geo.rFoot = P.bore_r_foot_m;
58+geo.boreLen = L;
59+
60+% Window (mouth) geometry.
61+geo.Sm = P.cutup_m * P.window_width_m; % window area
62+geo.deltaM = 0.7 * sqrt(geo.Sm); % mouth end correction
63+end
matlab/recorder_jet.madded+184−0View file
@@ -0,0 +1,184 @@
1+function [Zr, Zi, Qin, jinfo] = recorder_jet(P, basis, spanS, noteOrd, fsHi)
2+% RECORDER_JET Time-domain integration of the jet / bore feedback loop.
3+%
4+% This is the one stage of the model that cannot be closed-form: a recorder
5+% is a self-oscillator, so the loop is integrated sample by sample at fsHi.
6+% The loop, following the jet-drive model of Verge and Fabre:
7+%
8+% 1. Breath pressure p_f(t) (attack/release envelope, vibrato) sets the
9+% jet exit velocity Uj = sqrt(2 p_f / rho0) (Bernoulli through the flue).
10+% 2. The acoustic flow U through the mouth gives a transverse velocity
11+% v = U / S_m at the flue exit, which seeds a displacement wave on the
12+% jet. The wave convects at jet_conv Uj (so the transit delay varies
13+% with breath, which is what makes attacks speak the way they do),
14+% grows by exp(jet_growth W / h), and is low-passed (two poles) at
15+% jet_stcut Uj / W, the jet's own instability cutoff: frequencies the
16+% jet cannot amplify must not drive the loop, and a one-pole rolloff
17+% is shallow enough that an all-holes-closed fingering squeaks to its
18+% fifth resonance instead of speaking the fundamental. Turbulence
19+% enters here as seeded displacement noise scaling with Uj^1.5.
20+% 3. The displaced jet splits at the labium: the flow injected inside is
21+% Qin = w_j b Uj (1 + tanh((eta - y0) / b)), b the Bickley half-width
22+% grown over the cutup, y0 the labium offset.
23+% 4. The injection drives the bore as the jet-drive pressure source
24+% dp = (rho0 delta_d / S_m) dQin/dt. Combined with the modal admittance
25+% Y_n(s) = B_n s / (s^2 + g_n s + omega_n^2), the flow response of one
26+% mode to the injected jet flow is
27+% U_n(s) / Qin(s) = c_n s^2 / (s^2 + g_n s + omega_n^2),
28+% c_n = B_n rho0 delta_d / S_m,
29+% and this is discretized exactly as c_n (Qin - lowpass_n(Qin)) with
30+% one impulse-invariant complex pole per mode: no discrete derivative
31+% appears anywhere. (Approximating the two s factors by first
32+% differences at fsHi is Nyquist-unstable by three orders of magnitude,
33+% which is a property of the discretization, not of the physics.) The
34+% broadband c_n Qin part is the mouth spitting the injected flow
35+% straight back out through the window mass; the complex pole part is
36+% the resonant flow the apertures radiate mode by mode.
37+% 5. The nonlinear vortex loss at the window, a series resistance
38+% R = 0.5 rho0 K_v |v| / S_m, enters as amplitude-dependent extra
39+% damping B_n R on each mode; with the tanh it caps the amplitude.
40+% The total U closes the loop.
41+%
42+% Oscillation is not prescribed anywhere: it grows out of the noise when the
43+% loop gain crosses one, locks near a bore resonance where the phase closes
44+% (jet transit near a quarter period), and overblows to the second resonance
45+% when breath pressure carries the transit time out of the first's range.
46+%
47+% Returns the complex modal flow states Zr + i Zi (n x nT; the analytic
48+% signal 2 z_n of each modal flow, which recorder_radiate.m needs for its
49+% phase-correct aperture sum), the injected jet flow Qin (1 x nT, the direct
50+% broadband source), and diagnostics including the playing frequency
51+% estimated from zero crossings of the sustained half of the note.
52+
53+rho0 = P.air_density;
54+dt = 1 / fsHi;
55+nT = max(16, floor(spanS * fsHi));
56+n = basis.n;
57+
58+W = P.cutup_m;
59+h = P.flue_h_m;
60+Sm = basis.Sm;
61+bj = 0.41 * (h^2 * W)^(1 / 3); % Bickley jet half-width at the labium
62+growth = exp(min(2.5, P.jet_growth * W / h));
63+dScale = rho0 * (4 / pi) * sqrt(2 * h * W) / Sm;
64+lossC = 0.5 * rho0 * P.vortex_loss / Sm; % series resistance per unit |v|
65+y0 = P.labium_offset_m;
66+jg = P.jet_gain;
67+
68+% ---- breath pressure envelope and derived vectors (all precomputed) ----
69+t = (0:nT - 1) * dt;
70+tauA = max(1e-4, P.attack_ms * 1e-3) / 2.3;
71+tauR = max(1e-4, P.release_ms * 1e-3) / 2.3;
72+env = 1 - exp(-t / tauA);
73+rel = find(t > P.note_dur_s);
74+env(rel) = env(rel) .* exp(-(t(rel) - P.note_dur_s) / tauR);
75+vib = 1 + P.vibrato_depth * sin(2 * pi * P.vibrato_hz * t);
76+% Breath tracking: players blow higher notes harder. The jet transit time
77+% must stay near a quarter period of the note, so the pressure that centers
78+% the lock range grows as f^2; breath_track < 2 undershoots that slightly,
79+% which is what measured blowing pressures do.
80+trackFac = (basis.fTarget / 523.25)^P.breath_track;
81+pf = P.breath_pa * trackFac * env .* vib;
82+Uj = sqrt(2 * max(pf, 0) / rho0);
83+UjE = max(Uj, 1.0); % guarded: receptivity and cutoff need Uj > 0
84+tauS = min(0.005, W ./ (P.jet_conv * UjE)) * fsHi;
85+tauS = max(2, tauS);
86+aJet = 1 - exp(-2 * pi * (P.jet_stcut * UjE / W) * dt); % jet instability cutoff
87+etaAmp = P.noise_level * h * (UjE / 8).^1.5;
88+
89+% ---- seeded turbulence noise, filtered in the loop ----
90+state = mod(floor(abs(P.seed)) + 7919 * noteOrd, 4294967296);
91+[u, state] = recorder_rng(state, nT);
92+white = 2 * u.' - 1;
93+aN1 = 1 - exp(-2 * pi * 4000 * dt); % turbulence band: 400 Hz .. 4 kHz
94+aN2 = 1 - exp(-2 * pi * 400 * dt);
95+
96+% ---- per-mode discrete resonator coefficients (impulse invariant) ----
97+% U_n = c_n (Qin - 2 Re(rho_n z_n)), z_n' = lam_n z_n + Qin, with rho_n the
98+% residue of (g s + omega^2) / (s^2 + g s + omega^2) at lam_n.
99+omd = basis.omega .* sqrt(max(1e-12, 1 - 1 ./ (4 * basis.Q.^2)));
100+lam = -basis.g / 2 + 1i * omd;
101+Az = exp(lam * dt);
102+ar = real(Az);
103+ai = imag(Az);
104+cn = dScale * basis.B;
105+rho = (basis.g .* lam + basis.omega.^2) ./ (2i * omd);
106+% modal flow per state: U_res,n = wr_n zr_n + wi_n zi_n (the -2 c Re(rho z))
107+wr = -2 * cn .* real(rho);
108+wi = 2 * cn .* imag(rho);
109+cSum = sum(cn);
110+
111+% ---- the loop ----
112+Zr = zeros(n, nT);
113+Zi = zeros(n, nT);
114+Qin = zeros(1, nT);
115+vh = zeros(1, nT);
116+zr = zeros(n, 1);
117+zi = zeros(n, 1);
118+vprev = 0;
119+vlp = 0;
120+vlp2 = 0;
121+nlp1 = 0;
122+nlp2 = 0;
123+for k = 1:nT
124+ % turbulence: white noise bandpassed 400 Hz .. 4 kHz (two one-poles)
125+ nlp1 = nlp1 + aN1 * (white(k) - nlp1);
126+ nlp2 = nlp2 + aN2 * (nlp1 - nlp2);
127+ % transverse acoustic velocity at the flue exit, delayed by the jet
128+ % transit and low-passed at the jet's instability cutoff
129+ dk = k - tauS(k);
130+ vdel = 0;
131+ if dk >= 1
132+ i0 = floor(dk);
133+ wfr = dk - i0;
134+ vdel = vh(i0) * (1 - wfr) + vh(min(i0 + 1, nT)) * wfr;
135+ end
136+ vlp = vlp + aJet(k) * (vdel - vlp);
137+ vlp2 = vlp2 + aJet(k) * (vlp - vlp2);
138+ eta = jg * (h / UjE(k)) * growth * vlp2 + etaAmp(k) * (nlp1 - nlp2);
139+ qin = P.flue_w_m * bj * Uj(k) * (1 + tanh((eta - y0) / bj));
140+ damp = 1 - (dt / 2) * (lossC * abs(vprev)) * basis.B;
141+ zr2 = damp .* (ar .* zr - ai .* zi) + dt * qin;
142+ zi = damp .* (ar .* zi + ai .* zr);
143+ zr = zr2;
144+ Zr(:, k) = zr;
145+ Zi(:, k) = zi;
146+ U = cSum * qin + wr.' * zr + wi.' * zi;
147+ vprev = U / Sm;
148+ vh(k) = vprev;
149+ Qin(k) = qin;
150+end
151+
152+% ---- diagnostics: playing frequency from the sustained half ----
153+k0 = max(1, floor(0.35 * P.note_dur_s * fsHi));
154+k1 = min(nT, floor(0.95 * P.note_dur_s * fsHi));
155+seg = vh(k0:k1);
156+mu = sum(seg) / length(seg);
157+seg = seg - mu;
158+nzc = 0;
159+kFirst = 0;
160+kLast = 0;
161+for k = 2:length(seg)
162+ if seg(k - 1) <= 0 && seg(k) > 0
163+ nzc = nzc + 1;
164+ if kFirst == 0
165+ kFirst = k;
166+ end
167+ kLast = k;
168+ end
169+end
170+if nzc > 2
171+ jinfo.fPlay = (nzc - 1) * fsHi / (kLast - kFirst);
172+else
173+ jinfo.fPlay = 0;
174+end
175+jinfo.fTarget = basis.fTarget;
176+if jinfo.fPlay > 0
177+ jinfo.cents = 1200 * log(jinfo.fPlay / basis.fTarget) / log(2);
178+else
179+ jinfo.cents = 0;
180+end
181+jinfo.vacRms = sqrt(sum(seg.^2) / length(seg));
182+jinfo.ujMax = max(Uj);
183+jinfo.qinMean = sum(Qin) / nT;
184+end
matlab/recorder_notes.madded+34−0View file
@@ -0,0 +1,34 @@
1+function [noteIdx, tOn, spanS, ninfo] = recorder_notes(P)
2+% RECORDER_NOTES The gesture: one to three tongued notes.
3+%
4+% Each note starts note_dur_s + gap_ms after the previous one and is rendered
5+% over its own span: the sounding length plus enough ring for the release and
6+% the resonator decay (the last note gets the full tail_s). Notes are
7+% rendered independently and summed; the bore state does not carry across a
8+% note change, which is right for tongued (detached) articulation and is the
9+% model's stated limitation for legato.
10+%
11+% Returns row vectors: noteIdx (fingering chart rows), tOn (onset times, s),
12+% spanS (rendered span of each note, s), and ninfo with the total duration.
13+
14+K = min(3, max(1, round(P.note_count)));
15+notes = [P.note1, P.note2, P.note3];
16+gap = max(0, P.gap_ms) * 1e-3;
17+
18+noteIdx = zeros(1, K);
19+tOn = zeros(1, K);
20+spanS = zeros(1, K);
21+for i = 1:K
22+ noteIdx(i) = min(max(1, round(notes(i))), 14);
23+ tOn(i) = (i - 1) * (P.note_dur_s + gap);
24+ if i < K
25+ ring = min(P.tail_s, 0.30);
26+ else
27+ ring = P.tail_s;
28+ end
29+ spanS(i) = P.note_dur_s + ring;
30+end
31+
32+ninfo.count = K;
33+ninfo.durationS = tOn(K) + P.note_dur_s + P.tail_s;
34+end
matlab/recorder_output.madded+65−0View file
@@ -0,0 +1,65 @@
1+function [y, peakPa] = recorder_output(P, p, oversample)
2+% RECORDER_OUTPUT Anti-aliased decimation and output conditioning.
3+% Kaiser-windowed-sinc FIR lowpass at 0.45 * fs_lo, decimate by oversample,
4+% pre-roll, tail fade, then either normalize to peak_dbfs or divide by
5+% full_scale_pa. Identical in structure to the sibling instruments' output
6+% stages.
7+
8+peakPa = max(abs(p));
9+R = oversample;
10+T = length(p);
11+
12+if R > 1
13+ taps = 12 * R + 1;
14+ M2 = (taps - 1) / 2;
15+ fcNorm = 0.45 / R; % cycles per input sample
16+ m = (0:taps - 1) - M2;
17+ h = sin(2 * pi * fcNorm * m) ./ (pi * m);
18+ h(M2 + 1) = 2 * fcNorm; % sinc value at m = 0
19+ r = m / M2;
20+ beta = 8;
21+ h = h .* (rc_i0(beta * sqrt(max(0, 1 - r.^2))) / rc_i0(beta));
22+ h = h / sum(h); % unity DC gain
23+ pf = conv(p, h);
24+ outLen = floor(T / R);
25+ y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen); % centered pick
26+else
27+ y = p;
28+end
29+
30+% Pre-roll, then truncate to the nominal duration
31+padN = floor(P.pre_delay_ms * 1e-3 * P.sample_rate);
32+K = min(3, max(1, round(P.note_count)));
33+duration = K * P.note_dur_s + (K - 1) * max(0, P.gap_ms) * 1e-3 + P.tail_s;
34+total = floor(duration * P.sample_rate);
35+if padN > 0
36+ y = [zeros(1, padN) y];
37+end
38+y = y(1:min(length(y), total));
39+
40+% Squared fade on the tail
41+fade = min(2048, floor(length(y) / 8));
42+if fade > 1
43+ y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;
44+end
45+
46+if P.normalize ~= 0
47+ pk = max(abs(y));
48+ if pk > 0
49+ y = y / pk * 10^(P.peak_dbfs / 20);
50+ end
51+else
52+ y = y / P.full_scale_pa;
53+end
54+end
55+
56+function v = rc_i0(x)
57+% Modified Bessel function of the first kind, order zero (series expansion).
58+v = ones(size(x));
59+term = ones(size(x));
60+hx = x / 2;
61+for k = 1:40
62+ term = term .* (hx / k).^2;
63+ v = v + term;
64+end
65+end
matlab/recorder_radiate.madded+76−0View file
@@ -0,0 +1,76 @@
1+function [p, radinfo] = recorder_radiate(P, basis, Zr, Zi, Qin, fsHi)
2+% RECORDER_RADIATE Free-field pressure at the listener for one note.
3+%
4+% Two contributions:
5+%
6+% 1. Modal aperture radiation. The resonant part of each modal flow is the
7+% narrowband analytic signal a_n(t) = -2 c_n rho_n z_n(t) (see
8+% recorder_jet.m), and basis.apW holds each aperture's complex flow per
9+% unit loop flow at that mode's centre frequency. Every radiating
10+% aperture (the window, each open tone hole, the foot) is a point
11+% monopole, so per mode the apertures sum coherently with their relative
12+% retardations,
13+% G_n = sum_ap (i omega_n rho0 / (4 pi r_ap)) W_ap,n exp(-i k (r_ap - rcen)),
14+% and p(t) = sum_n Re(G_n a_n(t)). The common carrier exp(-i k rcen) is
15+% stripped: relative phase between the window, the open holes, and the
16+% foot is kept (this is where the directivity of a cross-fingered note
17+% lives), the overall propagation delay is not. Because each mode is
18+% narrowband, applying G_n at the centre frequency is the same
19+% approximation the sibling instruments make.
20+%
21+% 2. Direct broadband term. The injected jet flow minus the broadband part
22+% of the mouth response, (1 - sum c_n) Qin, leaves through the window and
23+% radiates as a monopole rho0 qdot / (4 pi r): the breath noise, the
24+% attack transient, and everything the modal truncation dropped. Scaled
25+% by direct_gain; this is the least principled term of the model, and
26+% direct_gain = 0 turns it off.
27+%
28+% The instrument: window at the origin, bore along -y (foot down), listener
29+% at (listener_x_m, listener_y_m, listener_z_m).
30+
31+rho0 = P.air_density;
32+c = P.speed_of_sound;
33+nT = size(Zr, 2);
34+n = basis.n;
35+
36+lx = P.listener_x_m;
37+ly = P.listener_y_m;
38+lz = P.listener_z_m;
39+
40+nap = length(basis.apX);
41+rap = zeros(nap, 1);
42+for a = 1:nap
43+ rap(a) = sqrt(lx^2 + (ly + basis.apX(a))^2 + lz^2);
44+end
45+listenerInside = min(rap) < 0.02;
46+rap = max(rap, 0.01);
47+rcen = rap(1); % the window
48+
49+% modal analytic-signal weights (must mirror recorder_jet.m)
50+dScale = rho0 * (4 / pi) * sqrt(2 * P.flue_h_m * P.cutup_m) / basis.Sm;
51+cn = dScale * basis.B;
52+omd = basis.omega .* sqrt(max(1e-12, 1 - 1 ./ (4 * basis.Q.^2)));
53+lam = -basis.g / 2 + 1i * omd;
54+rho = (basis.g .* lam + basis.omega.^2) ./ (2i * omd);
55+
56+H = zeros(1, n) + 0i;
57+for m = 1:n
58+ om = basis.omega(m);
59+ k = om / c;
60+ G = sum((1i * om * rho0 ./ (4 * pi * rap)) .* basis.apW(:, m) .* exp(-1i * k * (rap - rcen)));
61+ H(m) = G * (-2 * cn(m) * rho(m));
62+end
63+
64+p = real(H) * Zr - imag(H) * Zi; % 1 x nT
65+
66+% direct broadband monopole from the window
67+dq = zeros(1, nT);
68+if nT > 2
69+ dq(2:nT - 1) = (Qin(3:nT) - Qin(1:nT - 2)) * (fsHi / 2);
70+end
71+p = p + P.direct_gain * (1 - sum(cn)) * rho0 / (4 * pi * rcen) * dq;
72+
73+radinfo.gainAbs = abs(H);
74+radinfo.listenerInside = listenerInside;
75+radinfo.listenerDistM = rcen;
76+end
matlab/recorder_render.madded+49−0View file
@@ -0,0 +1,49 @@
1+function [y, basis, dinfo] = recorder_render(P)
2+% RECORDER_RENDER Full render: notes -> basis + jet -> radiate -> output.
3+% Convenience composition of the stage functions; the web app calls the
4+% stages separately so it can cache intermediates (the bore basis per note,
5+% and the jet integration, which is the expensive stage of this model).
6+% Returns the basis of the first note (the demo prints its resonances).
7+
8+[noteIdx, tOn, spanS, ninfo] = recorder_notes(P);
9+fsHi = P.sample_rate * P.oversample;
10+T = floor(ninfo.durationS * fsHi);
11+p = zeros(1, T);
12+
13+K = ninfo.count;
14+fPlay = zeros(1, K);
15+cents = zeros(1, K);
16+vacRms = zeros(1, K);
17+listenerInside = 0;
18+listenerDistM = 0;
19+basis = [];
20+for i = 1:K
21+ bi = recorder_basis(P, noteIdx(i));
22+ if i == 1
23+ basis = bi;
24+ end
25+ [Zr, Zi, Qin, jinfo] = recorder_jet(P, bi, spanS(i), i, fsHi);
26+ [pn, radinfo] = recorder_radiate(P, bi, Zr, Zi, Qin, fsHi);
27+ k0 = min(T, round(tOn(i) * fsHi) + 1);
28+ len = min(length(pn), T - k0 + 1);
29+ p(k0:k0 + len - 1) = p(k0:k0 + len - 1) + pn(1:len);
30+ fPlay(i) = jinfo.fPlay;
31+ cents(i) = jinfo.cents;
32+ vacRms(i) = jinfo.vacRms;
33+ listenerInside = listenerInside || radinfo.listenerInside;
34+ listenerDistM = radinfo.listenerDistM;
35+end
36+
37+[y, peakPa] = recorder_output(P, p, P.oversample);
38+
39+dinfo.noteCount = K;
40+dinfo.noteIdx = noteIdx;
41+dinfo.fPlay = fPlay;
42+dinfo.cents = cents;
43+dinfo.vacRms = vacRms;
44+dinfo.durationS = ninfo.durationS;
45+dinfo.peakPa = peakPa;
46+dinfo.peakSplDb = 20 * log10(peakPa / 20e-6 + 1e-30);
47+dinfo.listenerInside = listenerInside;
48+dinfo.listenerDistM = listenerDistM;
49+end
matlab/recorder_rng.madded+12−0View file
@@ -0,0 +1,12 @@
1+function [r, state] = recorder_rng(state, n)
2+% RECORDER_RNG Portable seeded uniform generator (32-bit LCG).
3+% Returns n uniforms in (0, 1) as a column vector, plus the advanced state.
4+% Deliberately not MATLAB's rand: the same seed must produce the identical
5+% breath noise in the numbl CLI, Node, and the browser, and every double
6+% here stays below 2^53 so the arithmetic is exact everywhere.
7+r = zeros(n, 1);
8+for i = 1:n
9+ state = mod(state * 1664525 + 1013904223, 4294967296);
10+ r(i) = (state + 0.5) / 4294967296;
11+end
12+end
matlab/recorder_validate.madded+167−0View file
@@ -0,0 +1,167 @@
1+% Validate the recorder model. There is no external reference implementation,
2+% so the checks are (a) closed-form identities the passive acoustics must
3+% satisfy, (b) internal identities of the modal fit, (c) emergent behavior
4+% the self-oscillating jet must reproduce (register selection, overblowing),
5+% (d) exact reproducibility of the seeded turbulence, and (e) pinned
6+% end-to-end diagnostics from the default configuration, recorded when the
7+% model was built and treated as a regression reference since.
8+% Run with: numbl run recorder_validate.m
9+
10+fails = 0;
11+P = recorder_defaults();
12+
13+% ---- 1. uniform open-open cylinder against the closed form ----
14+% With the tone holes shrunk to nothing and the taper removed, the bore
15+% admittance peaks must sit at n c / 2 (L + 0.61 a): the open-open pipe with
16+% the unflanged end correction at the foot (the window end is where Zbore is
17+% measured, so it carries no correction of its own). Viscothermal losses
18+% shift these by well under a percent.
19+Q = P;
20+Q.bore_r_head_m = 0.006;
21+Q.bore_r_foot_m = 0.006;
22+geo = recorder_geometry(Q);
23+geo.holeRad = 1e-9 * ones(1, 8);
24+f = 100:1:2500;
25+out = recorder_bore(Q, geo, zeros(1, 8), f);
26+mag = abs(1 ./ out.Zbore);
27+pk = [];
28+for i = 2:length(f) - 1
29+ if mag(i) > mag(i - 1) && mag(i) >= mag(i + 1)
30+ pk = [pk, f(i)];
31+ end
32+end
33+c = Q.speed_of_sound;
34+Leff = geo.boreLen + 0.61 * 0.006;
35+fClosed = (1:3) * c / (2 * Leff);
36+% viscothermal dispersion: the phase velocity is c / (1 + lossFac), with
37+% lossFac = (dv / 2a)(1 + (gam - 1) / sqrt(Pr)) evaluated at each resonance
38+dv = sqrt(2 * 1.5e-5 ./ (2 * pi * fClosed));
39+fClosed = fClosed ./ (1 + (dv / 0.012) * (1 + 0.4 / sqrt(0.71)));
40+errC = max(abs(pk(1:3) - fClosed) ./ fClosed);
41+ok = errC < 3e-3;
42+fails = fails + (~ok);
43+fprintf('open-open cylinder: %.1f %.1f %.1f Hz vs closed form %.1f %.1f %.1f (err %.1e) %s\n', ...
44+ pk(1), pk(2), pk(3), fClosed(1), fClosed(2), fClosed(3), errC, rv_verdict(ok));
45+
46+% ---- 2. the modal fit reconstructs the admittance it was fitted to ----
47+% At each fitted centre frequency the modal sum must reproduce the computed
48+% series-loop admittance. The check covers the first four resonances, which
49+% carry the oscillation; above them the per-peak fit overcounts by the
50+% skirts of neighbouring modes (up to ~20% on the weakest peaks), which is
51+% inherent to independent-peak fitting and acoustically negligible there.
52+basis = recorder_basis(P, 5);
53+geo = recorder_geometry(P);
54+bo = recorder_bore(P, geo, geo.fingerings(5, :), reshape(basis.f, 1, basis.n));
55+eFit = 0;
56+for m = 1:min(4, basis.n)
57+ om = basis.omega(m);
58+ Yfit = 0;
59+ for j = 1:basis.n
60+ Yfit = Yfit + 1i * om * basis.B(j) / (basis.omega(j)^2 - om^2 + 1i * om * basis.g(j));
61+ end
62+ eFit = max(eFit, abs(abs(Yfit) - abs(bo.Y(m))) / abs(bo.Y(m)));
63+end
64+ok = eFit < 0.06;
65+fails = fails + (~ok);
66+fprintf('modal fit: worst |Y| error over the first four peaks %.1e %s\n', eFit, rv_verdict(ok));
67+
68+% ---- 3. register selection is emergent ----
69+% Three fingerings that must each speak their own register at the default
70+% (tracked) breath: C5 all closed (the easiest fingering to squeak, so this
71+% is the sharpest test of the jet's instability cutoff), G5, and pinched E6,
72+% whose fundamental must be suppressed by the vent and not played.
73+Q = P;
74+Q.note_dur_s = 0.6;
75+Q.tail_s = 0.2;
76+fsHi = P.sample_rate * P.oversample;
77+centsAbs = zeros(1, 3);
78+notes3 = [1 5 11];
79+for i = 1:3
80+ b = recorder_basis(Q, notes3(i));
81+ [~, ~, ~, j] = recorder_jet(Q, b, 0.7, 1, fsHi);
82+ centsAbs(i) = abs(j.cents);
83+end
84+ok = max(centsAbs) < 60;
85+fails = fails + (~ok);
86+fprintf('registers: C5 %+.0f c, G5 %+.0f c, E6 %+.0f c of equal temperament %s\n', ...
87+ centsAbs(1), centsAbs(2), centsAbs(3), rv_verdict(ok));
88+
89+% ---- 4. overblowing is emergent ----
90+% Quadrupling the breath pressure on the G5 fingering must carry the jet
91+% transit time out of the fundamental's lock range and onto the second
92+% resonance: an octave (plus the resonator's inharmonicity) up. Nothing in
93+% the model prescribes this.
94+Q = P;
95+Q.note_dur_s = 0.6;
96+Q.tail_s = 0.2;
97+Q.breath_pa = 4 * P.breath_pa;
98+b = recorder_basis(Q, 5);
99+[~, ~, ~, j] = recorder_jet(Q, b, 0.7, 1, fsHi);
100+ratio = j.fPlay / b.f(1);
101+ok = ratio > 1.9 && ratio < 2.1;
102+fails = fails + (~ok);
103+fprintf('overblow at 4x breath: %.1f Hz = %.3f x f1 %s\n', j.fPlay, ratio, rv_verdict(ok));
104+
105+% ---- 5. the seeded turbulence is exactly reproducible ----
106+Q = P;
107+Q.note_dur_s = 0.4;
108+Q.tail_s = 0.15;
109+[y1, ~, ~] = recorder_render(Q);
110+[y2, ~, ~] = recorder_render(Q);
111+Q.seed = 2;
112+[y3, ~, ~] = recorder_render(Q);
113+same = max(abs(y1 - y2)) == 0;
114+diff = max(abs(y1 - y3)) > 0;
115+ok = same && diff;
116+fails = fails + (~ok);
117+fprintf('seed: identical rerun %d, different seed differs %d %s\n', same, diff, rv_verdict(ok));
118+
119+% ---- 6. far-field falloff: -6 dB per doubling of distance ----
120+% RMS over the render, both listeners on one ray from the window so every
121+% aperture monopole keeps its interference geometry.
122+Q = P;
123+Q.normalize = 0;
124+Q.note_dur_s = 0.5;
125+Q.tail_s = 0.15;
126+Q.listener_x_m = 0.5; Q.listener_y_m = 0.4; Q.listener_z_m = 1.0;
127+[y1, ~, ~] = recorder_render(Q);
128+Q.listener_x_m = 1.0; Q.listener_y_m = 0.8; Q.listener_z_m = 2.0;
129+[y2, ~, ~] = recorder_render(Q);
130+drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
131+ok = abs(drop + 6.02) < 0.35;
132+fails = fails + (~ok);
133+fprintf('falloff 1.2 m -> 2.4 m: %.2f dB rms (target -6.02) %s\n', drop, rv_verdict(ok));
134+
135+% ---- 7. pinned end-to-end diagnostics at the defaults ----
136+[y, basis, dinfo] = recorder_render(P);
137+ok = abs(dinfo.fPlay(1) - 781.88) < 1.5 ...
138+ && abs(dinfo.peakPa - 1.7982) < 0.05 ...
139+ && basis.n == 12 ...
140+ && abs(basis.f(1) - 783.94) < 0.1;
141+fails = fails + (~ok);
142+fprintf('pinned render: played %.2f Hz (target %.2f), peak %.4f Pa (target %.4f), %d modes %s\n', ...
143+ dinfo.fPlay(1), 781.88, dinfo.peakPa, 1.7982, basis.n, rv_verdict(ok));
144+
145+% ---- 8. the note stops when the breath stops ----
146+fs = P.sample_rate;
147+n = length(y);
148+sus = y(round(0.5 * fs):round(1.0 * fs));
149+tail = y(n - round(0.2 * fs):n);
150+dropE = 10 * log10(sum(tail.^2) / length(tail) / (sum(sus.^2) / length(sus)));
151+ok = dropE < -40;
152+fails = fails + (~ok);
153+fprintf('release: last 200 ms %.1f dB re sustain %s\n', dropE, rv_verdict(ok));
154+
155+if fails > 0
156+ fprintf('%d CHECK(S) FAILED\n', fails);
157+else
158+ fprintf('all checks passed\n');
159+end
160+
161+function s = rv_verdict(ok)
162+if ok
163+ s = 'ok';
164+else
165+ s = 'MISMATCH';
166+end
167+end
webapp/.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
webapp/README.mdadded+42−0View file
@@ -0,0 +1,42 @@
1+# Recorder Simulator — web app
2+
3+Interactive browser implementation of the physical model described in `../README.md`. The algorithm itself is written in MATLAB syntax, in `../matlab/*.m`, and executed client-side by [numbl](https://numbl.org) inside a Web Worker; the TypeScript here is only the UI, the plots, the WAV encoder, and the marshalling of parameters and tensors in and out of the numbl runtime.
4+
5+The signal chain: a transmission-line solve of the tapered bore with its tone hole lattice under the chosen fingering, a modal fit of the resonances, a sample-by-sample integration of the jet / bore feedback loop (the note is a self-sustained oscillation, not a prescribed one), a coherent sum of aperture monopoles to the listener point, and Kaiser-windowed anti-aliased decimation to 48 kHz.
6+
7+## Usage
8+
9+```
10+npm install
11+npm run dev # development server
12+npm run build # production build in dist/
13+npm run selftest # run ../matlab/recorder_validate.m through numbl in Node
14+```
15+
16+The self-test runs the closed-form and reproducibility checks, the emergent register and overblow checks, and the pinned diagnostics. The MATLAB code can also be run directly: `cd ../matlab && numbl run recorder_demo.m` (or `recorder_validate.m`).
17+
18+## Why rendering is tiered the way it is
19+
20+Unlike the sibling instruments, nothing here is closed-form between events: a recorder is a self-oscillator, so the jet loop is integrated sample by sample and a two-second note costs a few seconds. The worker therefore runs each stage as a separate numbl execution and keys its caches by tier: bore, hole, and window parameters invalidate the per-note basis (the transmission-line solve and modal fit); the breath, the jet, and the gesture re-run only the jet integration; the listener position and the direct gain re-run only the radiation, which is a matrix multiply; output conditioning is always re-run. Intermediate tensors stay inside the worker as numbl RuntimeValues. Changing any parameter triggers a debounced fast preview (oversample 1); the Render button runs the full-quality pass, and while the current result is a preview the Play and Download buttons and the status badge say so.
21+
22+The New seed button redraws the breath turbulence without touching any other parameter. The Loop button repeats the render continuously, swapping in each new render as parameters change. Parameters persist in localStorage; "Reset to defaults" restores the values in `src/params.ts` (which `../matlab/recorder_defaults.m` mirrors).
23+
24+## Reading and editing the code in the app
25+
26+The Code button opens the MATLAB source in a CodeMirror editor with MATLAB syntax highlighting. Every `.m` file in `../matlab` is listed, tagged by whether the app executes it or it is a CLI entry point. The sources are editable: an edit re-renders after the usual debounce, so changing the physics in the editor is heard and plotted the same way moving a slider is. Doubling `jet_conv`'s effect by halving the transit delay in `recorder_jet.m`, for instance, audibly sharpens every note and then flips the register.
27+
28+Edits live in memory only. They are deliberately not persisted, so a reload restores the code the app shipped with; "Revert file" and "Revert all" undo them sooner, and a dot marks the files that differ from the shipped source. A fault in edited code is reported with its file and line (numbl skips a file it cannot parse and then fails at the call site, so the worker captures the parse warning and reports that instead), the last successful render stays on screen, and reverting recovers.
29+
30+## Layout
31+
32+- `../matlab/` — the algorithm: `recorder_geometry.m`, `recorder_bore.m`, `recorder_basis.m`, `recorder_notes.m`, `recorder_jet.m`, `recorder_radiate.m`, `recorder_output.m`, `recorder_rng.m`, plus `recorder_render.m` / `recorder_demo.m` / `recorder_validate.m` for CLI use
33+- `src/params.ts` — parameter schema, defaults, tier assignment, persistence
34+- `src/mfiles.ts` — the `.m` sources, bundled at build time (`src/mfilesCore.ts` holds the hash used for cache keying)
35+- `src/worker/renderWorker.ts` — off-thread numbl execution with tiered caching
36+- `src/lib/instrument.ts` — the hole lattice and fingering chart, mirrored from `recorder_geometry.m` for the 3D scene
37+- `src/components/` — 3D scene (three.js), parameter panel, code panel and editor, plots (jet flow, waveform, spectrum, spectrogram, bore admittance, standing waves)
38+- `src/selftest.ts` — Node harness for the MATLAB validation script
39+
40+## A note on the seeded generator
41+
42+The turbulence noise draws from a 32-bit LCG in `recorder_rng.m` rather than MATLAB's `rand`, so a given seed produces the identical breath in the numbl CLI, in Node, and in the browser. That is what lets the self-test require two renders to be bit-identical.
webapp/index.htmladded+13−0View file
@@ -0,0 +1,13 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
6+ <link rel="icon" href="data:image/svg+xml,<svg xmlns=%22http://www.w3.org/2000/svg%22 viewBox=%220 0 100 100%22><text y=%22.9em%22 font-size=%2290%22>🪈</text></svg>" />
7+ <title>Recorder Simulator</title>
8+ </head>
9+ <body>
10+ <div id="root"></div>
11+ <script type="module" src="/src/main.tsx"></script>
12+ </body>
13+</html>
webapp/package-lock.jsonadded+4230−0View file
This diff is 4,235 lines long and is not shown.
webapp/package.jsonadded+34−0View file
@@ -0,0 +1,34 @@
1+{
2+ "name": "recorder-webapp",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "scripts": {
7+ "dev": "vite",
8+ "build": "tsc --noEmit && vite build",
9+ "preview": "vite preview",
10+ "selftest": "tsx src/selftest.ts"
11+ },
12+ "dependencies": {
13+ "@codemirror/commands": "^6.11.0",
14+ "@codemirror/language": "^6.12.4",
15+ "@codemirror/legacy-modes": "^6.5.3",
16+ "@codemirror/state": "^6.7.1",
17+ "@codemirror/theme-one-dark": "^6.1.3",
18+ "@codemirror/view": "^6.43.9",
19+ "codemirror": "^6.0.2",
20+ "numbl": "^0.4.19",
21+ "react": "^18.3.1",
22+ "react-dom": "^18.3.1",
23+ "three": "^0.168.0"
24+ },
25+ "devDependencies": {
26+ "@types/react": "^18.3.5",
27+ "@types/react-dom": "^18.3.0",
28+ "@types/three": "^0.168.0",
29+ "@vitejs/plugin-react": "^4.3.1",
30+ "tsx": "^4.19.0",
31+ "typescript": "^5.5.4",
32+ "vite": "^5.4.3"
33+ }
34+}
webapp/src/App.tsxadded+271−0View file
@@ -0,0 +1,271 @@
1+import { useCallback, useEffect, useRef, useState } from 'react';
2+import { DEFAULTS, loadParams, saveParams, type RecorderParams } from './params';
3+import { DEFAULT_MFILES, type MFile } from './mfiles';
4+import { useRenderer } from './hooks/useRenderer';
5+import { encodeWav } from './lib/wav';
6+import ParamPanel from './components/ParamPanel';
7+import CodePanel from './components/CodePanel';
8+import Scene3D from './components/Scene3D';
9+import DiagnosticsPanel from './components/DiagnosticsPanel';
10+import WaveformPlot from './components/plots/WaveformPlot';
11+import SpectrumPlot from './components/plots/SpectrumPlot';
12+import SpectrogramPlot from './components/plots/SpectrogramPlot';
13+import BreathPlot from './components/plots/BreathPlot';
14+import AdmittancePlot from './components/plots/AdmittancePlot';
15+import BoreProfile from './components/plots/BoreProfile';
16+
17+export default function App() {
18+ const [params, setParams] = useState<RecorderParams>(() => loadParams());
19+ // Editable copy of the MATLAB implementation. Deliberately not persisted:
20+ // a reload restores the code the app shipped with.
21+ const [mfiles, setMfiles] = useState<MFile[]>(() => DEFAULT_MFILES.map((f) => ({ ...f })));
22+ const [codeOpen, setCodeOpen] = useState(false);
23+ const { result, error, busy, stage, progress, renderFull } = useRenderer(params, mfiles);
24+ const audioCtxRef = useRef<AudioContext | null>(null);
25+
26+ useEffect(() => saveParams(params), [params]);
27+
28+ const onChange = useCallback((key: keyof RecorderParams, value: number | boolean) => {
29+ setParams((p) => ({ ...p, [key]: value }));
30+ }, []);
31+
32+ const reset = useCallback(() => setParams({ ...DEFAULTS }), []);
33+
34+ // A new seed redraws the breath turbulence without touching anything else.
35+ const reseed = useCallback(
36+ () => setParams((p) => ({ ...p, seed: 1 + Math.floor(Math.random() * 9999) })),
37+ [],
38+ );
39+
40+ const onChangeFile = useCallback((name: string, source: string) => {
41+ setMfiles((fs) => fs.map((f) => (f.name === name ? { ...f, source } : f)));
42+ }, []);
43+
44+ const onResetFile = useCallback((name: string) => {
45+ setMfiles((fs) =>
46+ fs.map((f) =>
47+ f.name === name
48+ ? { ...f, source: DEFAULT_MFILES.find((d) => d.name === name)?.source ?? f.source }
49+ : f,
50+ ),
51+ );
52+ }, []);
53+
54+ const onResetAllFiles = useCallback(() => setMfiles(DEFAULT_MFILES.map((f) => ({ ...f }))), []);
55+ const codeEdited = mfiles.some(
56+ (f) => f.source !== DEFAULT_MFILES.find((d) => d.name === f.name)?.source,
57+ );
58+
59+ const [looping, setLooping] = useState(false);
60+ const loopSrcRef = useRef<AudioBufferSourceNode | null>(null);
61+
62+ const getAudioCtx = useCallback((sampleRate: number) => {
63+ let ctx = audioCtxRef.current;
64+ if (!ctx || ctx.state === 'closed') {
65+ ctx = new AudioContext({ sampleRate });
66+ audioCtxRef.current = ctx;
67+ }
68+ void ctx.resume();
69+ return ctx;
70+ }, []);
71+
72+ const makeBuffer = useCallback((ctx: AudioContext, r: NonNullable<typeof result>) => {
73+ const buf = ctx.createBuffer(1, r.y.length, r.sampleRate);
74+ buf.copyToChannel(new Float32Array(r.y), 0);
75+ return buf;
76+ }, []);
77+
78+ const play = useCallback(() => {
79+ if (!result) return;
80+ const ctx = getAudioCtx(result.sampleRate);
81+ const src = ctx.createBufferSource();
82+ src.buffer = makeBuffer(ctx, result);
83+ src.connect(ctx.destination);
84+ src.start();
85+ }, [result, getAudioCtx, makeBuffer]);
86+
87+ // Loop mode: repeat the note or figure continuously, swapping in each new
88+ // render as parameters change so edits are audible without pressing Play.
89+ useEffect(() => {
90+ if (!looping || !result) return;
91+ const ctx = getAudioCtx(result.sampleRate);
92+ const src = ctx.createBufferSource();
93+ src.buffer = makeBuffer(ctx, result);
94+ src.loop = true;
95+ src.connect(ctx.destination);
96+ src.start();
97+ loopSrcRef.current = src;
98+ return () => {
99+ loopSrcRef.current = null;
100+ try {
101+ src.stop();
102+ } catch {
103+ // already stopped
104+ }
105+ src.disconnect();
106+ };
107+ }, [looping, result, getAudioCtx, makeBuffer]);
108+
109+ const download = useCallback(() => {
110+ if (!result) return;
111+ const blob = encodeWav(result.y, result.sampleRate, params.bit_depth === 16 ? 16 : 24);
112+ const url = URL.createObjectURL(blob);
113+ const a = document.createElement('a');
114+ a.href = url;
115+ a.download = result.quality === 'full' ? 'recorder.wav' : 'recorder_preview.wav';
116+ a.click();
117+ setTimeout(() => URL.revokeObjectURL(url), 5000);
118+ }, [result, params.bit_depth]);
119+
120+ return (
121+ <div className="app">
122+ <div className="topbar">
123+ <h1>Recorder Simulator</h1>
124+ <button className="primary" onClick={renderFull} disabled={busy}>
125+ Render
126+ </button>
127+ <button
128+ className={`play ${result?.quality === 'preview' ? 'preview' : ''}`}
129+ onClick={play}
130+ disabled={!result}
131+ title={
132+ result?.quality === 'preview'
133+ ? 'Playing the fast preview render; press Render for full quality'
134+ : 'Play the rendered note or figure'
135+ }
136+ >
137+ ▶ Play{result?.quality === 'preview' ? ' (preview)' : ''}
138+ </button>
139+ <button
140+ className={`loop ${looping ? 'active' : ''}`}
141+ onClick={() => setLooping((v) => !v)}
142+ disabled={!result && !looping}
143+ title="Repeat the render continuously so parameter changes are audible immediately"
144+ >
145+ ⟳ Loop{looping ? ' on' : ''}
146+ </button>
147+ <button onClick={reseed} title="Draw new breath turbulence from a new random seed">
148+ ⚄ New seed
149+ </button>
150+ <button onClick={download} disabled={!result}>
151+ Download .wav{result?.quality === 'preview' ? ' (preview)' : ''}
152+ </button>
153+ <button onClick={reset}>Reset to defaults</button>
154+ <button
155+ className={codeOpen ? 'code active' : 'code'}
156+ onClick={() => setCodeOpen((v) => !v)}
157+ title="View and edit the MATLAB source that numbl executes"
158+ >
159+ {'{ }'} Code{codeEdited ? ' •' : ''}
160+ </button>
161+ <div className="status">
162+ {error && <span className="badge err">code error</span>}
163+ {result && (
164+ <span className={`badge ${result.quality}`}>
165+ {result.quality === 'full' ? 'full quality' : 'preview quality — press Render'}
166+ </span>
167+ )}
168+ {busy ? (
169+ <>
170+ <span>{stage}…</span>
171+ <div className="progressbar">
172+ <div style={{ width: `${Math.round(progress * 100)}%` }} />
173+ </div>
174+ </>
175+ ) : (
176+ <span>ready</span>
177+ )}
178+ </div>
179+ </div>
180+ <div className="main">
181+ <div className="sidebar">
182+ <ParamPanel params={params} onChange={onChange} />
183+ </div>
184+ <div className="content">
185+ {error && (
186+ <div className="card wide error-card">
187+ <h2>MATLAB error</h2>
188+ <div>
189+ <strong>
190+ {error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}
191+ </strong>{' '}
192+ {error.message}
193+ </div>
194+ <div className="code-hint">
195+ Showing the last successful render. Open Code to fix the source, or revert it.
196+ </div>
197+ </div>
198+ )}
199+ <div className="card wide">
200+ <h2>Geometry — bore, tone holes, microphone</h2>
201+ <Scene3D params={params} />
202+ </div>
203+ {result && (
204+ <>
205+ <div className="card">
206+ <h2>Diagnostics</h2>
207+ <DiagnosticsPanel result={result} />
208+ </div>
209+ <div className="card">
210+ <h2>Jet flow</h2>
211+ <BreathPlot
212+ qPlot={result.qPlot}
213+ qPlotRate={result.qPlotRate}
214+ noteOn={result.noteOn}
215+ fPlay={result.fPlay}
216+ durationS={result.diag.durationS}
217+ />
218+ </div>
219+ <div className="card">
220+ <h2>Waveform</h2>
221+ <WaveformPlot y={result.y} sampleRate={result.sampleRate} />
222+ </div>
223+ <div className="card">
224+ <h2>Spectrum</h2>
225+ <SpectrumPlot
226+ y={result.y}
227+ sampleRate={result.sampleRate}
228+ modeF={result.modeF}
229+ fPlay={result.fPlay}
230+ />
231+ </div>
232+ <div className="card">
233+ <h2>Spectrogram</h2>
234+ <SpectrogramPlot y={result.y} sampleRate={result.sampleRate} />
235+ </div>
236+ <div className="card">
237+ <h2>Bore admittance</h2>
238+ <AdmittancePlot
239+ gridF={result.gridF}
240+ gridY={result.gridY}
241+ modeF={result.modeF}
242+ fPlay={result.fPlay.length > 0 ? result.fPlay[0] : 0}
243+ />
244+ </div>
245+ <div className="card wide">
246+ <h2>Standing waves</h2>
247+ <BoreProfile
248+ profX={result.profX}
249+ profP={result.profP}
250+ modeF={result.modeF}
251+ />
252+ </div>
253+ </>
254+ )}
255+ </div>
256+ </div>
257+ {codeOpen && (
258+ <CodePanel
259+ mfiles={mfiles}
260+ onChangeFile={onChangeFile}
261+ onResetFile={onResetFile}
262+ onResetAll={onResetAllFiles}
263+ onClose={() => setCodeOpen(false)}
264+ onRenderFull={renderFull}
265+ error={error}
266+ busy={busy}
267+ />
268+ )}
269+ </div>
270+ );
271+}
webapp/src/components/CodeEditor.tsxadded+91−0View file
@@ -0,0 +1,91 @@
1+// CodeMirror 6 editor for one .m file. The parent gives each instance a key
2+// of the file name, so switching files remounts rather than syncing docs.
3+
4+import { useEffect, useRef } from 'react';
5+import { EditorState } from '@codemirror/state';
6+import { EditorView, keymap, lineNumbers, highlightActiveLine } from '@codemirror/view';
7+import { defaultKeymap, history, historyKeymap, indentWithTab } from '@codemirror/commands';
8+import { StreamLanguage, syntaxHighlighting, HighlightStyle } from '@codemirror/language';
9+import type { StreamParser } from '@codemirror/language';
10+import { octave } from '@codemirror/legacy-modes/mode/octave';
11+import { tags as t } from '@lezer/highlight';
12+
13+// The octave mode's keyword list is missing these MATLAB keywords, so promote
14+// them after it has tokenized the word.
15+const EXTRA_KEYWORDS = new Set(['function', 'break', 'parfor', 'spmd']);
16+
17+const matlabMode: StreamParser<unknown> = {
18+ ...octave,
19+ token(stream, state) {
20+ const start = stream.pos;
21+ const tok = octave.token(stream, state);
22+ if (EXTRA_KEYWORDS.has(stream.string.slice(start, stream.pos))) return 'keyword';
23+ return tok;
24+ },
25+};
26+
27+const matlabHighlight = HighlightStyle.define([
28+ { tag: t.keyword, color: '#e8a33d' },
29+ { tag: t.comment, color: '#6b7a86', fontStyle: 'italic' },
30+ { tag: t.string, color: '#a5d6a7' },
31+ { tag: t.number, color: '#4dd0e1' },
32+ { tag: t.operator, color: '#c9a0dc' },
33+ { tag: t.variableName, color: '#d6dae2' },
34+ { tag: t.standard(t.variableName), color: '#82b1ff' }, // builtins: zeros, sqrt, ...
35+ { tag: t.function(t.variableName), color: '#82b1ff' },
36+ { tag: t.atom, color: '#4dd0e1' },
37+ { tag: t.meta, color: '#c9a0dc' },
38+ { tag: t.invalid, color: '#e57373' },
39+]);
40+
41+const editorTheme = EditorView.theme(
42+ {
43+ '&': { color: '#d6dae2', backgroundColor: '#181b20', height: '100%' },
44+ '.cm-content': { fontFamily: 'ui-monospace, SFMono-Regular, Menlo, monospace', fontSize: '12.5px' },
45+ '.cm-gutters': { backgroundColor: '#181b20', color: '#4d5560', border: 'none' },
46+ '.cm-activeLine': { backgroundColor: '#20242b' },
47+ '.cm-cursor': { borderLeftColor: '#e8a33d' },
48+ '&.cm-focused': { outline: 'none' },
49+ '.cm-scroller': { overflow: 'auto' },
50+ '.cm-selectionBackground, ::selection': { backgroundColor: '#31465a' },
51+ '&.cm-focused .cm-selectionBackground': { backgroundColor: '#31465a' },
52+ },
53+ { dark: true },
54+);
55+
56+interface Props {
57+ source: string;
58+ onChange: (source: string) => void;
59+}
60+
61+export default function CodeEditor({ source, onChange }: Props) {
62+ const hostRef = useRef<HTMLDivElement | null>(null);
63+ const onChangeRef = useRef(onChange);
64+ onChangeRef.current = onChange;
65+
66+ useEffect(() => {
67+ const host = hostRef.current!;
68+ const view = new EditorView({
69+ parent: host,
70+ state: EditorState.create({
71+ doc: source,
72+ extensions: [
73+ lineNumbers(),
74+ highlightActiveLine(),
75+ history(),
76+ keymap.of([...defaultKeymap, ...historyKeymap, indentWithTab]),
77+ StreamLanguage.define(matlabMode),
78+ syntaxHighlighting(matlabHighlight),
79+ editorTheme,
80+ EditorView.updateListener.of((u) => {
81+ if (u.docChanged) onChangeRef.current(u.state.doc.toString());
82+ }),
83+ ],
84+ }),
85+ });
86+ return () => view.destroy();
87+ // eslint-disable-next-line react-hooks/exhaustive-deps
88+ }, []);
89+
90+ return <div className="code-editor" ref={hostRef} />;
91+}
webapp/src/components/CodePanel.tsxadded+109−0View file
@@ -0,0 +1,109 @@
1+// Overlay for reading and editing the MATLAB implementation. Edits take
2+// effect on the next debounced render and live only as long as the page.
3+
4+import { useState } from 'react';
5+import { DEFAULT_MFILES, EXECUTED_BY_APP, type MFile } from '../mfiles';
6+import type { RenderError } from '../hooks/useRenderer';
7+import CodeEditor from './CodeEditor';
8+
9+interface Props {
10+ mfiles: MFile[];
11+ onChangeFile: (name: string, source: string) => void;
12+ onResetFile: (name: string) => void;
13+ onResetAll: () => void;
14+ onClose: () => void;
15+ onRenderFull: () => void;
16+ error: RenderError | null;
17+ busy: boolean;
18+}
19+
20+function defaultSourceOf(name: string): string {
21+ return DEFAULT_MFILES.find((f) => f.name === name)?.source ?? '';
22+}
23+
24+export default function CodePanel({
25+ mfiles,
26+ onChangeFile,
27+ onResetFile,
28+ onResetAll,
29+ onClose,
30+ onRenderFull,
31+ error,
32+ busy,
33+}: Props) {
34+ const [selected, setSelected] = useState(
35+ mfiles.find((f) => EXECUTED_BY_APP.has(f.name))?.name ?? mfiles[0]?.name ?? '',
36+ );
37+ const file = mfiles.find((f) => f.name === selected) ?? mfiles[0];
38+ const dirty = (f: MFile) => f.source !== defaultSourceOf(f.name);
39+ const anyDirty = mfiles.some(dirty);
40+
41+ return (
42+ <div className="modal-backdrop" onClick={onClose}>
43+ <div className="code-panel" onClick={(e) => e.stopPropagation()}>
44+ <div className="code-header">
45+ <h2>MATLAB source</h2>
46+ <span className="code-note">
47+ Executed by numbl. Edits re-render as you type and are discarded on reload.
48+ </span>
49+ <button onClick={onRenderFull} disabled={busy}>
50+ Render full quality
51+ </button>
52+ <button onClick={onResetAll} disabled={!anyDirty}>
53+ Revert all
54+ </button>
55+ <button onClick={onClose}>Close</button>
56+ </div>
57+ <div className="code-body">
58+ <div className="code-files">
59+ {mfiles.map((f) => (
60+ <button
61+ key={f.name}
62+ className={f.name === file?.name ? 'active' : ''}
63+ onClick={() => setSelected(f.name)}
64+ >
65+ <span className="fname">
66+ {f.name}
67+ {dirty(f) ? <span className="dot" title="edited" /> : null}
68+ </span>
69+ <span className="ftag">{EXECUTED_BY_APP.has(f.name) ? 'in app' : 'CLI'}</span>
70+ </button>
71+ ))}
72+ </div>
73+ <div className="code-main">
74+ {file && (
75+ <>
76+ <div className="code-toolbar">
77+ <span>{file.name}</span>
78+ {!EXECUTED_BY_APP.has(file.name) && (
79+ <span className="code-hint">
80+ entry point for the numbl CLI; the app does not run this file
81+ </span>
82+ )}
83+ <button
84+ onClick={() => onResetFile(file.name)}
85+ disabled={!dirty(file)}
86+ style={{ marginLeft: 'auto' }}
87+ >
88+ Revert file
89+ </button>
90+ </div>
91+ <CodeEditor
92+ key={file.name}
93+ source={file.source}
94+ onChange={(src) => onChangeFile(file.name, src)}
95+ />
96+ </>
97+ )}
98+ </div>
99+ </div>
100+ {error && (
101+ <div className="code-error">
102+ <strong>{error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}</strong>{' '}
103+ {error.message}
104+ </div>
105+ )}
106+ </div>
107+ </div>
108+ );
109+}
webapp/src/components/DiagnosticsPanel.tsxadded+58−0View file
@@ -0,0 +1,58 @@
1+import { NOTE_NAMES } from '../params';
2+import type { RenderResponse } from '../worker/renderWorker';
3+
4+export default function DiagnosticsPanel({ result }: { result: RenderResponse }) {
5+ const d = result.diag;
6+ const modes = Array.from(result.modeF);
7+ const shown = modes.slice(0, 8);
8+ return (
9+ <div className="diag">
10+ <table>
11+ <tbody>
12+ {Array.from(result.noteIdx).map((idx, i) => (
13+ <tr key={i}>
14+ <th>{d.noteCount === 1 ? 'Note' : `Note ${i + 1}`}</th>
15+ <td>
16+ {NOTE_NAMES[idx - 1] ?? '?'} ·{' '}
17+ {result.fPlay[i] > 0
18+ ? `played ${result.fPlay[i].toFixed(1)} Hz (${result.cents[i] >= 0 ? '+' : ''}${result.cents[i].toFixed(0)} cents of ${result.fTarget[i].toFixed(1)})`
19+ : 'did not speak'}
20+ </td>
21+ </tr>
22+ ))}
23+ <tr>
24+ <th>Bore modes</th>
25+ <td className="partials">
26+ {shown.map((f) => f.toFixed(0)).join(', ')}
27+ {modes.length > shown.length ? `, … (${modes.length} total)` : ''} Hz
28+ </td>
29+ </tr>
30+ <tr>
31+ <th>Listener</th>
32+ <td>{(d.listenerDistM * 100).toFixed(0)} cm from the window</td>
33+ </tr>
34+ <tr>
35+ <th>Peak pressure</th>
36+ <td>
37+ {d.peakPa.toFixed(3)} Pa ({d.peakSplDb.toFixed(1)} dB SPL)
38+ </td>
39+ </tr>
40+ <tr>
41+ <th>Length</th>
42+ <td>
43+ {d.durationS.toFixed(2)} s · {result.quality} · oversample {d.oversample}
44+ </td>
45+ </tr>
46+ </tbody>
47+ </table>
48+ <div className="cachenote">
49+ bore solve: {d.basisCached ? 'cached' : `${(d.basisMs / 1000).toFixed(2)} s`} · jet:{' '}
50+ {d.jetCached ? 'cached' : `${(d.jetMs / 1000).toFixed(2)} s`} · radiation:{' '}
51+ {(d.radiateMs / 1000).toFixed(2)} s
52+ </div>
53+ {d.warnings.map((w, i) => (
54+ <div key={i} className="warning">⚠ {w}</div>
55+ ))}
56+ </div>
57+ );
58+}
webapp/src/components/ParamPanel.tsxadded+154−0View file
@@ -0,0 +1,154 @@
1+import { NOTE_NAMES, PARAM_GROUPS, type ParamSpec, type RecorderParams } from '../params';
2+
3+interface Props {
4+ params: RecorderParams;
5+ onChange: (key: keyof RecorderParams, value: number | boolean) => void;
6+}
7+
8+const SLIDER_STEPS = 1000;
9+
10+function sliderFromValue(spec: ParamSpec, v: number): number {
11+ if (spec.log) {
12+ const t = (Math.log(v) - Math.log(spec.min)) / (Math.log(spec.max) - Math.log(spec.min));
13+ return Math.round(t * SLIDER_STEPS);
14+ }
15+ return Math.round(((v - spec.min) / (spec.max - spec.min)) * SLIDER_STEPS);
16+}
17+
18+function valueFromSlider(spec: ParamSpec, s: number): number {
19+ const t = s / SLIDER_STEPS;
20+ let v: number;
21+ if (spec.log) {
22+ v = Math.exp(Math.log(spec.min) + t * (Math.log(spec.max) - Math.log(spec.min)));
23+ // keep three significant digits so the displayed number stays readable
24+ const mag = 10 ** Math.floor(Math.log10(Math.abs(v)));
25+ v = Math.round((v / mag) * 100) / 100 * mag;
26+ } else {
27+ v = spec.min + t * (spec.max - spec.min);
28+ if (spec.step) v = Math.round(v / spec.step) * spec.step;
29+ }
30+ if (spec.integer) v = Math.round(v);
31+ return Math.min(spec.max, Math.max(spec.min, v));
32+}
33+
34+function formatValue(v: number): string {
35+ if (v === 0) return '0';
36+ const a = Math.abs(v);
37+ if (a >= 1e5 || a < 1e-3) return v.toExponential(2);
38+ return String(Math.round(v * 1e6) / 1e6);
39+}
40+
41+function ParamRow({ spec, params, onChange }: Props & { spec: ParamSpec }) {
42+ const v = params[spec.key] as number;
43+ if (spec.note) {
44+ return (
45+ <div className="param-row" title={spec.help ?? ''}>
46+ <div className="param-label">
47+ <span className="name">{spec.label}</span>
48+ </div>
49+ <div className="param-controls">
50+ <select
51+ value={Math.min(NOTE_NAMES.length, Math.max(1, Math.round(v)))}
52+ onChange={(e) => onChange(spec.key, Number(e.target.value))}
53+ >
54+ {NOTE_NAMES.map((name, i) => (
55+ <option key={name} value={i + 1}>
56+ {name}
57+ {i >= 10 ? ' (pinched)' : ''}
58+ </option>
59+ ))}
60+ </select>
61+ </div>
62+ </div>
63+ );
64+ }
65+ return (
66+ <div className="param-row" title={spec.help ?? ''}>
67+ <div className="param-label">
68+ <span className="name">
69+ {spec.label}
70+ {spec.unit ? <span style={{ color: 'var(--muted)' }}> ({spec.unit})</span> : null}
71+ </span>
72+ </div>
73+ <div className="param-controls">
74+ <input
75+ type="range"
76+ min={0}
77+ max={SLIDER_STEPS}
78+ value={sliderFromValue(spec, v)}
79+ onChange={(e) => onChange(spec.key, valueFromSlider(spec, Number(e.target.value)))}
80+ />
81+ <input
82+ type="number"
83+ value={formatValue(v)}
84+ step={spec.step ?? 'any'}
85+ onChange={(e) => {
86+ const nv = Number(e.target.value);
87+ if (Number.isFinite(nv)) onChange(spec.key, spec.integer ? Math.round(nv) : nv);
88+ }}
89+ />
90+ </div>
91+ </div>
92+ );
93+}
94+
95+export default function ParamPanel({ params, onChange }: Props) {
96+ return (
97+ <div>
98+ {PARAM_GROUPS.map((g, gi) => (
99+ <details key={g.name} className="group" open={gi < 2}>
100+ <summary>{g.name}</summary>
101+ <div className="group-body">
102+ {g.specs
103+ .filter(
104+ (spec) =>
105+ !(
106+ (spec.key === 'note2' && params.note_count < 2) ||
107+ (spec.key === 'note3' && params.note_count < 3) ||
108+ (spec.key === 'gap_ms' && params.note_count < 2)
109+ ),
110+ )
111+ .map((spec) => (
112+ <ParamRow key={spec.key} spec={spec} params={params} onChange={onChange} />
113+ ))}
114+ {g.name === 'Output' && (
115+ <>
116+ <div className="param-row">
117+ <div className="param-label">
118+ <span className="name">Normalize peak</span>
119+ </div>
120+ <div className="param-controls">
121+ <input
122+ type="checkbox"
123+ checked={params.normalize}
124+ onChange={(e) => onChange('normalize', e.target.checked)}
125+ />
126+ <span className="param-help">
127+ {params.normalize
128+ ? `scaled to ${params.peak_dbfs} dBFS`
129+ : `absolute: ${params.full_scale_pa} Pa = 0 dBFS`}
130+ </span>
131+ </div>
132+ </div>
133+ <div className="param-row">
134+ <div className="param-label">
135+ <span className="name">WAV bit depth</span>
136+ </div>
137+ <div className="param-controls">
138+ <select
139+ value={params.bit_depth}
140+ onChange={(e) => onChange('bit_depth', Number(e.target.value))}
141+ >
142+ <option value={16}>16-bit</option>
143+ <option value={24}>24-bit</option>
144+ </select>
145+ </div>
146+ </div>
147+ </>
148+ )}
149+ </div>
150+ </details>
151+ ))}
152+ </div>
153+ );
154+}
webapp/src/components/Scene3D.tsxadded+254−0View file
@@ -0,0 +1,254 @@
1+// 3D view of the geometry: the recorder hanging vertically (window at the
2+// origin, foot down along -y, tone holes opening toward +z, thumb hole
3+// toward -z), the fingering of the first note shown as finger pads on the
4+// closed holes, and the listener. three.js is y-up, so the mapping is the
5+// identity.
6+
7+import { useEffect, useRef } from 'react';
8+import * as THREE from 'three';
9+import { OrbitControls } from 'three/examples/jsm/controls/OrbitControls.js';
10+import type { RecorderParams } from '../params';
11+import { FINGERINGS, HOLE_FRAC, HOLE_RAD } from '../lib/instrument';
12+
13+function makeLabel(text: string, color: string, scale = 1): THREE.Sprite {
14+ const canvas = document.createElement('canvas');
15+ canvas.width = 320;
16+ canvas.height = 64;
17+ const ctx = canvas.getContext('2d')!;
18+ ctx.font = '600 36px system-ui, sans-serif';
19+ ctx.fillStyle = color;
20+ ctx.textAlign = 'center';
21+ ctx.textBaseline = 'middle';
22+ ctx.fillText(text, 160, 32);
23+ const tex = new THREE.CanvasTexture(canvas);
24+ const mat = new THREE.SpriteMaterial({ map: tex, depthTest: false, transparent: true });
25+ const sprite = new THREE.Sprite(mat);
26+ sprite.scale.set(0.11 * scale, 0.022 * scale, 1);
27+ return sprite;
28+}
29+
30+interface Props {
31+ params: RecorderParams;
32+}
33+
34+export default function Scene3D({ params }: Props) {
35+ const mountRef = useRef<HTMLDivElement | null>(null);
36+ const sceneRef = useRef<{
37+ scene: THREE.Scene;
38+ group: THREE.Group | null;
39+ renderer: THREE.WebGLRenderer;
40+ camera: THREE.PerspectiveCamera;
41+ } | null>(null);
42+
43+ // one-time scene setup
44+ useEffect(() => {
45+ const mount = mountRef.current!;
46+ const renderer = new THREE.WebGLRenderer({ antialias: true });
47+ renderer.setPixelRatio(window.devicePixelRatio);
48+ renderer.setClearColor(0x14161a);
49+ mount.appendChild(renderer.domElement);
50+
51+ const scene = new THREE.Scene();
52+ const camera = new THREE.PerspectiveCamera(45, 1, 0.005, 20);
53+ camera.position.set(0.30, 0.06, 0.42);
54+
55+ const controls = new OrbitControls(camera, renderer.domElement);
56+ controls.target.set(0, -0.14, 0);
57+ controls.enableDamping = true;
58+
59+ scene.add(new THREE.AmbientLight(0xffffff, 0.55));
60+ const dir = new THREE.DirectionalLight(0xffffff, 1.1);
61+ dir.position.set(0.5, 1.0, 0.6);
62+ scene.add(dir);
63+ const dir2 = new THREE.DirectionalLight(0x8899ff, 0.3);
64+ dir2.position.set(-0.6, 0.3, -0.5);
65+ scene.add(dir2);
66+
67+ const grid = new THREE.GridHelper(1.0, 32, 0x3a4048, 0x262b33);
68+ grid.position.y = -0.42;
69+ scene.add(grid);
70+
71+ // model axes: x red, y green (along the bore, up), z blue (front)
72+ const axisLen = 0.06;
73+ const axes: [THREE.Vector3, number, string][] = [
74+ [new THREE.Vector3(axisLen, 0, 0), 0xe57373, '+x'],
75+ [new THREE.Vector3(0, axisLen, 0), 0x81c784, '+y'],
76+ [new THREE.Vector3(0, 0, axisLen), 0x64b5f6, '+z'],
77+ ];
78+ for (const [v, color, label] of axes) {
79+ const geo = new THREE.BufferGeometry().setFromPoints([new THREE.Vector3(0, 0, 0), v]);
80+ scene.add(new THREE.Line(geo, new THREE.LineBasicMaterial({ color })));
81+ const s = makeLabel(label, '#' + color.toString(16).padStart(6, '0'));
82+ s.position.copy(v).multiplyScalar(1.15);
83+ scene.add(s);
84+ }
85+
86+ sceneRef.current = { scene, group: null, renderer, camera };
87+
88+ let raf = 0;
89+ const animate = () => {
90+ raf = requestAnimationFrame(animate);
91+ controls.update();
92+ renderer.render(scene, camera);
93+ };
94+ animate();
95+
96+ const resize = () => {
97+ const w = mount.clientWidth;
98+ const h = mount.clientHeight;
99+ renderer.setSize(w, h, false);
100+ camera.aspect = w / h;
101+ camera.updateProjectionMatrix();
102+ };
103+ resize();
104+ const ro = new ResizeObserver(resize);
105+ ro.observe(mount);
106+
107+ return () => {
108+ cancelAnimationFrame(raf);
109+ ro.disconnect();
110+ controls.dispose();
111+ renderer.dispose();
112+ mount.removeChild(renderer.domElement);
113+ sceneRef.current = null;
114+ };
115+ }, []);
116+
117+ // rebuild the geometry group whenever relevant parameters change
118+ useEffect(() => {
119+ const ctx = sceneRef.current;
120+ if (!ctx) return;
121+ if (ctx.group) {
122+ ctx.scene.remove(ctx.group);
123+ ctx.group.traverse((o) => {
124+ const mesh = o as THREE.Mesh;
125+ if (mesh.geometry) mesh.geometry.dispose();
126+ });
127+ }
128+
129+ const g = new THREE.Group();
130+ const L = params.bore_length_m;
131+ const rHead = params.bore_r_head_m + params.wall_m;
132+ const rFoot = params.bore_r_foot_m + params.wall_m;
133+ const wood = new THREE.MeshStandardMaterial({
134+ color: 0xa8763e,
135+ roughness: 0.55,
136+ metalness: 0.05,
137+ });
138+ const darkWood = new THREE.MeshStandardMaterial({
139+ color: 0x6e4a26,
140+ roughness: 0.6,
141+ metalness: 0.05,
142+ });
143+ const holeMat = new THREE.MeshStandardMaterial({ color: 0x17130e, roughness: 0.9 });
144+
145+ // body: window (y = 0) down to the foot
146+ const body = new THREE.Mesh(new THREE.CylinderGeometry(rHead, rFoot, L, 28), wood);
147+ body.position.set(0, -L / 2, 0);
148+ g.add(body);
149+
150+ // head above the window: beak with the windway
151+ const headLen = 0.045;
152+ const head = new THREE.Mesh(new THREE.CylinderGeometry(rHead * 1.06, rHead, headLen, 28), darkWood);
153+ head.position.set(0, headLen / 2 + params.cutup_m, 0);
154+ g.add(head);
155+ const beak = new THREE.Mesh(
156+ new THREE.SphereGeometry(rHead * 1.06, 20, 12, 0, Math.PI * 2, 0, Math.PI / 2),
157+ darkWood,
158+ );
159+ beak.position.set(0, headLen + params.cutup_m, 0);
160+ g.add(beak);
161+
162+ // window: dark opening on the +z side between flue exit and labium
163+ const win = new THREE.Mesh(
164+ new THREE.BoxGeometry(params.window_width_m, params.cutup_m, 0.0015),
165+ holeMat,
166+ );
167+ win.position.set(0, params.cutup_m / 2, rHead - 0.0002);
168+ g.add(win);
169+
170+ // tone holes; a fingering pad covers each closed hole
171+ const fing = FINGERINGS[Math.min(14, Math.max(1, Math.round(params.note1))) - 1];
172+ const pad = new THREE.MeshStandardMaterial({ color: 0xd8a37a, roughness: 0.8 });
173+ for (let j = 0; j < 8; j++) {
174+ const y = -HOLE_FRAC[j] * L;
175+ const r = HOLE_RAD[j] * params.hole_scale;
176+ const frac = HOLE_FRAC[j];
177+ const rAt = rHead + (rFoot - rHead) * frac;
178+ const zdir = j === 0 ? -1 : 1; // thumb hole on the back
179+ const hole = new THREE.Mesh(new THREE.CylinderGeometry(r, r, 0.004, 16), holeMat);
180+ hole.rotation.x = Math.PI / 2;
181+ hole.position.set(0, y, zdir * rAt);
182+ g.add(hole);
183+ const state = fing[j];
184+ if (state < 0.75) {
185+ // closed or pinched: a finger pad (smaller for the pinch)
186+ const pr = state < 0.25 ? r * 1.9 : r * 1.1;
187+ const cap = new THREE.Mesh(new THREE.SphereGeometry(pr, 14, 10), pad);
188+ cap.scale.z = 0.45;
189+ cap.position.set(0, y, zdir * (rAt + 0.001));
190+ g.add(cap);
191+ }
192+ }
193+
194+ // foot opening
195+ const foot = new THREE.Mesh(
196+ new THREE.CylinderGeometry(params.bore_r_foot_m, params.bore_r_foot_m, 0.002, 20),
197+ holeMat,
198+ );
199+ foot.position.set(0, -L - 0.0005, 0);
200+ g.add(foot);
201+
202+ // listener: microphone marker + line from the window
203+ const lis = new THREE.Vector3(params.listener_x_m, params.listener_y_m, params.listener_z_m);
204+ const micHead = new THREE.Mesh(
205+ new THREE.SphereGeometry(0.016, 20, 16),
206+ new THREE.MeshStandardMaterial({ color: 0x4dd0e1, emissive: 0x0c3b42, roughness: 0.35 }),
207+ );
208+ micHead.position.copy(lis);
209+ g.add(micHead);
210+ const micBody = new THREE.Mesh(
211+ new THREE.CylinderGeometry(0.008, 0.008, 0.05, 12),
212+ new THREE.MeshStandardMaterial({ color: 0x455a64, roughness: 0.6 }),
213+ );
214+ micBody.position.copy(lis).add(new THREE.Vector3(0, -0.033, 0));
215+ g.add(micBody);
216+ const windowPos = new THREE.Vector3(0, 0, 0);
217+ const lineGeo = new THREE.BufferGeometry().setFromPoints([windowPos, lis]);
218+ const line = new THREE.Line(
219+ lineGeo,
220+ new THREE.LineDashedMaterial({ color: 0x4dd0e1, dashSize: 0.012, gapSize: 0.009 }),
221+ );
222+ line.computeLineDistances();
223+ g.add(line);
224+ const dist = windowPos.distanceTo(lis);
225+ const micLabel = makeLabel(`mic ${(dist * 100).toFixed(0)} cm`, '#4dd0e1');
226+ micLabel.position.copy(lis).add(new THREE.Vector3(0, 0.045, 0));
227+ g.add(micLabel);
228+
229+ const winLabel = makeLabel('window', '#e8a33d', 0.8);
230+ winLabel.position.set(0.05, 0.012, rHead + 0.01);
231+ g.add(winLabel);
232+
233+ ctx.scene.add(g);
234+ ctx.group = g;
235+ }, [
236+ params.bore_length_m,
237+ params.bore_r_head_m,
238+ params.bore_r_foot_m,
239+ params.wall_m,
240+ params.hole_scale,
241+ params.cutup_m,
242+ params.window_width_m,
243+ params.note1,
244+ params.listener_x_m,
245+ params.listener_y_m,
246+ params.listener_z_m,
247+ ]);
248+
249+ return (
250+ <div className="scene3d" ref={mountRef}>
251+ <div className="scene-hint">drag to orbit · scroll to zoom · right-drag to pan</div>
252+ </div>
253+ );
254+}
webapp/src/components/plots/AdmittancePlot.tsxadded+111−0View file
@@ -0,0 +1,111 @@
1+// The passive resonator of the first note's fingering: |Y(f)| of the series
2+// loop (window mass + bore input impedance) on a log-log axis, with the
3+// fitted modal centre frequencies overlaid and the frequency the jet
4+// actually locked to marked. The playable notes are the peaks; where the
5+// played line sits relative to the first peak is the intonation.
6+
7+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
8+
9+interface Props {
10+ gridF: Float32Array;
11+ gridY: Float32Array;
12+ modeF: Float32Array;
13+ fPlay: number;
14+}
15+
16+const DB_RANGE = 60;
17+
18+export default function AdmittancePlot({ gridF, gridY, modeF, fPlay }: Props) {
19+ const ref = useCanvas(
20+ (ctx, w, h) => {
21+ drawFrame(ctx, w, h);
22+ if (gridF.length < 2) return;
23+ const padL = 34;
24+ const padB = 18;
25+ const pw = w - padL - 6;
26+ const ph = h - padB - 8;
27+ const fLo = Math.max(80, gridF[0]);
28+ const fHi = gridF[gridF.length - 1];
29+ const fx = (f: number) =>
30+ padL + ((Math.log10(f) - Math.log10(fLo)) / (Math.log10(fHi) - Math.log10(fLo))) * pw;
31+
32+ ctx.font = '10px system-ui';
33+ ctx.textAlign = 'center';
34+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000]) {
35+ if (f < fLo || f > fHi) continue;
36+ const x = fx(f);
37+ ctx.strokeStyle = PLOT_COLORS.grid;
38+ ctx.beginPath();
39+ ctx.moveTo(x, 8);
40+ ctx.lineTo(x, 8 + ph);
41+ ctx.stroke();
42+ ctx.fillStyle = PLOT_COLORS.axis;
43+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 5);
44+ }
45+
46+ let yPk = 0;
47+ for (let i = 0; i < gridY.length; i++) yPk = Math.max(yPk, gridY[i]);
48+ if (yPk <= 0) return;
49+
50+ ctx.textAlign = 'right';
51+ for (let db = 0; db >= -DB_RANGE; db -= 20) {
52+ const yy = 8 + (-db / DB_RANGE) * ph;
53+ ctx.strokeStyle = PLOT_COLORS.grid;
54+ ctx.beginPath();
55+ ctx.moveTo(padL, yy);
56+ ctx.lineTo(padL + pw, yy);
57+ ctx.stroke();
58+ ctx.fillStyle = PLOT_COLORS.axis;
59+ ctx.fillText(`${db}`, padL - 4, yy + 3);
60+ }
61+
62+ // fitted modes, then the played frequency, then the trace on top
63+ ctx.setLineDash([4, 3]);
64+ for (let i = 0; i < modeF.length; i++) {
65+ const f = modeF[i];
66+ if (f < fLo || f > fHi) continue;
67+ ctx.strokeStyle = 'rgba(232, 163, 61, 0.45)';
68+ ctx.beginPath();
69+ ctx.moveTo(fx(f), 8);
70+ ctx.lineTo(fx(f), 8 + ph);
71+ ctx.stroke();
72+ }
73+ ctx.setLineDash([]);
74+ if (fPlay > fLo && fPlay < fHi) {
75+ ctx.strokeStyle = 'rgba(111, 191, 115, 0.9)';
76+ ctx.lineWidth = 1.5;
77+ ctx.beginPath();
78+ ctx.moveTo(fx(fPlay), 8);
79+ ctx.lineTo(fx(fPlay), 8 + ph);
80+ ctx.stroke();
81+ ctx.lineWidth = 1;
82+ }
83+
84+ ctx.strokeStyle = PLOT_COLORS.trace;
85+ ctx.beginPath();
86+ let started = false;
87+ for (let i = 0; i < gridF.length; i++) {
88+ if (gridF[i] < fLo) continue;
89+ const x = fx(gridF[i]);
90+ const db = 20 * Math.log10(gridY[i] / yPk + 1e-12);
91+ const yy = 8 + Math.min(ph, (-db / DB_RANGE) * ph);
92+ if (!started) {
93+ ctx.moveTo(x, yy);
94+ started = true;
95+ } else {
96+ ctx.lineTo(x, yy);
97+ }
98+ }
99+ ctx.stroke();
100+
101+ ctx.textAlign = 'left';
102+ ctx.fillStyle = 'rgba(232, 163, 61, 0.9)';
103+ ctx.fillText('fitted modes', padL + 4, 18);
104+ ctx.fillStyle = 'rgba(111, 191, 115, 0.9)';
105+ ctx.fillText('played', padL + 74, 18);
106+ },
107+ [gridF, gridY, modeF, fPlay],
108+ 200,
109+ );
110+ return <canvas ref={ref} />;
111+}
webapp/src/components/plots/BoreProfile.tsxadded+74−0View file
@@ -0,0 +1,74 @@
1+// Small multiples of the pressure standing wave |p(x)| along the bore for
2+// the lowest resonances of the first note's fingering. The window is at the
3+// left, the foot at the right; an open tone hole shows up as the profile
4+// being pinned near zero beyond it, which is how a fingering shortens the
5+// instrument without shortening the tube.
6+
7+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
8+
9+interface Props {
10+ profX: Float32Array;
11+ profP: Float32Array; // nsta x n, column-major
12+ modeF: Float32Array;
13+}
14+
15+const SHOWN = 8;
16+const CELL_H = 100;
17+
18+export default function BoreProfile({ profX, profP, modeF }: Props) {
19+ const nShown = Math.min(SHOWN, modeF.length);
20+ const rows = Math.ceil(nShown / 4);
21+ const ref = useCanvas(
22+ (ctx, w, h) => {
23+ drawFrame(ctx, w, h);
24+ const nsta = profX.length;
25+ if (nsta === 0 || nShown === 0) return;
26+ const cols = 4;
27+ const cw = w / cols;
28+ const ch = h / rows;
29+ const xMax = profX[nsta - 1];
30+
31+ for (let m = 0; m < nShown; m++) {
32+ const x0 = (m % cols) * cw + 10;
33+ const y0 = Math.floor(m / cols) * ch + 12;
34+ const iw = cw - 20;
35+ const ih = ch - 34;
36+ const off = m * nsta;
37+
38+ let pMax = 1e-12;
39+ for (let i = 0; i < nsta; i++) pMax = Math.max(pMax, profP[off + i]);
40+
41+ // bore axis
42+ ctx.strokeStyle = 'rgba(139, 147, 161, 0.3)';
43+ ctx.beginPath();
44+ ctx.moveTo(x0, y0 + ih);
45+ ctx.lineTo(x0 + iw, y0 + ih);
46+ ctx.stroke();
47+
48+ ctx.strokeStyle = PLOT_COLORS.trace;
49+ ctx.lineWidth = 1.6;
50+ ctx.beginPath();
51+ for (let i = 0; i < nsta; i++) {
52+ const x = x0 + (profX[i] / xMax) * iw;
53+ const y = y0 + ih - (profP[off + i] / pMax) * ih;
54+ if (i === 0) ctx.moveTo(x, y);
55+ else ctx.lineTo(x, y);
56+ }
57+ ctx.stroke();
58+ ctx.lineWidth = 1;
59+
60+ ctx.font = '10px system-ui';
61+ ctx.textAlign = 'center';
62+ ctx.fillStyle = PLOT_COLORS.accent;
63+ ctx.fillText(`${modeF[m].toFixed(0)} Hz`, x0 + iw / 2, y0 + ih + 14);
64+ }
65+ ctx.textAlign = 'left';
66+ ctx.fillStyle = PLOT_COLORS.axis;
67+ ctx.font = '10px system-ui';
68+ ctx.fillText('window → foot', 10, 9);
69+ },
70+ [profX, profP, modeF, nShown, rows],
71+ rows * CELL_H,
72+ );
73+ return <canvas ref={ref} />;
74+}
webapp/src/components/plots/BreathPlot.tsxadded+92−0View file
@@ -0,0 +1,92 @@
1+// The jet inflow envelope over the whole gesture: the mean flow rises with
2+// the breath envelope, the oscillation rides on it once the note locks, and
3+// the attack chiff is visible as the ragged region before lock. Note onsets
4+// are marked with the frequency each note locked to.
5+
6+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
7+
8+interface Props {
9+ qPlot: Float32Array;
10+ qPlotRate: number;
11+ noteOn: Float32Array;
12+ fPlay: Float32Array;
13+ durationS: number;
14+}
15+
16+export default function BreathPlot({ qPlot, qPlotRate, noteOn, fPlay, durationS }: Props) {
17+ const ref = useCanvas(
18+ (ctx, w, h) => {
19+ drawFrame(ctx, w, h);
20+ const padL = 40;
21+ const padB = 16;
22+ const pw = w - padL - 8;
23+ const ph = h - padB - 12;
24+ const top = 10;
25+
26+ let qMax = 0;
27+ for (let i = 0; i < qPlot.length; i++) qMax = Math.max(qMax, qPlot[i]);
28+ if (qMax <= 0) qMax = 1;
29+
30+ const tx = (t: number) => padL + (t / durationS) * pw;
31+
32+ ctx.font = '10px system-ui';
33+ ctx.textAlign = 'center';
34+ const tStep = durationS > 4 ? 1.0 : durationS > 1.5 ? 0.5 : 0.1;
35+ for (let t = 0; t <= durationS + 1e-9; t += tStep) {
36+ const x = tx(t);
37+ ctx.strokeStyle = PLOT_COLORS.grid;
38+ ctx.beginPath();
39+ ctx.moveTo(x, top);
40+ ctx.lineTo(x, top + ph);
41+ ctx.stroke();
42+ ctx.fillStyle = PLOT_COLORS.axis;
43+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
44+ }
45+ ctx.textAlign = 'right';
46+ ctx.fillStyle = PLOT_COLORS.axis;
47+ ctx.fillText(`${(qMax * 1000).toFixed(2)} L/s`, padL - 4, top + 4);
48+ ctx.fillText('0', padL - 4, top + ph);
49+
50+ // per-pixel min/max envelope of the decimated jet flow
51+ ctx.strokeStyle = PLOT_COLORS.trace;
52+ ctx.lineWidth = 1;
53+ for (let px = 0; px < pw; px++) {
54+ const i0 = Math.floor((px / pw) * durationS * qPlotRate);
55+ const i1 = Math.min(qPlot.length, Math.floor(((px + 1) / pw) * durationS * qPlotRate) + 1);
56+ if (i0 >= qPlot.length) break;
57+ let lo = Infinity;
58+ let hi = -Infinity;
59+ for (let i = i0; i < i1; i++) {
60+ lo = Math.min(lo, qPlot[i]);
61+ hi = Math.max(hi, qPlot[i]);
62+ }
63+ if (lo > hi) continue;
64+ const y0 = top + ph - (hi / qMax) * ph * 0.92;
65+ const y1 = top + ph - (lo / qMax) * ph * 0.92;
66+ ctx.beginPath();
67+ ctx.moveTo(padL + px + 0.5, y0);
68+ ctx.lineTo(padL + px + 0.5, Math.max(y0 + 0.5, y1));
69+ ctx.stroke();
70+ }
71+
72+ // note onsets with locked frequency
73+ ctx.textAlign = 'left';
74+ for (let i = 0; i < noteOn.length; i++) {
75+ const x = tx(noteOn[i]);
76+ ctx.strokeStyle = 'rgba(232, 163, 61, 0.6)';
77+ ctx.setLineDash([4, 3]);
78+ ctx.beginPath();
79+ ctx.moveTo(x, top);
80+ ctx.lineTo(x, top + ph);
81+ ctx.stroke();
82+ ctx.setLineDash([]);
83+ ctx.fillStyle = PLOT_COLORS.accent;
84+ const f = fPlay[i];
85+ ctx.fillText(f > 0 ? `${f.toFixed(0)} Hz` : 'did not speak', x + 3, top + 10);
86+ }
87+ },
88+ [qPlot, qPlotRate, noteOn, fPlay, durationS],
89+ 150,
90+ );
91+ return <canvas ref={ref} />;
92+}
webapp/src/components/plots/SpectrogramPlot.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+// Time-frequency view of the rendered note: the attack chiff, the harmonic
2+// comb of the locked oscillation, and the breath-noise floor between them.
3+
4+import { stft } from '../../lib/fft';
5+import { useCanvas, drawFrame, PLOT_COLORS, colormap } from './canvasUtils';
6+
7+interface Props {
8+ y: Float32Array;
9+ sampleRate: number;
10+}
11+
12+export default function SpectrogramPlot({ y, sampleRate }: Props) {
13+ const ref = useCanvas(
14+ (ctx, w, h) => {
15+ drawFrame(ctx, w, h);
16+ const padL = 34;
17+ const padB = 16;
18+ const pw = w - padL - 6;
19+ const ph = h - padB - 6;
20+
21+ const { frames, bins, db, floorDb } = stft(y, 512, 128);
22+
23+ const img = ctx.createImageData(frames, bins);
24+ for (let f = 0; f < frames; f++) {
25+ for (let b = 0; b < bins; b++) {
26+ const v = db[f * bins + b];
27+ const t = Math.min(1, Math.max(0, (v - floorDb) / -floorDb));
28+ const [r, g, bl] = colormap(t);
29+ // flip vertically: low frequency at the bottom
30+ const row = bins - 1 - b;
31+ const idx = (row * frames + f) * 4;
32+ img.data[idx] = r;
33+ img.data[idx + 1] = g;
34+ img.data[idx + 2] = bl;
35+ img.data[idx + 3] = 255;
36+ }
37+ }
38+ const off = new OffscreenCanvas(frames, bins);
39+ const octx = off.getContext('2d')!;
40+ octx.putImageData(img, 0, 0);
41+ ctx.imageSmoothingEnabled = true;
42+ ctx.drawImage(off, padL, 6, pw, ph);
43+
44+ // axes
45+ ctx.fillStyle = PLOT_COLORS.axis;
46+ ctx.font = '10px system-ui';
47+ ctx.textAlign = 'right';
48+ const nyq = sampleRate / 2;
49+ for (let kf = 0; kf <= nyq; kf += 4000) {
50+ const yy = 6 + ph - (kf / nyq) * ph;
51+ ctx.fillText(kf >= 1000 ? `${kf / 1000}k` : `${kf}`, padL - 4, yy + 3);
52+ }
53+ ctx.textAlign = 'center';
54+ const durS = y.length / sampleRate;
55+ for (let t = 0; t <= durS + 1e-9; t += 0.1) {
56+ const x = padL + (t / durS) * pw;
57+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
58+ }
59+ },
60+ [y, sampleRate],
61+ 200,
62+ );
63+ return <canvas ref={ref} />;
64+}
webapp/src/components/plots/SpectrumPlot.tsxadded+120−0View file
@@ -0,0 +1,120 @@
1+// Magnitude spectrum of the rendered note, log-frequency axis, with the
2+// fitted bore resonances and the harmonics of the played fundamental
3+// overlaid. A recorder's spectrum is a harmonic comb over a breath-noise
4+// floor: the harmonics come from the jet nonlinearity, and only the first
5+// one or two sit on bore resonances.
6+
7+import { magnitudeSpectrum } from '../../lib/fft';
8+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
9+
10+interface Props {
11+ y: Float32Array;
12+ sampleRate: number;
13+ modeF: Float32Array;
14+ fPlay: Float32Array;
15+}
16+
17+const F_LO = 100;
18+const F_HI = 24000;
19+const DB_RANGE = 80;
20+
21+export default function SpectrumPlot({ y, sampleRate, modeF, fPlay }: Props) {
22+ const ref = useCanvas(
23+ (ctx, w, h) => {
24+ drawFrame(ctx, w, h);
25+ const padL = 34;
26+ const padB = 18;
27+ const pw = w - padL - 6;
28+ const ph = h - padB - 8;
29+ const fx = (f: number) =>
30+ padL + ((Math.log10(f) - Math.log10(F_LO)) / (Math.log10(F_HI) - Math.log10(F_LO))) * pw;
31+
32+ ctx.font = '10px system-ui';
33+ ctx.textAlign = 'center';
34+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000, 20000]) {
35+ const x = fx(f);
36+ ctx.strokeStyle = PLOT_COLORS.grid;
37+ ctx.beginPath();
38+ ctx.moveTo(x, 8);
39+ ctx.lineTo(x, 8 + ph);
40+ ctx.stroke();
41+ ctx.fillStyle = PLOT_COLORS.axis;
42+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 5);
43+ }
44+ ctx.textAlign = 'right';
45+ for (let db = 0; db >= -DB_RANGE; db -= 20) {
46+ const yy = 8 + (-db / DB_RANGE) * ph;
47+ ctx.strokeStyle = PLOT_COLORS.grid;
48+ ctx.beginPath();
49+ ctx.moveTo(padL, yy);
50+ ctx.lineTo(padL + pw, yy);
51+ ctx.stroke();
52+ ctx.fillStyle = PLOT_COLORS.axis;
53+ ctx.fillText(`${db}`, padL - 4, yy + 3);
54+ }
55+
56+ const { mag, size } = magnitudeSpectrum(y, 32768);
57+ let peak = 0;
58+ for (let i = 1; i < mag.length; i++) peak = Math.max(peak, mag[i]);
59+ if (peak <= 0) return;
60+
61+ // overlays first so the trace draws on top
62+ ctx.setLineDash([4, 3]);
63+ for (let i = 0; i < modeF.length; i++) {
64+ const f = modeF[i];
65+ if (f < F_LO || f > F_HI) continue;
66+ ctx.strokeStyle = 'rgba(232, 163, 61, 0.4)';
67+ ctx.beginPath();
68+ ctx.moveTo(fx(f), 8);
69+ ctx.lineTo(fx(f), 8 + ph);
70+ ctx.stroke();
71+ }
72+ if (fPlay.length > 0 && fPlay[0] > 0) {
73+ for (let hnum = 1; hnum * fPlay[0] < F_HI; hnum++) {
74+ const f = hnum * fPlay[0];
75+ if (f < F_LO) continue;
76+ ctx.strokeStyle = 'rgba(111, 191, 115, 0.35)';
77+ ctx.beginPath();
78+ ctx.moveTo(fx(f), 8);
79+ ctx.lineTo(fx(f), 8 + ph);
80+ ctx.stroke();
81+ }
82+ }
83+ ctx.setLineDash([]);
84+
85+ ctx.strokeStyle = PLOT_COLORS.trace;
86+ ctx.beginPath();
87+ let started = false;
88+ for (let px = 0; px <= pw; px++) {
89+ const f = 10 ** (Math.log10(F_LO) + (px / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
90+ const bin = Math.min(mag.length - 1, Math.round((f / sampleRate) * size));
91+ const fNext =
92+ 10 ** (Math.log10(F_LO) + ((px + 1) / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
93+ const binNext = Math.min(
94+ mag.length - 1,
95+ Math.max(bin + 1, Math.round((fNext / sampleRate) * size)),
96+ );
97+ let m = 0;
98+ for (let b = bin; b < binNext; b++) m = Math.max(m, mag[b]);
99+ const db = 20 * Math.log10(m / peak + 1e-12);
100+ const yy = 8 + Math.min(ph, (-db / DB_RANGE) * ph);
101+ if (!started) {
102+ ctx.moveTo(padL + px, yy);
103+ started = true;
104+ } else {
105+ ctx.lineTo(padL + px, yy);
106+ }
107+ }
108+ ctx.stroke();
109+
110+ ctx.textAlign = 'left';
111+ ctx.fillStyle = 'rgba(232, 163, 61, 0.9)';
112+ ctx.fillText('bore resonances', padL + 4, 18);
113+ ctx.fillStyle = 'rgba(111, 191, 115, 0.9)';
114+ ctx.fillText('harmonics of the played note', padL + 94, 18);
115+ },
116+ [y, sampleRate, modeF, fPlay],
117+ 200,
118+ );
119+ return <canvas ref={ref} />;
120+}
webapp/src/components/plots/WaveformPlot.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
2+
3+interface Props {
4+ y: Float32Array;
5+ sampleRate: number;
6+}
7+
8+export default function WaveformPlot({ y, sampleRate }: Props) {
9+ const ref = useCanvas(
10+ (ctx, w, h) => {
11+ drawFrame(ctx, w, h);
12+ const padL = 34;
13+ const padB = 16;
14+ const pw = w - padL - 6;
15+ const ph = h - padB - 6;
16+ const midY = 6 + ph / 2;
17+
18+ ctx.strokeStyle = PLOT_COLORS.grid;
19+ ctx.beginPath();
20+ ctx.moveTo(padL, midY);
21+ ctx.lineTo(padL + pw, midY);
22+ ctx.stroke();
23+
24+ // time gridlines every 0.1 s
25+ const durS = y.length / sampleRate;
26+ ctx.fillStyle = PLOT_COLORS.axis;
27+ ctx.font = '10px system-ui';
28+ ctx.textAlign = 'center';
29+ for (let t = 0; t <= durS + 1e-9; t += 0.1) {
30+ const x = padL + (t / durS) * pw;
31+ ctx.strokeStyle = PLOT_COLORS.grid;
32+ ctx.beginPath();
33+ ctx.moveTo(x, 6);
34+ ctx.lineTo(x, 6 + ph);
35+ ctx.stroke();
36+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
37+ }
38+ ctx.textAlign = 'right';
39+ ctx.fillText('+1', padL - 4, 12);
40+ ctx.fillText('-1', padL - 4, 6 + ph);
41+
42+ // min/max envelope per pixel column
43+ ctx.strokeStyle = PLOT_COLORS.trace;
44+ ctx.beginPath();
45+ const perPx = y.length / pw;
46+ for (let px = 0; px < pw; px++) {
47+ const i0 = Math.floor(px * perPx);
48+ const i1 = Math.min(y.length, Math.max(i0 + 1, Math.floor((px + 1) * perPx)));
49+ let mn = Infinity, mx = -Infinity;
50+ for (let i = i0; i < i1; i++) {
51+ if (y[i] < mn) mn = y[i];
52+ if (y[i] > mx) mx = y[i];
53+ }
54+ if (!Number.isFinite(mn)) continue;
55+ ctx.moveTo(padL + px + 0.5, midY - mx * (ph / 2));
56+ ctx.lineTo(padL + px + 0.5, midY - mn * (ph / 2) + 0.5);
57+ }
58+ ctx.stroke();
59+ },
60+ [y, sampleRate],
61+ 170,
62+ );
63+ return <canvas ref={ref} />;
64+}
webapp/src/components/plots/canvasUtils.tsadded+66−0View file
@@ -0,0 +1,66 @@
1+import { useEffect, useRef } from 'react';
2+
3+// Canvas hook: resizes to element width at devicePixelRatio and calls draw.
4+export function useCanvas(
5+ draw: (ctx: CanvasRenderingContext2D, w: number, h: number) => void,
6+ deps: unknown[],
7+ height: number,
8+) {
9+ const ref = useRef<HTMLCanvasElement | null>(null);
10+
11+ useEffect(() => {
12+ const canvas = ref.current;
13+ if (!canvas) return;
14+ const render = () => {
15+ const dpr = window.devicePixelRatio || 1;
16+ const w = canvas.clientWidth;
17+ canvas.width = Math.max(1, Math.round(w * dpr));
18+ canvas.height = Math.max(1, Math.round(height * dpr));
19+ canvas.style.height = `${height}px`;
20+ const ctx = canvas.getContext('2d')!;
21+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
22+ ctx.clearRect(0, 0, w, height);
23+ draw(ctx, w, height);
24+ };
25+ render();
26+ const ro = new ResizeObserver(render);
27+ ro.observe(canvas);
28+ return () => ro.disconnect();
29+ // eslint-disable-next-line react-hooks/exhaustive-deps
30+ }, deps);
31+
32+ return ref;
33+}
34+
35+export const PLOT_COLORS = {
36+ bg: '#181b20',
37+ grid: '#2a2f38',
38+ axis: '#8b93a1',
39+ trace: '#4dd0e1',
40+ accent: '#e8a33d',
41+ good: '#6fbf73',
42+ warn: '#e57373',
43+};
44+
45+export function drawFrame(ctx: CanvasRenderingContext2D, w: number, h: number) {
46+ ctx.fillStyle = PLOT_COLORS.bg;
47+ ctx.fillRect(0, 0, w, h);
48+}
49+
50+// Viridis-like colormap for the spectrogram, t in [0,1].
51+export function colormap(t: number): [number, number, number] {
52+ const stops: [number, number, number][] = [
53+ [13, 8, 60], [84, 39, 143], [140, 41, 129], [196, 66, 90],
54+ [237, 105, 37], [251, 155, 6], [247, 209, 61], [252, 255, 164],
55+ ];
56+ const x = Math.min(0.9999, Math.max(0, t)) * (stops.length - 1);
57+ const i = Math.floor(x);
58+ const f = x - i;
59+ const a = stops[i];
60+ const b = stops[i + 1];
61+ return [
62+ Math.round(a[0] + f * (b[0] - a[0])),
63+ Math.round(a[1] + f * (b[1] - a[1])),
64+ Math.round(a[2] + f * (b[2] - a[2])),
65+ ];
66+}
webapp/src/hooks/useRenderer.tsadded+114−0View file
@@ -0,0 +1,114 @@
1+// Hook owning the render worker: debounced auto-preview whenever the
2+// parameters or the .m sources change, explicit full-quality render, and a
3+// latest-wins request queue.
4+
5+import { useCallback, useEffect, useRef, useState } from 'react';
6+import type { RecorderParams } from '../params';
7+import type { MFile } from '../mfiles';
8+import type { ErrorMessage, ProgressMessage, RenderResponse } from '../worker/renderWorker';
9+
10+export interface RenderError {
11+ message: string;
12+ file?: string;
13+ line?: number;
14+}
15+
16+export interface RenderState {
17+ result: RenderResponse | null;
18+ error: RenderError | null;
19+ busy: boolean;
20+ stage: string;
21+ progress: number; // 0..1 of the current stage
22+}
23+
24+interface Request {
25+ params: RecorderParams;
26+ mfiles: MFile[];
27+ quality: 'preview' | 'full';
28+}
29+
30+export function useRenderer(params: RecorderParams, mfiles: MFile[]) {
31+ const workerRef = useRef<Worker | null>(null);
32+ const nextId = useRef(1);
33+ const inFlight = useRef<number | null>(null);
34+ const queued = useRef<Request | null>(null);
35+ const sendRef = useRef<(r: Request) => void>(() => {});
36+ const [state, setState] = useState<RenderState>({
37+ result: null,
38+ error: null,
39+ busy: false,
40+ stage: '',
41+ progress: 0,
42+ });
43+
44+ const send = useCallback((req: Request) => {
45+ const worker = workerRef.current;
46+ if (!worker) return;
47+ if (inFlight.current !== null) {
48+ queued.current = req; // latest wins
49+ return;
50+ }
51+ const id = nextId.current++;
52+ inFlight.current = id;
53+ setState((s) => ({ ...s, busy: true, stage: 'starting', progress: 0 }));
54+ worker.postMessage({ type: 'render', id, ...req });
55+ }, []);
56+ sendRef.current = send;
57+
58+ useEffect(() => {
59+ const worker = new Worker(new URL('../worker/renderWorker.ts', import.meta.url), {
60+ type: 'module',
61+ });
62+ workerRef.current = worker;
63+
64+ const finish = (update: (s: RenderState) => RenderState) => {
65+ inFlight.current = null;
66+ const next = queued.current;
67+ queued.current = null;
68+ setState((s) => {
69+ const done = update(s);
70+ return next ? { ...done, busy: true, stage: 'queued', progress: 0 } : done;
71+ });
72+ if (next) sendRef.current(next);
73+ };
74+
75+ worker.onmessage = (ev: MessageEvent<RenderResponse | ProgressMessage | ErrorMessage>) => {
76+ const msg = ev.data;
77+ if (msg.type === 'progress') {
78+ setState((s) => ({ ...s, stage: msg.stage, progress: msg.frac }));
79+ return;
80+ }
81+ if (msg.type === 'error') {
82+ // Keep the last good result on screen; the code panel shows the fault.
83+ finish((s) => ({
84+ ...s,
85+ error: { message: msg.message, file: msg.file, line: msg.line },
86+ busy: false,
87+ stage: '',
88+ progress: 0,
89+ }));
90+ return;
91+ }
92+ finish(() => ({ result: msg, error: null, busy: false, stage: '', progress: 0 }));
93+ };
94+ return () => {
95+ worker.terminate();
96+ workerRef.current = null;
97+ inFlight.current = null;
98+ queued.current = null;
99+ };
100+ }, []);
101+
102+ const renderFull = useCallback(
103+ () => send({ params, mfiles, quality: 'full' }),
104+ [params, mfiles, send],
105+ );
106+
107+ // Debounced auto-preview on any change to the parameters or the code.
108+ useEffect(() => {
109+ const t = setTimeout(() => send({ params, mfiles, quality: 'preview' }), 250);
110+ return () => clearTimeout(t);
111+ }, [params, mfiles, send]);
112+
113+ return { ...state, renderFull };
114+}
webapp/src/lib/fft.tsadded+91−0View file
@@ -0,0 +1,91 @@
1+// Minimal iterative radix-2 complex FFT, in-place on separate re/im arrays.
2+
3+export function fft(re: Float64Array, im: Float64Array): void {
4+ const n = re.length;
5+ if ((n & (n - 1)) !== 0) throw new Error('FFT length must be a power of 2');
6+ // bit-reversal permutation
7+ for (let i = 1, j = 0; i < n; i++) {
8+ let bit = n >> 1;
9+ for (; j & bit; bit >>= 1) j ^= bit;
10+ j ^= bit;
11+ if (i < j) {
12+ const tr = re[i]; re[i] = re[j]; re[j] = tr;
13+ const ti = im[i]; im[i] = im[j]; im[j] = ti;
14+ }
15+ }
16+ for (let len = 2; len <= n; len <<= 1) {
17+ const ang = (-2 * Math.PI) / len;
18+ const wr = Math.cos(ang);
19+ const wi = Math.sin(ang);
20+ for (let i = 0; i < n; i += len) {
21+ let cwr = 1, cwi = 0;
22+ for (let k = 0; k < len / 2; k++) {
23+ const ur = re[i + k], ui = im[i + k];
24+ const vr = re[i + k + len / 2] * cwr - im[i + k + len / 2] * cwi;
25+ const vi = re[i + k + len / 2] * cwi + im[i + k + len / 2] * cwr;
26+ re[i + k] = ur + vr;
27+ im[i + k] = ui + vi;
28+ re[i + k + len / 2] = ur - vr;
29+ im[i + k + len / 2] = ui - vi;
30+ const nwr = cwr * wr - cwi * wi;
31+ cwi = cwr * wi + cwi * wr;
32+ cwr = nwr;
33+ }
34+ }
35+ }
36+}
37+
38+// Magnitude spectrum of a real signal with a Hann window, zero-padded to the
39+// next power of two (at least minSize). Returns bins 0..N/2.
40+export function magnitudeSpectrum(x: Float32Array | Float64Array, minSize = 16384): {
41+ mag: Float64Array;
42+ size: number;
43+} {
44+ let size = minSize;
45+ while (size < x.length) size <<= 1;
46+ const re = new Float64Array(size);
47+ const im = new Float64Array(size);
48+ const n = x.length;
49+ for (let i = 0; i < n; i++) {
50+ const w = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (n - 1)));
51+ re[i] = x[i] * w;
52+ }
53+ fft(re, im);
54+ const half = size / 2;
55+ const mag = new Float64Array(half + 1);
56+ for (let i = 0; i <= half; i++) mag[i] = Math.hypot(re[i], im[i]);
57+ return { mag, size };
58+}
59+
60+// Simple STFT magnitude in dB for the spectrogram, frames x bins.
61+export function stft(
62+ x: Float32Array,
63+ win = 512,
64+ hop = 128,
65+): { frames: number; bins: number; db: Float32Array; floorDb: number } {
66+ const nFrames = Math.max(1, Math.floor((x.length - win) / hop) + 1);
67+ const bins = win / 2 + 1;
68+ const db = new Float32Array(nFrames * bins);
69+ const re = new Float64Array(win);
70+ const im = new Float64Array(win);
71+ const hann = new Float64Array(win);
72+ for (let i = 0; i < win; i++) hann[i] = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (win - 1)));
73+ let maxDb = -Infinity;
74+ for (let f = 0; f < nFrames; f++) {
75+ const off = f * hop;
76+ for (let i = 0; i < win; i++) {
77+ re[i] = (off + i < x.length ? x[off + i] : 0) * hann[i];
78+ im[i] = 0;
79+ }
80+ fft(re, im);
81+ for (let b = 0; b < bins; b++) {
82+ const m = Math.hypot(re[b], im[b]);
83+ const v = 20 * Math.log10(m + 1e-12);
84+ db[f * bins + b] = v;
85+ if (v > maxDb) maxDb = v;
86+ }
87+ }
88+ // normalize so 0 dB is the frame-bin maximum
89+ for (let i = 0; i < db.length; i++) db[i] -= maxDb;
90+ return { frames: nFrames, bins, db, floorDb: -90 };
91+}
webapp/src/lib/instrument.tsadded+25−0View file
@@ -0,0 +1,25 @@
1+// The tone hole lattice and fingering chart, mirrored from
2+// ../matlab/recorder_geometry.m for the 3D scene. Hole 0 is the thumb hole
3+// (back of the instrument); 1..7 are the finger holes from the top down.
4+// Positions are fractions of the bore length below the window.
5+
6+export const HOLE_FRAC = [0.3699, 0.4271, 0.464, 0.5292, 0.6204, 0.6519, 0.7424, 0.828];
7+export const HOLE_RAD = [0.0027, 0.0032, 0.003, 0.0029, 0.0033, 0.0026, 0.0032, 0.0025];
8+
9+// 0 closed, 1 open, 0.5 pinched; rows are the fingering chart notes.
10+export const FINGERINGS: number[][] = [
11+ [0, 0, 0, 0, 0, 0, 0, 0], // C5
12+ [0, 0, 0, 0, 0, 0, 0, 1], // D5
13+ [0, 0, 0, 0, 0, 0, 1, 1], // E5
14+ [0, 0, 0, 0, 0, 1, 0, 0], // F5 (fork)
15+ [0, 0, 0, 0, 1, 1, 1, 1], // G5
16+ [0, 0, 0, 1, 1, 1, 1, 1], // A5
17+ [0, 0, 1, 1, 1, 1, 1, 1], // B5
18+ [0, 1, 0, 1, 1, 1, 1, 1], // C6
19+ [1, 0, 0, 1, 1, 1, 1, 1], // C#6
20+ [1, 1, 0, 1, 1, 1, 1, 1], // D6
21+ [0.5, 0, 0, 0, 0, 0, 1, 1], // E6
22+ [0.5, 0, 0, 0, 0, 1, 0, 1], // F6
23+ [0.5, 0, 0, 0, 1, 1, 1, 1], // G6
24+ [0.5, 0, 0, 1, 1, 1, 1, 1], // A6
25+];
webapp/src/lib/wav.tsadded+40−0View file
@@ -0,0 +1,40 @@
1+// Mono PCM WAV encoding (16- or 24-bit).
2+
3+export function encodeWav(y: Float32Array, sampleRate: number, bitDepth: 16 | 24): Blob {
4+ const bytesPerSample = bitDepth / 8;
5+ const dataSize = y.length * bytesPerSample;
6+ const buf = new ArrayBuffer(44 + dataSize);
7+ const dv = new DataView(buf);
8+ const writeStr = (off: number, s: string) => {
9+ for (let i = 0; i < s.length; i++) dv.setUint8(off + i, s.charCodeAt(i));
10+ };
11+ writeStr(0, 'RIFF');
12+ dv.setUint32(4, 36 + dataSize, true);
13+ writeStr(8, 'WAVE');
14+ writeStr(12, 'fmt ');
15+ dv.setUint32(16, 16, true);
16+ dv.setUint16(20, 1, true); // PCM
17+ dv.setUint16(22, 1, true); // mono
18+ dv.setUint32(24, sampleRate, true);
19+ dv.setUint32(28, sampleRate * bytesPerSample, true);
20+ dv.setUint16(32, bytesPerSample, true);
21+ dv.setUint16(34, bitDepth, true);
22+ writeStr(36, 'data');
23+ dv.setUint32(40, dataSize, true);
24+
25+ let off = 44;
26+ for (let i = 0; i < y.length; i++) {
27+ const v = Math.max(-1, Math.min(1, y[i]));
28+ if (bitDepth === 16) {
29+ dv.setInt16(off, Math.round(v * 32767), true);
30+ off += 2;
31+ } else {
32+ const iv = Math.round(v * 8388607);
33+ dv.setUint8(off, iv & 0xff);
34+ dv.setUint8(off + 1, (iv >> 8) & 0xff);
35+ dv.setUint8(off + 2, (iv >> 16) & 0xff);
36+ off += 3;
37+ }
38+ }
39+ return new Blob([buf], { type: 'audio/wav' });
40+}
webapp/src/main.tsxadded+10−0View file
@@ -0,0 +1,10 @@
1+import React from 'react';
2+import { createRoot } from 'react-dom/client';
3+import App from './App';
4+import './style.css';
5+
6+createRoot(document.getElementById('root')!).render(
7+ <React.StrictMode>
8+ <App />
9+ </React.StrictMode>,
10+);
webapp/src/mfiles.tsadded+32−0View file
@@ -0,0 +1,32 @@
1+// The MATLAB implementation, loaded from ../matlab at build time. The app
2+// keeps an editable in-memory copy of these sources (see App.tsx); edits feed
3+// straight into the render worker and are deliberately not persisted, so a
4+// reload restores the shipped code.
5+
6+export type { MFile } from './mfilesCore';
7+export { hashMFiles } from './mfilesCore';
8+import type { MFile } from './mfilesCore';
9+
10+const modules = import.meta.glob('../../matlab/*.m', {
11+ eager: true,
12+ query: '?raw',
13+ import: 'default',
14+}) as Record<string, string>;
15+
16+export const DEFAULT_MFILES: MFile[] = Object.entries(modules)
17+ .map(([path, source]) => ({ name: path.split('/').pop()!, source }))
18+ .sort((a, b) => a.name.localeCompare(b.name));
19+
20+// Files the app itself executes. The others are entry points for running the
21+// model from the numbl CLI, and editing them changes nothing in the browser.
22+export const EXECUTED_BY_APP = new Set([
23+ 'recorder_geometry.m',
24+ 'recorder_bore.m',
25+ 'recorder_basis.m',
26+ 'recorder_notes.m',
27+ 'recorder_jet.m',
28+ 'recorder_rng.m',
29+ 'recorder_radiate.m',
30+ 'recorder_output.m',
31+]);
32+
webapp/src/mfilesCore.tsadded+21−0View file
@@ -0,0 +1,21 @@
1+// The MFile shape and the content hash, kept free of Vite-specific APIs so
2+// the render worker (and the Node harnesses) can import them without
3+// touching import.meta.glob.
4+
5+export interface MFile {
6+ name: string;
7+ source: string;
8+}
9+
10+// Cheap content hash (FNV-1a) so the worker can key its caches on the code.
11+export function hashMFiles(files: MFile[]): string {
12+ let h = 0x811c9dc5;
13+ for (const f of files) {
14+ const s = f.name + ' ' + f.source + '';
15+ for (let i = 0; i < s.length; i++) {
16+ h ^= s.charCodeAt(i);
17+ h = Math.imul(h, 0x01000193);
18+ }
19+ }
20+ return (h >>> 0).toString(36);
21+}
webapp/src/params.tsadded+284−0View file
@@ -0,0 +1,284 @@
1+// Parameter definitions for the recorder model, mirroring recorder_defaults.m.
2+// Each parameter is assigned a tier that determines how much recomputation a
3+// change requires:
4+// basis -> bore transmission-line solve + modal fit (per note), then jet,
5+// radiation, and output
6+// jet -> the time-domain jet integration (the expensive stage of this
7+// model), then radiation and output
8+// radiate -> radiation + output only (moving the listener is cheap)
9+// output -> output stage only
10+//
11+// Note that the render duration is not a parameter: it is
12+// note_count * note_dur_s + (note_count - 1) * gap_ms + tail_s.
13+
14+export interface RecorderParams {
15+ // Bore and tone holes
16+ bore_length_m: number;
17+ bore_r_head_m: number;
18+ bore_r_foot_m: number;
19+ wall_m: number;
20+ hole_scale: number;
21+ thumb_pinch: number;
22+ pinch_loss: number;
23+ // Window, windway, and jet
24+ cutup_m: number;
25+ window_width_m: number;
26+ flue_h_m: number;
27+ flue_w_m: number;
28+ labium_offset_m: number;
29+ jet_gain: number;
30+ jet_growth: number;
31+ jet_conv: number;
32+ jet_stcut: number;
33+ vortex_loss: number;
34+ noise_level: number;
35+ direct_gain: number;
36+ // Breath and gesture
37+ breath_pa: number;
38+ breath_track: number;
39+ attack_ms: number;
40+ release_ms: number;
41+ vibrato_hz: number;
42+ vibrato_depth: number;
43+ note_count: number;
44+ note1: number;
45+ note2: number;
46+ note3: number;
47+ note_dur_s: number;
48+ gap_ms: number;
49+ seed: number;
50+ // Modal fit
51+ max_mode_hz: number;
52+ max_modes: number;
53+ grid_hz: number;
54+ // Radiation / listener
55+ listener_x_m: number;
56+ listener_y_m: number;
57+ listener_z_m: number;
58+ speed_of_sound: number;
59+ air_density: number;
60+ // Rendering
61+ sample_rate: number;
62+ oversample: number;
63+ tail_s: number;
64+ pre_delay_ms: number;
65+ normalize: boolean;
66+ peak_dbfs: number;
67+ full_scale_pa: number;
68+ bit_depth: number;
69+}
70+
71+export const DEFAULTS: RecorderParams = {
72+ bore_length_m: 0.3015,
73+ bore_r_head_m: 0.0067,
74+ bore_r_foot_m: 0.0052,
75+ wall_m: 0.003,
76+ hole_scale: 1.0,
77+ thumb_pinch: 0.35,
78+ pinch_loss: 0.5,
79+ cutup_m: 0.0045,
80+ window_width_m: 0.01,
81+ flue_h_m: 0.001,
82+ flue_w_m: 0.01,
83+ labium_offset_m: 0.0004,
84+ jet_gain: 1.0,
85+ jet_growth: 0.4,
86+ jet_conv: 0.55,
87+ jet_stcut: 0.45,
88+ vortex_loss: 1.0,
89+ noise_level: 0.03,
90+ direct_gain: 1.0,
91+ breath_pa: 350.0,
92+ breath_track: 1.6,
93+ attack_ms: 12.0,
94+ release_ms: 25.0,
95+ vibrato_hz: 5.0,
96+ vibrato_depth: 0.0,
97+ note_count: 1,
98+ note1: 5,
99+ note2: 6,
100+ note3: 7,
101+ note_dur_s: 1.6,
102+ gap_ms: 60.0,
103+ seed: 1,
104+ max_mode_hz: 10000.0,
105+ max_modes: 12,
106+ grid_hz: 2.5,
107+ listener_x_m: 0.25,
108+ listener_y_m: -0.1,
109+ listener_z_m: 0.35,
110+ speed_of_sound: 343.0,
111+ air_density: 1.204,
112+ sample_rate: 48000,
113+ oversample: 2,
114+ tail_s: 0.4,
115+ pre_delay_ms: 3.0,
116+ normalize: true,
117+ peak_dbfs: -3.0,
118+ full_scale_pa: 0.5,
119+ bit_depth: 24,
120+};
121+
122+export const NOTE_NAMES = [
123+ 'C5', 'D5', 'E5', 'F5', 'G5', 'A5', 'B5', 'C6', 'C♯6', 'D6',
124+ 'E6', 'F6', 'G6', 'A6',
125+];
126+
127+export type Tier = 'basis' | 'jet' | 'radiate' | 'output';
128+
129+// Keys whose change invalidates the cached bore basis of every note. The
130+// jet, the breath, and the listener are deliberately absent: changing how
131+// the instrument is blown or where the microphone stands reuses the bore.
132+export const BASIS_KEYS: (keyof RecorderParams)[] = [
133+ 'bore_length_m', 'bore_r_head_m', 'bore_r_foot_m', 'wall_m',
134+ 'hole_scale', 'thumb_pinch', 'pinch_loss',
135+ 'cutup_m', 'window_width_m',
136+ 'max_mode_hz', 'max_modes', 'grid_hz',
137+ 'speed_of_sound', 'air_density',
138+];
139+
140+// Keys that re-run the time-domain jet integration (per note).
141+export const JET_KEYS: (keyof RecorderParams)[] = [
142+ 'flue_h_m', 'flue_w_m', 'labium_offset_m',
143+ 'jet_gain', 'jet_growth', 'jet_conv', 'jet_stcut', 'vortex_loss',
144+ 'noise_level',
145+ 'breath_pa', 'breath_track', 'attack_ms', 'release_ms',
146+ 'vibrato_hz', 'vibrato_depth',
147+ 'note_dur_s', 'tail_s', 'seed',
148+ 'sample_rate',
149+];
150+
151+export const RADIATE_KEYS: (keyof RecorderParams)[] = [
152+ 'listener_x_m', 'listener_y_m', 'listener_z_m', 'direct_gain', 'gap_ms',
153+];
154+
155+export const OUTPUT_KEYS: (keyof RecorderParams)[] = [
156+ 'pre_delay_ms', 'normalize', 'peak_dbfs', 'full_scale_pa',
157+];
158+
159+export interface ParamSpec {
160+ key: keyof RecorderParams;
161+ label: string;
162+ unit?: string;
163+ min: number;
164+ max: number;
165+ step?: number;
166+ log?: boolean;
167+ integer?: boolean;
168+ note?: boolean; // render as a note-name dropdown
169+ help?: string;
170+}
171+
172+export interface ParamGroup {
173+ name: string;
174+ specs: ParamSpec[];
175+}
176+
177+export const PARAM_GROUPS: ParamGroup[] = [
178+ {
179+ name: 'Notes',
180+ specs: [
181+ { key: 'note1', label: 'Note', min: 1, max: 14, step: 1, integer: true, note: true, help: 'Fingering from the baroque chart. E6 and above are pinched-thumb second-register fingerings.' },
182+ { key: 'note_count', label: 'Notes played', min: 1, max: 3, step: 1, integer: true, help: '1 plays a single note; 2 or 3 play a tongued figure.' },
183+ { key: 'note2', label: 'Second note', min: 1, max: 14, step: 1, integer: true, note: true },
184+ { key: 'note3', label: 'Third note', min: 1, max: 14, step: 1, integer: true, note: true },
185+ { key: 'note_dur_s', label: 'Note length', unit: 's', min: 0.2, max: 4.0, step: 0.05 },
186+ { key: 'gap_ms', label: 'Gap', unit: 'ms', min: 0, max: 300, step: 5, help: 'Silence between tongued notes.' },
187+ ],
188+ },
189+ {
190+ name: 'Breath',
191+ specs: [
192+ { key: 'breath_pa', label: 'Breath pressure', unit: 'Pa', min: 80, max: 3000, log: true, help: 'Blowing pressure for C5; ~350 Pa speaks the first register. Blow harder and the note sharpens, then overblows to the octave.' },
193+ { key: 'breath_track', label: 'Breath tracking', min: 0, max: 2.5, step: 0.05, help: 'How much harder higher notes are blown, as (f/523 Hz)^track. 0 blows every fingering identically; high notes then underblow, as on the instrument.' },
194+ { key: 'attack_ms', label: 'Attack', unit: 'ms', min: 2, max: 60, step: 1, help: 'Breath rise time. The attack chiff is the jet sweeping through its lock range as the pressure builds.' },
195+ { key: 'release_ms', label: 'Release', unit: 'ms', min: 5, max: 120, step: 1 },
196+ { key: 'vibrato_hz', label: 'Vibrato rate', unit: 'Hz', min: 0.5, max: 8, step: 0.1 },
197+ { key: 'vibrato_depth', label: 'Vibrato depth', min: 0, max: 0.2, step: 0.005, help: 'Breath pressure modulation; a recorder vibrato is breath, not pitch, so it is mostly loudness and a little pitch.' },
198+ { key: 'noise_level', label: 'Breathiness', min: 0, max: 0.15, step: 0.005, help: 'Turbulence displacement noise on the jet.' },
199+ { key: 'seed', label: 'Noise seed', min: 1, max: 9999, step: 1, integer: true, help: 'Any change redraws the turbulence; the same seed reproduces it exactly.' },
200+ ],
201+ },
202+ {
203+ name: 'Windway & jet',
204+ specs: [
205+ { key: 'cutup_m', label: 'Cutup', unit: 'm', min: 0.003, max: 0.007, step: 0.0001, help: 'Flue exit to labium. Longer needs more breath and darkens the tone.' },
206+ { key: 'flue_h_m', label: 'Flue height', unit: 'm', min: 0.0005, max: 0.002, step: 0.00005, help: 'Windway exit slot height; sets the jet thickness.' },
207+ { key: 'flue_w_m', label: 'Flue width', unit: 'm', min: 0.006, max: 0.014, step: 0.0005 },
208+ { key: 'labium_offset_m', label: 'Labium offset', unit: 'm', min: -0.001, max: 0.0015, step: 0.00005, help: 'Labium position relative to the jet centerline; sets the even-harmonic content and the ease of speech.' },
209+ { key: 'jet_gain', label: 'Jet receptivity', min: 0.3, max: 3.0, step: 0.05, help: 'Trim on how strongly the mouth field deflects the jet. Too low and notes fail to speak.' },
210+ { key: 'jet_growth', label: 'Jet growth', min: 0.2, max: 0.7, step: 0.01, help: 'Instability amplification per flue height of travel.' },
211+ { key: 'jet_conv', label: 'Jet wave speed', min: 0.3, max: 0.7, step: 0.01, help: 'Displacement wave speed as a fraction of the jet velocity; sets which breath pressure centers each note.' },
212+ { key: 'jet_stcut', label: 'Jet cutoff', min: 0.25, max: 0.7, step: 0.01, help: 'Instability cutoff (Strouhal-like). Lower favors low registers; raise it and low fingerings squeak.' },
213+ { key: 'vortex_loss', label: 'Vortex loss', min: 0, max: 3, step: 0.05, help: 'Nonlinear loss at the window; with the jet saturation it caps the amplitude.' },
214+ ],
215+ },
216+ {
217+ name: 'Bore & holes',
218+ specs: [
219+ { key: 'bore_length_m', label: 'Bore length', unit: 'm', min: 0.2, max: 0.45, step: 0.001, help: 'Window to foot. The hole lattice scales with it, so this transposes the instrument.' },
220+ { key: 'bore_r_head_m', label: 'Bore radius (head)', unit: 'm', min: 0.004, max: 0.01, step: 0.0001 },
221+ { key: 'bore_r_foot_m', label: 'Bore radius (foot)', unit: 'm', min: 0.003, max: 0.009, step: 0.0001, help: 'The taper tunes the octaves against each other.' },
222+ { key: 'wall_m', label: 'Wall thickness', unit: 'm', min: 0.002, max: 0.005, step: 0.0001, help: 'Tone hole chimney height.' },
223+ { key: 'hole_scale', label: 'Hole size scale', min: 0.7, max: 1.3, step: 0.01, help: 'Scales every tone hole radius. Smaller holes flatten and darken the scale.' },
224+ { key: 'thumb_pinch', label: 'Thumb pinch', min: 0.15, max: 0.7, step: 0.01, help: 'Effective radius fraction of the pinched thumb hole (the second-register vent).' },
225+ { key: 'pinch_loss', label: 'Pinch loss', min: 0, max: 2, step: 0.05, help: 'Viscous resistance of the thumb crack; damps the fundamental so the second register speaks.' },
226+ ],
227+ },
228+ {
229+ name: 'Listener',
230+ specs: [
231+ { key: 'listener_x_m', label: 'Listener x', unit: 'm', min: -2.0, max: 2.0, step: 0.005, help: 'To the player’s side. Radiation-only: no re-solve needed.' },
232+ { key: 'listener_y_m', label: 'Listener y', unit: 'm', min: -2.0, max: 2.0, step: 0.005, help: 'Along the instrument; the foot points down (-y).' },
233+ { key: 'listener_z_m', label: 'Listener z', unit: 'm', min: -2.0, max: 2.0, step: 0.005, help: 'In front of the instrument, the side the holes open toward.' },
234+ { key: 'direct_gain', label: 'Direct breath gain', min: 0.0, max: 3.0, step: 0.05, help: 'Broadband monopole from the window: breath noise and attack transient. 0 leaves only the resonant modes.' },
235+ ],
236+ },
237+ {
238+ name: 'Output',
239+ specs: [
240+ { key: 'tail_s', label: 'Tail', unit: 's', min: 0.1, max: 2.0, step: 0.05, help: 'Rendered after the last release. A recorder stops quickly; half a second is plenty.' },
241+ { key: 'pre_delay_ms', label: 'Pre-delay', unit: 'ms', min: 0, max: 50, step: 0.5 },
242+ { key: 'peak_dbfs', label: 'Peak level', unit: 'dBFS', min: -24, max: 0, step: 0.5, help: 'Used when normalize is on.' },
243+ { key: 'full_scale_pa', label: 'Full scale', unit: 'Pa', min: 0.05, max: 10, step: 0.05, help: 'Used when normalize is off: this many pascals maps to 0 dBFS.' },
244+ ],
245+ },
246+ {
247+ name: 'Advanced / quality',
248+ specs: [
249+ { key: 'max_mode_hz', label: 'Max mode freq', unit: 'Hz', min: 4000, max: 14000, step: 250, help: 'Bore resonances above this are dropped; their broadband remainder rides on the direct term.' },
250+ { key: 'max_modes', label: 'Max modes', min: 4, max: 16, step: 1, integer: true },
251+ { key: 'grid_hz', label: 'Impedance grid', unit: 'Hz', min: 1, max: 5, step: 0.25, help: 'Frequency step of the transmission-line solve the modal fit reads.' },
252+ { key: 'oversample', label: 'Oversample', min: 1, max: 4, step: 1, integer: true, help: 'Jet integration rate as a multiple of the sample rate; 2 for render, 1 for preview.' },
253+ { key: 'speed_of_sound', label: 'Speed of sound', unit: 'm/s', min: 300, max: 400, step: 1 },
254+ { key: 'air_density', label: 'Air density', unit: 'kg/m³', min: 0.9, max: 1.5, step: 0.01 },
255+ ],
256+ },
257+];
258+
259+const STORAGE_KEY = 'recorder-webapp-params-v1';
260+
261+export function loadParams(): RecorderParams {
262+ try {
263+ const raw = localStorage.getItem(STORAGE_KEY);
264+ if (!raw) return { ...DEFAULTS };
265+ const parsed = JSON.parse(raw);
266+ const out = { ...DEFAULTS };
267+ for (const k of Object.keys(DEFAULTS) as (keyof RecorderParams)[]) {
268+ if (typeof parsed[k] === typeof DEFAULTS[k]) {
269+ (out as Record<string, unknown>)[k] = parsed[k];
270+ }
271+ }
272+ return out;
273+ } catch {
274+ return { ...DEFAULTS };
275+ }
276+}
277+
278+export function saveParams(p: RecorderParams): void {
279+ try {
280+ localStorage.setItem(STORAGE_KEY, JSON.stringify(p));
281+ } catch {
282+ // storage unavailable; persistence is a convenience only
283+ }
284+}
webapp/src/selftest.tsadded+26−0View file
@@ -0,0 +1,26 @@
1+// Runs the MATLAB validation script (../matlab/recorder_validate.m) through
2+// numbl in Node: closed-form identities of the passive acoustics, emergent
3+// register selection and overblowing, exact reproducibility of the seeded
4+// turbulence, and the pinned end-to-end diagnostics.
5+// Run with: npm run selftest
6+
7+import { readFileSync, readdirSync } from 'fs';
8+import { join, dirname } from 'path';
9+import { fileURLToPath } from 'url';
10+import { executeCode } from 'numbl';
11+
12+const matlabDir = join(dirname(fileURLToPath(import.meta.url)), '..', '..', 'matlab');
13+const workspaceFiles = readdirSync(matlabDir)
14+ .filter((f) => f.endsWith('.m'))
15+ .map((f) => ({ name: f, source: readFileSync(join(matlabDir, f), 'utf8') }));
16+
17+const main = workspaceFiles.find((f) => f.name === 'recorder_validate.m')!;
18+const t0 = performance.now();
19+const result = executeCode(main.source, { displayResults: false }, workspaceFiles, 'recorder_validate.m');
20+const text = result.output.join('');
21+process.stdout.write(text);
22+console.log(`total time : ${((performance.now() - t0) / 1000).toFixed(2)} s`);
23+
24+if (text.includes('MISMATCH') || text.includes('FAILED')) {
25+ process.exit(1);
26+}
webapp/src/style.cssadded+358−0View file
@@ -0,0 +1,358 @@
1+:root {
2+ --bg: #14161a;
3+ --panel: #1d2026;
4+ --panel2: #23272f;
5+ --border: #31363f;
6+ --text: #d6dae2;
7+ --muted: #8b93a1;
8+ --accent: #e8a33d;
9+ --accent2: #4dd0e1;
10+ --good: #6fbf73;
11+ --warn: #e57373;
12+ font-family: system-ui, -apple-system, 'Segoe UI', Roboto, sans-serif;
13+}
14+
15+* { box-sizing: border-box; }
16+
17+html, body, #root {
18+ margin: 0;
19+ height: 100%;
20+ background: var(--bg);
21+ color: var(--text);
22+}
23+
24+.app {
25+ display: flex;
26+ flex-direction: column;
27+ height: 100vh;
28+}
29+
30+.topbar {
31+ display: flex;
32+ align-items: center;
33+ gap: 10px;
34+ padding: 8px 14px;
35+ background: var(--panel);
36+ border-bottom: 1px solid var(--border);
37+ flex-wrap: wrap;
38+}
39+
40+.topbar h1 {
41+ font-size: 16px;
42+ margin: 0 12px 0 0;
43+ font-weight: 600;
44+ color: var(--accent);
45+ white-space: nowrap;
46+}
47+
48+button {
49+ background: var(--panel2);
50+ color: var(--text);
51+ border: 1px solid var(--border);
52+ border-radius: 6px;
53+ padding: 6px 14px;
54+ font-size: 13px;
55+ cursor: pointer;
56+}
57+button:hover { border-color: var(--muted); }
58+button:disabled { opacity: 0.45; cursor: default; }
59+button.primary {
60+ background: var(--accent);
61+ color: #1a1204;
62+ border-color: var(--accent);
63+ font-weight: 600;
64+}
65+button.play {
66+ background: var(--good);
67+ color: #06210a;
68+ border-color: var(--good);
69+ font-weight: 600;
70+}
71+button.play.preview {
72+ background: transparent;
73+ color: var(--accent);
74+ border-color: var(--accent);
75+}
76+button.loop.active {
77+ background: var(--accent2);
78+ color: #062a2f;
79+ border-color: var(--accent2);
80+ font-weight: 600;
81+}
82+button.code { font-family: ui-monospace, SFMono-Regular, Menlo, monospace; }
83+button.code.active { border-color: var(--accent2); color: var(--accent2); }
84+
85+.status {
86+ margin-left: auto;
87+ font-size: 12px;
88+ color: var(--muted);
89+ display: flex;
90+ align-items: center;
91+ gap: 8px;
92+ min-width: 180px;
93+ justify-content: flex-end;
94+}
95+
96+.badge {
97+ font-size: 11px;
98+ border-radius: 10px;
99+ padding: 2px 8px;
100+ border: 1px solid var(--border);
101+}
102+.badge.preview { color: var(--accent); border-color: var(--accent); }
103+.badge.full { color: var(--good); border-color: var(--good); }
104+.badge.err { color: var(--warn); border-color: var(--warn); }
105+
106+.progressbar {
107+ width: 90px;
108+ height: 6px;
109+ background: var(--panel2);
110+ border-radius: 3px;
111+ overflow: hidden;
112+}
113+.progressbar > div {
114+ height: 100%;
115+ background: var(--accent2);
116+ transition: width 0.1s linear;
117+}
118+
119+.main {
120+ display: flex;
121+ flex: 1;
122+ min-height: 0;
123+}
124+
125+.sidebar {
126+ width: 340px;
127+ min-width: 340px;
128+ overflow-y: auto;
129+ background: var(--panel);
130+ border-right: 1px solid var(--border);
131+ padding: 8px;
132+}
133+
134+.content {
135+ flex: 1;
136+ overflow-y: auto;
137+ padding: 10px;
138+ display: grid;
139+ grid-template-columns: 1fr 1fr;
140+ gap: 10px;
141+ align-content: start;
142+}
143+
144+.card {
145+ background: var(--panel);
146+ border: 1px solid var(--border);
147+ border-radius: 8px;
148+ padding: 8px 10px 10px;
149+ min-width: 0;
150+}
151+.card.wide { grid-column: 1 / -1; }
152+.card h2 {
153+ font-size: 12px;
154+ text-transform: uppercase;
155+ letter-spacing: 0.06em;
156+ color: var(--muted);
157+ margin: 0 0 6px;
158+ font-weight: 600;
159+}
160+.card canvas { display: block; width: 100%; }
161+
162+.scene3d { height: 380px; position: relative; }
163+.scene3d canvas { height: 100%; border-radius: 6px; }
164+.scene-hint {
165+ position: absolute;
166+ bottom: 10px;
167+ right: 14px;
168+ font-size: 11px;
169+ color: var(--muted);
170+ pointer-events: none;
171+}
172+
173+/* Parameter panel */
174+details.group {
175+ border: 1px solid var(--border);
176+ border-radius: 6px;
177+ margin-bottom: 6px;
178+ background: var(--panel2);
179+}
180+details.group summary {
181+ cursor: pointer;
182+ padding: 6px 10px;
183+ font-size: 13px;
184+ font-weight: 600;
185+ user-select: none;
186+}
187+details.group[open] summary { border-bottom: 1px solid var(--border); }
188+.group-body { padding: 6px 10px 8px; }
189+
190+.param-row { margin-bottom: 8px; }
191+.param-label {
192+ display: flex;
193+ justify-content: space-between;
194+ font-size: 12px;
195+ margin-bottom: 2px;
196+ align-items: baseline;
197+}
198+.param-label .name { color: var(--text); }
199+.param-controls { display: flex; gap: 8px; align-items: center; }
200+.param-controls input[type='range'] { flex: 1; accent-color: var(--accent2); }
201+.param-controls input[type='number'] {
202+ width: 88px;
203+ background: var(--bg);
204+ color: var(--text);
205+ border: 1px solid var(--border);
206+ border-radius: 4px;
207+ padding: 3px 5px;
208+ font-size: 12px;
209+}
210+.param-controls select, .param-controls input[type='checkbox'] {
211+ background: var(--bg);
212+ color: var(--text);
213+ border: 1px solid var(--border);
214+ border-radius: 4px;
215+ font-size: 12px;
216+}
217+.param-help { font-size: 11px; color: var(--muted); margin-top: 1px; }
218+
219+/* Diagnostics */
220+.diag table { border-collapse: collapse; font-size: 12px; width: 100%; }
221+.diag td, .diag th {
222+ padding: 2px 8px 2px 0;
223+ text-align: left;
224+ color: var(--text);
225+ font-weight: 400;
226+}
227+.diag th { color: var(--muted); white-space: nowrap; }
228+.diag .partials { font-variant-numeric: tabular-nums; }
229+.warning { color: var(--warn); font-size: 12px; margin-top: 4px; }
230+.cachenote { color: var(--muted); font-size: 11px; margin-top: 6px; }
231+
232+.mode-list { display: flex; gap: 6px; flex-wrap: wrap; margin-bottom: 6px; }
233+.mode-list button { padding: 3px 10px; font-size: 12px; }
234+.mode-list button.active { border-color: var(--accent2); color: var(--accent2); }
235+
236+/* Code panel */
237+.modal-backdrop {
238+ position: fixed;
239+ inset: 0;
240+ background: rgba(8, 10, 13, 0.72);
241+ display: flex;
242+ align-items: center;
243+ justify-content: center;
244+ z-index: 50;
245+}
246+
247+.code-panel {
248+ width: min(1180px, 94vw);
249+ height: min(860px, 92vh);
250+ background: var(--panel);
251+ border: 1px solid var(--border);
252+ border-radius: 10px;
253+ display: flex;
254+ flex-direction: column;
255+ overflow: hidden;
256+}
257+
258+.code-header {
259+ display: flex;
260+ align-items: center;
261+ gap: 10px;
262+ padding: 9px 12px;
263+ border-bottom: 1px solid var(--border);
264+}
265+.code-header h2 {
266+ margin: 0;
267+ font-size: 13px;
268+ text-transform: uppercase;
269+ letter-spacing: 0.06em;
270+ color: var(--accent);
271+}
272+.code-header button:last-child { margin-left: 0; }
273+.code-note {
274+ font-size: 11.5px;
275+ color: var(--muted);
276+ margin-right: auto;
277+}
278+
279+.code-body { display: flex; flex: 1; min-height: 0; }
280+
281+.code-files {
282+ width: 236px;
283+ min-width: 236px;
284+ border-right: 1px solid var(--border);
285+ overflow-y: auto;
286+ padding: 6px;
287+}
288+.code-files button {
289+ display: flex;
290+ align-items: center;
291+ justify-content: space-between;
292+ gap: 6px;
293+ width: 100%;
294+ text-align: left;
295+ background: transparent;
296+ border: 1px solid transparent;
297+ border-radius: 5px;
298+ padding: 5px 8px;
299+ margin-bottom: 2px;
300+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
301+ font-size: 12px;
302+}
303+.code-files button:hover { background: var(--panel2); }
304+.code-files button.active {
305+ background: var(--panel2);
306+ border-color: var(--accent2);
307+ color: var(--accent2);
308+}
309+.code-files .fname { display: flex; align-items: center; gap: 5px; }
310+.code-files .dot {
311+ width: 6px;
312+ height: 6px;
313+ border-radius: 50%;
314+ background: var(--accent);
315+}
316+.code-files .ftag {
317+ font-size: 9.5px;
318+ color: var(--muted);
319+ text-transform: uppercase;
320+ letter-spacing: 0.05em;
321+ font-family: system-ui, sans-serif;
322+}
323+
324+.code-main { flex: 1; min-width: 0; display: flex; flex-direction: column; }
325+.code-toolbar {
326+ display: flex;
327+ align-items: center;
328+ gap: 10px;
329+ padding: 6px 10px;
330+ border-bottom: 1px solid var(--border);
331+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
332+ font-size: 12px;
333+ color: var(--text);
334+}
335+.code-hint {
336+ font-size: 11px;
337+ color: var(--muted);
338+ font-family: system-ui, sans-serif;
339+}
340+.code-editor { flex: 1; min-height: 0; overflow: hidden; }
341+
342+.code-error {
343+ border-top: 1px solid var(--warn);
344+ background: #2a1a1a;
345+ color: #ffb4a9;
346+ padding: 7px 12px;
347+ font-size: 12px;
348+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
349+}
350+
351+.error-card { border-color: var(--warn); }
352+.error-card h2 { color: var(--warn); }
353+.error-card div { font-size: 12px; }
354+.error-card strong { font-family: ui-monospace, SFMono-Regular, Menlo, monospace; }
355+
356+@media (max-width: 1100px) {
357+ .content { grid-template-columns: 1fr; }
358+}
webapp/src/vite-env.d.tsadded+1−0View file
@@ -0,0 +1 @@
1+/// <reference types="vite/client" />
webapp/src/worker/renderWorker.tsadded+436−0View file
@@ -0,0 +1,436 @@
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).
15+
16+import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
17+import { hashMFiles, type MFile } from '../mfilesCore';
18+import { BASIS_KEYS, JET_KEYS, type RecorderParams } from '../params';
19+
20+export interface RenderRequest {
21+ type: 'render';
22+ id: number;
23+ params: RecorderParams;
24+ quality: 'preview' | 'full';
25+ mfiles: MFile[];
26+}
27+
28+export 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+}
63+
64+export interface ProgressMessage {
65+ type: 'progress';
66+ id: number;
67+ stage: 'basis' | 'jet' | 'radiate';
68+ frac: number;
69+}
70+
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.
73+export interface ErrorMessage {
74+ type: 'error';
75+ id: number;
76+ message: string;
77+ file?: string;
78+ line?: number;
79+}
80+
81+// Basis fields that cross from the basis stage into the jet and radiate
82+// stages as numbl RuntimeValues.
83+const BASIS_FIELDS = ['f', 'omega', 'Q', 'g', 'B', 'apX', 'apW', 'Sm', 'fTarget'];
84+
85+interface 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+}
94+
95+interface 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+}
106+
107+const basisCache = new Map<string, BasisEntry>();
108+const jetCache = new Map<string, JetEntry>();
109+const BASIS_CACHE_MAX = 6;
110+const JET_CACHE_MAX = 4;
111+
112+function 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+}
120+
121+function 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+}
129+
130+function keyFor(params: RecorderParams, keys: readonly (keyof RecorderParams)[]): string {
131+ return keys.map((k) => String(params[k])).join('|');
132+}
133+
134+// MATLAB source lines assigning every parameter onto a struct P.
135+function 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+}
143+
144+function 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+}
152+
153+function 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+}
158+
159+function 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+}
166+
167+function 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+}
173+
174+const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
175+
176+self.onmessage = (ev: MessageEvent<RenderRequest>) => {
177+ const msg = ev.data;
178+ if (msg.type !== 'render') return;
179+
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+ };
188+
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+};
212+
213+function render(msg: RenderRequest): void {
214+ const { id, params, mfiles } = msg;
215+ let quality = msg.quality;
216+
217+ const codeKey = hashMFiles(mfiles);
218+ const basisParamsKey = codeKey + '|' + keyFor(params, BASIS_KEYS);
219+ const jetParamsKey = keyFor(params, JET_KEYS);
220+
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;
232+
233+ const jetKeyAt = (i: number, os: number) =>
234+ basisParamsKey + `|n=${notes[i]}` + '||' + jetParamsKey + `|span=${spans[i]}|ord=${i + 1}|os=${os}`;
235+
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;
242+
243+ const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
244+
245+ let basisMs = 0;
246+ let jetMs = 0;
247+ let basisCached = true;
248+ let jetCached = true;
249+ const basisEntries: BasisEntry[] = [];
250+ const jetEntries: JetEntry[] = [];
251+
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)}
262+basis = recorder_basis(P, ${notes[i]});
263+${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
264+b_profP = basis.profP;
265+b_profX = basis.profX;
266+b_gridF = basis.gridF;
267+b_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);
284+
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')}
296+basis.n = length(basis.f);
297+[Zr, Zi, Qin, jinfo] = recorder_jet(P, basis, ${spans[i]}, ${i + 1}, ${fsHi});
298+fPlay = jinfo.fPlay;
299+fTarget = jinfo.fTarget;
300+cents = jinfo.cents;
301+qPlot = 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+ }
319+
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')}
332+basis.n = length(basis.f);
333+[pn, radinfo] = recorder_radiate(P, basis, n${i + 1}_Zr, n${i + 1}_Zi, n${i + 1}_Qin, ${fsHi});
334+k0 = min(T, round(${tOn[i]} * ${fsHi}) + 1);
335+len = min(length(pn), T - k0 + 1);
336+p(k0:k0 + len - 1) = p(k0:k0 + len - 1) + pn(1:len);
337+inside = max(inside, radinfo.listenerInside);
338+distM = radinfo.listenerDistM;
339+`);
340+ }
341+ const rOut = runStage(mfiles, `
342+${pLiteral(params)}
343+T = floor(${durationS} * ${fsHi});
344+p = zeros(1, T);
345+inside = 0;
346+distM = 0;
347+${noteBlocks.join('\n')}
348+[y, peakPa] = recorder_output(P, p, ${oversample});
349+`, inputs);
350+ const radiateMs = performance.now() - t0;
351+
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);
356+
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+ }
365+
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+ }
385+
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+}
webapp/tsconfig.jsonadded+20−0View file
@@ -0,0 +1,20 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "useDefineForClassFields": true,
5+ "lib": ["ES2022", "DOM", "DOM.Iterable", "WebWorker"],
6+ "module": "ESNext",
7+ "skipLibCheck": true,
8+ "moduleResolution": "bundler",
9+ "resolveJsonModule": true,
10+ "isolatedModules": true,
11+ "noEmit": true,
12+ "jsx": "react-jsx",
13+ "strict": true,
14+ "noUnusedLocals": true,
15+ "noUnusedParameters": true,
16+ "noFallthroughCasesInSwitch": true
17+ },
18+ "include": ["src"],
19+ "exclude": ["src/selftest.ts"]
20+}
webapp/vite.config.tsadded+9−0View file
@@ -0,0 +1,9 @@
1+import { defineConfig } from 'vite';
2+import react from '@vitejs/plugin-react';
3+
4+export default defineConfig({
5+ plugins: [react()],
6+ base: './',
7+ worker: { format: 'es' },
8+ server: { fs: { allow: ['..'] } }, // the .m sources live in ../matlab
9+});