1function [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.
53rho0 = P.air_density;
54dt = 1 / fsHi;
55nT = max(16, floor(spanS * fsHi));
56n = basis.n;
58W = P.cutup_m;
59h = P.flue_h_m;
60Sm = basis.Sm;
61bj = 0.41 * (h^2 * W)^(1 / 3); % Bickley jet half-width at the labium
62growth = exp(min(2.5, P.jet_growth * W / h));
63dScale = rho0 * (4 / pi) * sqrt(2 * h * W) / Sm;
64lossC = 0.5 * rho0 * P.vortex_loss / Sm; % series resistance per unit |v|
65y0 = P.labium_offset_m;
66jg = P.jet_gain;
68% ---- breath pressure envelope and derived vectors (all precomputed) ----
69t = (0:nT - 1) * dt;
70tauA = max(1e-4, P.attack_ms * 1e-3) / 2.3;
71tauR = max(1e-4, P.release_ms * 1e-3) / 2.3;
72env = 1 - exp(-t / tauA);
73rel = find(t > P.note_dur_s);
74env(rel) = env(rel) .* exp(-(t(rel) - P.note_dur_s) / tauR);
75vib = 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.
80trackFac = (basis.fTarget / 523.25)^P.breath_track;
81pf = P.breath_pa * trackFac * env .* vib;
82Uj = sqrt(2 * max(pf, 0) / rho0);
83UjE = max(Uj, 1.0); % guarded: receptivity and cutoff need Uj > 0
84tauS = min(0.005, W ./ (P.jet_conv * UjE)) * fsHi;
85tauS = max(2, tauS);
86aJet = 1 - exp(-2 * pi * (P.jet_stcut * UjE / W) * dt); % jet instability cutoff
87etaAmp = P.noise_level * h * (UjE / 8).^1.5;
89% ---- seeded turbulence noise, filtered in the loop ----
90state = mod(floor(abs(P.seed)) + 7919 * noteOrd, 4294967296);
91[u, state] = recorder_rng(state, nT);
92white = 2 * u.' - 1;
93aN1 = 1 - exp(-2 * pi * 4000 * dt); % turbulence band: 400 Hz .. 4 kHz
94aN2 = 1 - exp(-2 * pi * 400 * dt);
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.
99omd = basis.omega .* sqrt(max(1e-12, 1 - 1 ./ (4 * basis.Q.^2)));
100lam = -basis.g / 2 + 1i * omd;
101Az = exp(lam * dt);
102ar = real(Az);
103ai = imag(Az);
104cn = dScale * basis.B;
105rho = (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))
107wr = -2 * cn .* real(rho);
108wi = 2 * cn .* imag(rho);
109cSum = sum(cn);
111% ---- the loop ----
112Zr = zeros(n, nT);
113Zi = zeros(n, nT);
114Qin = zeros(1, nT);
115vh = zeros(1, nT);
116zr = zeros(n, 1);
117zi = zeros(n, 1);
118vprev = 0;
119vlp = 0;
120vlp2 = 0;
121nlp1 = 0;
122nlp2 = 0;
123for 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;
150end
152% ---- diagnostics: playing frequency from the sustained half ----
153k0 = max(1, floor(0.35 * P.note_dur_s * fsHi));
154k1 = min(nT, floor(0.95 * P.note_dur_s * fsHi));
155seg = vh(k0:k1);
156mu = sum(seg) / length(seg);
157seg = seg - mu;
158nzc = 0;
159kFirst = 0;
160kLast = 0;
161for 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
169end
170if nzc > 2
171 jinfo.fPlay = (nzc - 1) * fsHi / (kLast - kFirst);
172else
173 jinfo.fPlay = 0;
174end
175jinfo.fTarget = basis.fTarget;
176if jinfo.fPlay > 0
177 jinfo.cents = 1200 * log(jinfo.fPlay / basis.fTarget) / log(2);
178else
179 jinfo.cents = 0;
180end
181jinfo.vacRms = sqrt(sum(seg.^2) / length(seg));
182jinfo.ujMax = max(Uj);
183jinfo.qinMean = sum(Qin) / nT;
184end