5623b11Recorder simulator: MATLAB jet-drive model and interactive web appJeremy Magland 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
10fails = 0;
11P = recorder_defaults();
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.
19Q = P;
20Q.bore_r_head_m = 0.006;
21Q.bore_r_foot_m = 0.006;
22geo = recorder_geometry(Q);
23geo.holeRad = 1e-9 * ones(1, 8);
24f = 100:1:2500;
25out = recorder_bore(Q, geo, zeros(1, 8), f);
26mag = abs(1 ./ out.Zbore);
27pk = [];
28for i = 2:length(f) - 1
29 if mag(i) > mag(i - 1) && mag(i) >= mag(i + 1)
30 pk = [pk, f(i)];
31 end
32end
33c = Q.speed_of_sound;
34Leff = geo.boreLen + 0.61 * 0.006;
35fClosed = (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
38dv = sqrt(2 * 1.5e-5 ./ (2 * pi * fClosed));
39fClosed = fClosed ./ (1 + (dv / 0.012) * (1 + 0.4 / sqrt(0.71)));
40errC = max(abs(pk(1:3) - fClosed) ./ fClosed);
41ok = errC < 3e-3;
42fails = fails + (~ok);
43fprintf('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));
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.
52basis = recorder_basis(P, 5);
53geo = recorder_geometry(P);
54bo = recorder_bore(P, geo, geo.fingerings(5, :), reshape(basis.f, 1, basis.n));
55eFit = 0;
56for 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)));
63end
64ok = eFit < 0.06;
65fails = fails + (~ok);
66fprintf('modal fit: worst |Y| error over the first four peaks %.1e %s\n', eFit, rv_verdict(ok));
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.
73Q = P;
74Q.note_dur_s = 0.6;
75Q.tail_s = 0.2;
76fsHi = P.sample_rate * P.oversample;
77centsAbs = zeros(1, 3);
78notes3 = [1 5 11];
79for 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);
83end
84ok = max(centsAbs) < 60;
85fails = fails + (~ok);
86fprintf('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));
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.
94Q = P;
95Q.note_dur_s = 0.6;
96Q.tail_s = 0.2;
97Q.breath_pa = 4 * P.breath_pa;
98b = recorder_basis(Q, 5);
99[~, ~, ~, j] = recorder_jet(Q, b, 0.7, 1, fsHi);
100ratio = j.fPlay / b.f(1);
101ok = ratio > 1.9 && ratio < 2.1;
102fails = fails + (~ok);
103fprintf('overblow at 4x breath: %.1f Hz = %.3f x f1 %s\n', j.fPlay, ratio, rv_verdict(ok));
105% ---- 5. the seeded turbulence is exactly reproducible ----
106Q = P;
107Q.note_dur_s = 0.4;
108Q.tail_s = 0.15;
109[y1, ~, ~] = recorder_render(Q);
110[y2, ~, ~] = recorder_render(Q);
111Q.seed = 2;
112[y3, ~, ~] = recorder_render(Q);
113same = max(abs(y1 - y2)) == 0;
114diff = max(abs(y1 - y3)) > 0;
115ok = same && diff;
116fails = fails + (~ok);
117fprintf('seed: identical rerun %d, different seed differs %d %s\n', same, diff, rv_verdict(ok));
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.
122Q = P;
123Q.normalize = 0;
124Q.note_dur_s = 0.5;
125Q.tail_s = 0.15;
126Q.listener_x_m = 0.5; Q.listener_y_m = 0.4; Q.listener_z_m = 1.0;
127[y1, ~, ~] = recorder_render(Q);
128Q.listener_x_m = 1.0; Q.listener_y_m = 0.8; Q.listener_z_m = 2.0;
129[y2, ~, ~] = recorder_render(Q);
130drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
131ok = abs(drop + 6.02) < 0.35;
132fails = fails + (~ok);
133fprintf('falloff 1.2 m -> 2.4 m: %.2f dB rms (target -6.02) %s\n', drop, rv_verdict(ok));
135% ---- 7. pinned end-to-end diagnostics at the defaults ----
136[y, basis, dinfo] = recorder_render(P);
137ok = 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;
141fails = fails + (~ok);
142fprintf('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));
145% ---- 8. the note stops when the breath stops ----
146fs = P.sample_rate;
147n = length(y);
148sus = y(round(0.5 * fs):round(1.0 * fs));
149tail = y(n - round(0.2 * fs):n);
150dropE = 10 * log10(sum(tail.^2) / length(tail) / (sum(sus.^2) / length(sus)));
151ok = dropE < -40;
152fails = fails + (~ok);
153fprintf('release: last 200 ms %.1f dB re sustain %s\n', dropE, rv_verdict(ok));
155if fails > 0
156 fprintf('%d CHECK(S) FAILED\n', fails);
157else
158 fprintf('all checks passed\n');
159end
161function s = rv_verdict(ok)
162if ok
163 s = 'ok';
164else
165 s = 'MISMATCH';
166end
167end