5623b11Recorder simulator: MATLAB jet-drive model and interactive web appJeremy Magland 1function [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).
31rho0 = P.air_density;
32c = P.speed_of_sound;
33nT = size(Zr, 2);
34n = basis.n;
36lx = P.listener_x_m;
37ly = P.listener_y_m;
38lz = P.listener_z_m;
40nap = length(basis.apX);
41rap = zeros(nap, 1);
42for a = 1:nap
43 rap(a) = sqrt(lx^2 + (ly + basis.apX(a))^2 + lz^2);
44end
45listenerInside = min(rap) < 0.02;
46rap = max(rap, 0.01);
47rcen = rap(1); % the window
49% modal analytic-signal weights (must mirror recorder_jet.m)
50dScale = rho0 * (4 / pi) * sqrt(2 * P.flue_h_m * P.cutup_m) / basis.Sm;
51cn = dScale * basis.B;
52omd = basis.omega .* sqrt(max(1e-12, 1 - 1 ./ (4 * basis.Q.^2)));
53lam = -basis.g / 2 + 1i * omd;
54rho = (basis.g .* lam + basis.omega.^2) ./ (2i * omd);
56H = zeros(1, n) + 0i;
57for 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));
62end
64p = real(H) * Zr - imag(H) * Zi; % 1 x nT
66% direct broadband monopole from the window
67dq = zeros(1, nT);
68if nT > 2
69 dq(2:nT - 1) = (Qin(3:nT) - Qin(1:nT - 2)) * (fsHi / 2);
70end
71p = p + P.direct_gain * (1 - sum(cn)) * rho0 / (4 * pi * rcen) * dq;
73radinfo.gainAbs = abs(H);
74radinfo.listenerInside = listenerInside;
75radinfo.listenerDistM = rcen;
76end