1function basis = bass_basis(P, str, stop)
2% BASS_BASIS Coupled string-body mode set, one eigensolve per fingering.
3%
4% For each fingering of the gesture (a string and a stop, deduplicated so a
5% walking line that revisits a note reuses its solve), the vibrating string
6% and the modal body of bass_body.m are assembled into ONE symmetric
7% generalized eigenproblem and diagonalized together, so the coupling
8% through the bridge is two-way: the string drives the body, and the moving
9% bridge detunes and damps the string in return. This is where wolf-note
10% physics lives; a note whose fundamental lands on a strong body resonance
11% comes out as a split, fast-beating pair of coupled modes rather than one.
12%
13% The displacement ansatz is
14% u(x) = sum_k q_k phi_k(x) + y_b x/Lv, y_b = sum_j a_j xi_j,
15% with phi_k the mass-normalized pinned-pinned modes of the stiff string
16% (om_k^2 = (k pi/Lv)^2 T/mu + (k pi/Lv)^4 EI/mu; the bending stiffness EI
17% comes from the inharmonicity parameter, EI = B T L0^2/pi^2, and is what
18% stretches the high partials sharp) and x/Lv the static shape the string
19% takes when the bridge end moves. With this split the elastic coupling
20% integrals vanish (the ramp has no curvature and phi_k vanishes at both
21% ends), so all the coupling sits in the mass matrix, plus the ramp's own
22% inertia and the tension load T/Lv on the body block; bass_assemble.m
23% builds the matrices. Both are symmetric and M is positive definite, so
24% the problem reduces through the Cholesky factor of M and is diagonalized
25% by svd, which keeps near-degenerate coupled pairs orthonormal where a
26% nonsymmetric eigensolver does not. The shapes come back mass-normalized.
27%
28% Damping is assigned in the uncoupled coordinates and projected,
29% zeta_m = (Phi' C0 Phi)_mm / (2 om_m): string internal loss (loss_factor +
30% loss_slope om), the finger termination on stopped notes (finger_loss), the
31% resistive residual bridge admittance above the retained body modes
32% (residual_y, high-passed so it does not double-count the modal region),
33% and the body Qs. Radiation damping from each coupled mode's own monopole
34% and dipole strengths is added on top, scaled by radiation_loss.
35%
36% Stored per coupled mode: frequency, damping, monopole volume Vm and dipole
37% volume Dd per unit modal coordinate (what bass_radiate.m radiates), the
38% bridge force coefficient Fcoef (what drives the direct term), and the
39% projection matrix PM = (M Phi)' restricted to the string block, which
40% turns a plucked initial shape into coupled initial amplitudes. Fields are
41% concatenated across fingerings with offset indexing (numbl has no struct
42% arrays).
44rho0 = P.air_density;
45cair = P.speed_of_sound;
47% ---- deduplicate fingerings ----
48K = length(str);
49fkey = str * 100 + stop;
50uKey = zeros(1, K);
51fIdx = zeros(1, K);
52nf = 0;
53for i = 1:K
54 hit = 0;
55 for u = 1:nf
56 if uKey(u) == fkey(i)
57 hit = u;
58 end
59 end
60 if hit == 0
61 nf = nf + 1;
62 uKey(nf) = fkey(i);
63 hit = nf;
64 end
65 fIdx(i) = hit;
66end
68% ---- per-fingering eigensolve, concatenated storage ----
69f = [];
70zeta = [];
71Vm = [];
72Dd = [];
73Fcoef = [];
74ybAmp = [];
75bodyFrac = [];
76PMflat = [];
77off = zeros(nf, 1);
78nModes = zeros(nf, 1);
79pmOff = zeros(nf, 1);
80nsArr = zeros(nf, 1);
81LvArr = zeros(nf, 1);
82muArr = zeros(nf, 1);
83TArr = zeros(nf, 1);
84f0Arr = zeros(nf, 1);
85strArr = zeros(nf, 1);
86stopArr = zeros(nf, 1);
87body = bass_body(P);
89for u = 1:nf
90 s = floor(uKey(u) / 100);
91 st = uKey(u) - 100 * s;
92 [M, Kmat, omS, zS, meta] = bass_assemble(P, s, st);
93 ns = meta.ns;
94 nb = meta.nb;
95 n = meta.n;
96 Lv = meta.Lv;
97 mu = meta.mu;
98 T = meta.T;
99 kv = (1:ns)';
101 Rc = chol(M);
102 Am = (Rc') \ (Kmat / Rc);
103 Am = (Am + Am') / 2;
104 [Uv, Sv, ~] = svd(Am);
105 omc2 = flipud(diag(Sv)); % ascending
106 Vv = fliplr(Uv);
107 omc = sqrt(max(1e-6, omc2));
108 Phi = Rc \ Vv; % mass-normalized coupled shapes, n x n
110 % projected damping plus per-mode radiation damping
111 cvec = 2 * [zS .* omS; body.zeta .* body.omega];
112 zc = ((cvec') * (Phi.^2))' ./ (2 * omc);
113 PhiB = Phi(ns + 1:n, :);
114 bodyDisp = (body.a * ones(1, n)) .* PhiB; % nb x n physical
115 vm = ((body.S .* body.mono)' * bodyDisp)'; % n x 1 monopole volume
116 dd = P.dipole_length_m * ((body.S .* (1 - body.mono))' * bodyDisp)';
117 zRad = P.radiation_loss * (rho0 * omc.^4 .* vm.^2 / (8 * pi * cair) ...
118 + rho0 * omc.^6 .* dd.^2 / (24 * pi * cair^3)) ./ omc.^3;
119 zc = zc + zRad;
121 yb = ((body.a') * PhiB)'; % n x 1 bridge displacement
122 phiSlope = sqrt(2 / (mu * Lv)) * (kv * pi / Lv) .* ((-1).^kv);
123 fc = T * (((phiSlope') * Phi(1:ns, :))' + yb / Lv);
124 bfrac = ((ones(1, nb) * (PhiB.^2)))';
126 PM = (M * Phi);
127 PM = (PM(1:ns, :))'; % n x ns projection
129 off(u) = length(f);
130 nModes(u) = n;
131 pmOff(u) = length(PMflat);
132 nsArr(u) = ns;
133 LvArr(u) = Lv;
134 muArr(u) = mu;
135 TArr(u) = T;
136 f0Arr(u) = meta.f0;
137 strArr(u) = s;
138 stopArr(u) = st;
140 f = [f; omc / (2 * pi)];
141 zeta = [zeta; zc];
142 Vm = [Vm; vm];
143 Dd = [Dd; dd];
144 Fcoef = [Fcoef; fc];
145 ybAmp = [ybAmp; yb];
146 bodyFrac = [bodyFrac; bfrac];
147 PMflat = [PMflat; reshape(PM, n * ns, 1)];
148end
150basis.nf = nf;
151basis.fIdx = fIdx;
152basis.off = off;
153basis.nModes = nModes;
154basis.pmOff = pmOff;
155basis.ns = nsArr;
156basis.Lv = LvArr;
157basis.mu = muArr;
158basis.T = TArr;
159basis.f0 = f0Arr;
160basis.str = strArr;
161basis.stop = stopArr;
162basis.f = f;
163basis.zeta = zeta;
164basis.Vm = Vm;
165basis.Dd = Dd;
166basis.Fcoef = Fcoef;
167basis.ybAmp = ybAmp;
168basis.bodyFrac = bodyFrac;
169basis.PMflat = PMflat;
170basis.n = length(f);
171basis.bodyF = body.f;
172basis.bodyZeta = body.zeta;
173basis.bodyA = body.a;
174basis.bodyS = body.S;
175basis.bodyMono = body.mono;
176end