0e92d35Pizzicato standup bass: coupled string-body modal model, numbl MATLAB reference, and interactive web appJeremy Magland 1function [p, radinfo] = bass_radiate(P, basis, tOn, str, stop, force, fsHi)
2% BASS_RADIATE Pluck projection, closed-form synthesis, and free-field
3% pressure at the listener.
4%
5% The pluck is an initial condition, not a force history. Releasing a
6% string pulled with force F at fraction pluck_pos of the vibrating length
7% from the bridge leaves the classic triangular shape, whose sine
8% coefficients are b_k = 2 F Lv sin(k pi xp/Lv) / (T pi^2 k^2); projected
9% through PM = (M Phi)' from bass_basis.m these become initial amplitudes
10% of the coupled modes, with zero initial velocity. The finite release from
11% the fingertip is a first-order rolloff 1/sqrt(1 + (om tau)^2) applied to
12% the coupled amplitudes (finger_release_ms), which is what separates a
13% fleshy jazz pizzicato from a hard snap. From there every coupled mode is a
14% damped exponential; nothing is integrated. When the next note lands, the
15% previous note is choked: its poles pick up an extra decay 1/tau_choke
16% (choke_ms), standing in for the player stopping the string, and the last
17% note rings through the full tail.
18%
19% Radiation, per coupled mode at its own centre frequency (the modes are
20% narrow-band):
21% - the monopole volume Vm radiates from the two f-holes, half each, as
22% point monopoles with their own retardations;
23% - the dipole volume Dd radiates from the top-plate centre as a point
24% dipole along the plate normal (+z), with the near-field 1/r term kept;
25% - the direct term: above the retained body modes the bridge sees the
26% resistive residual admittance residual_y, so the bridge force (Fcoef
27% per unit modal amplitude) moves the plate broadband, and that motion
28% radiates as a small dipole. This is how the high partials reach the
29% listener at all, and it is the least principled term of the model;
30% direct_gain scales it and 0 turns it off.
31% The common carrier exp(-i k rc) to the bridge is stripped, so relative
32% phase between the f-holes, the plate, and the listener position is kept
33% while the overall propagation delay is not double counted.
34%
35% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.
37c = P.speed_of_sound;
38rho0 = P.air_density;
39K = length(tOn);
40beat = 60 / max(20, P.tempo_bpm);
41duration = (K - 1) * beat + P.tail_s;
42T = floor(duration * fsHi);
43dt = 1 / fsHi;
44tauRel = P.finger_release_ms * 1e-3;
45tauChoke = max(1e-3, P.choke_ms * 1e-3);
47% ---- source and listener geometry (origin at the bridge) ----
48lx = P.listener_x_m;
49ly = P.listener_y_m;
50lz = P.listener_z_m;
51rc = sqrt(lx^2 + ly^2 + lz^2);
52fhX = 0.095; fhY = -0.02; fhZ = 0.11; % the two f-holes
53dpX = 0.0; dpY = -0.10; dpZ = 0.11; % top-plate dipole centre
54r1 = max(0.05, sqrt((lx - fhX)^2 + (ly - fhY)^2 + (lz - fhZ)^2));
55r2 = max(0.05, sqrt((lx + fhX)^2 + (ly - fhY)^2 + (lz - fhZ)^2));
56rd = max(0.05, sqrt((lx - dpX)^2 + (ly - dpY)^2 + (lz - dpZ)^2));
57cosg = (lz - dpZ) / rd; % dipole axis is the plate normal +z
58listenerInside = rc < 0.35;
cfafae8Retilt the body radiativity toward the low registerJeremy Magland 59Sdir = 0.004 * P.rad_area_scale;
0e92d35Pizzicato standup bass: coupled string-body modal model, numbl MATLAB reference, and interactive web appJeremy Magland 60
61% ---- per-fingering radiation transfer, poles, and unit pluck amplitudes ----
62nf = basis.nf;
63ntot = basis.n;
64G = zeros(ntot, 1) + 0i;
65lam = zeros(ntot, 1) + 0i;
66eta0u = zeros(ntot, 1); % per unit pluck force
67for u = 1:nf
68 n = basis.nModes(u);
69 ns = basis.ns(u);
70 Lv = basis.Lv(u);
71 Tt = basis.T(u);
72 mu = basis.mu(u);
73 rows = basis.off(u) + (1:n);
74 om = 2 * pi * basis.f(rows);
75 zt = basis.zeta(rows);
77 ph1 = exp(-1i * (om / c) * (r1 - rc));
78 ph2 = exp(-1i * (om / c) * (r2 - rc));
79 phd = exp(-1i * (om / c) * (rd - rc));
80 near = 1i * (om / c) + 1 / rd;
81 Gm = -(om.^2 * rho0 / (4 * pi)) .* (basis.Vm(rows) / 2) .* (ph1 / r1 + ph2 / r2);
82 Gd = -(om.^2 * rho0 / (4 * pi)) .* basis.Dd(rows) .* near .* phd * (cosg / rd);
83 Gr = (1i * om * rho0 / (4 * pi)) * (Sdir * P.residual_y * P.dipole_length_m) ...
84 .* basis.Fcoef(rows) .* near .* phd * (cosg / rd);
85 G(rows) = Gm + Gd + P.direct_gain * Gr;
87 omd = om .* sqrt(max(1e-12, 1 - zt.^2));
88 lam(rows) = -zt .* om + 1i * omd;
90 % unit-force pluck projection through the coupled shapes
91 xp = Lv * (1 - min(0.95, max(0.05, P.pluck_pos)));
92 kv = (1:ns)';
93 bk = (2 * Lv / (Tt * pi^2)) * sin(kv * pi * xp / Lv) ./ kv.^2;
94 q0 = bk * sqrt(mu * Lv / 2);
95 PM = reshape(basis.PMflat(basis.pmOff(u) + (1:n * ns)), n, ns);
96 e0 = PM * q0;
97 e0 = e0 ./ sqrt(1 + (om * tauRel).^2); % finite finger release, first order
98 eta0u(rows) = e0;
99end
101% ---- per-note synthesis: damped exponentials, choked at the next onset ----
102ks = min(T, max(1, round(tOn * fsHi) + 1));
103p = zeros(1, T);
104CH = 32768;
105chokeSpan = ceil(10 * tauChoke * fsHi);
106peakEta = 0;
107for i = 1:K
108 u = basis.fIdx(i);
109 n = basis.nModes(u);
110 rows = basis.off(u) + (1:n);
111 om = 2 * pi * basis.f(rows);
112 zt = basis.zeta(rows);
113 omd = om .* sqrt(max(1e-12, 1 - zt.^2));
114 cAmp = (force(i) * eta0u(rows)) .* (1 - 1i * (zt .* om ./ omd));
115 s = G(rows) .* cAmp;
116 lamN = lam(rows);
117 peakEta = max(peakEta, force(i) * max(abs(eta0u(rows))));
119 if i < K
120 segEnd = max(ks(i), ks(i + 1));
121 stopAt = min(T + 1, segEnd + chokeSpan);
122 else
123 segEnd = T + 1;
124 stopAt = T + 1;
125 end
127 pos = ks(i);
128 while pos < segEnd
129 len = min(CH, segEnd - pos);
130 seg = real((s.') * exp(lamN * ((0:len - 1) * dt)));
131 p(pos:pos + len - 1) = p(pos:pos + len - 1) + seg;
132 s = s .* exp(lamN * (len * dt));
133 pos = pos + len;
134 end
135 if pos < stopAt
136 lamC = lamN - 1 / tauChoke;
137 while pos < stopAt
138 len = min(CH, stopAt - pos);
139 seg = real((s.') * exp(lamC * ((0:len - 1) * dt)));
140 p(pos:pos + len - 1) = p(pos:pos + len - 1) + seg;
141 s = s .* exp(lamC * (len * dt));
142 pos = pos + len;
143 end
144 end
145end
147% ---- diagnostics: the fundamental of the first note ----
148u1 = basis.fIdx(1);
149rows1 = basis.off(u1) + (1:basis.nModes(u1));
150a1 = abs(eta0u(rows1));
151best = 1;
152bv = 0;
153for m = 1:length(rows1)
154 if a1(m) > bv
155 bv = a1(m);
156 best = m;
157 end
158end
159f1 = basis.f(rows1(best));
160z1 = basis.zeta(rows1(best));
161radinfo.fundHz = f1;
162radinfo.fundT60 = log(1000) / max(1e-9, z1 * 2 * pi * f1);
163radinfo.listenerDistM = rc;
164radinfo.listenerInside = listenerInside;
165radinfo.gainAbs = abs(G);
166radinfo.peakEta = peakEta;
167end