/ simulated_instruments / bass
Sign in
simulated_instruments / bass
d71e860
Branches
main
bass / matlab / bass_validate.m
197 lines · 7.1 KBCodeBlameHistory
0e92d35Pizzicato standup bass: coupled string-body modal model, numbl MATLAB reference, and interactive web appJeremy Magland 1% Validate the bass model. There is no external reference implementation,
2% so the checks are (a) closed-form identities the physics must satisfy,
3% (b) structural identities of the coupled eigensolve, (c) exact
4% reproducibility of the seeded walking line, and (d) pinned end-to-end
5% diagnostics from the default configuration, recorded when the model was
6% built and treated as a regression reference since.
7% Run with: numbl run bass_validate.m
9fails = 0;
11% ---- 1. decoupled stiff string against the closed form ----
12% With the bridge mobility driven to zero the body cannot move, and the
13% coupled solve must hand back the textbook pinned-pinned stiff string,
14% f_k = sqrt((k pi/Lv)^2 T/mu + (k pi/Lv)^4 EI/mu) / 2 pi, including the
15% inharmonic stretch. Checked on the open A with an exaggerated B.
16P = bass_defaults();
17Q = P;
18Q.bridge_mobility = 1e-9;
19Q.inharmonicity = 2e-4;
20basis = bass_basis(Q, 3, 0);
21[~, ~, omS, ~, meta] = bass_assemble(Q, 3, 0);
22errS = 0;
23for k = 1:10
24 fClosed = omS(k) / (2 * pi);
25 fGot = bv_nearest(basis.f, fClosed);
26 errS = max(errS, abs(fGot - fClosed) / fClosed);
27end
28stretch = (omS(10) / (10 * omS(1))) - 1;
29ok = errS < 1e-6 && stretch > 0.005;
30fails = fails + (~ok);
31fprintf('stiff string: err %.1e vs closed form, partial-10 stretch %.3f %s\n', ...
32 errS, stretch, bv_verdict(ok));
34% ---- 2. eigensolve identities on the full coupled system ----
35% The Cholesky reduction must hand back mass-orthonormal shapes that
36% diagonalize the stiffness: Phi' M Phi = I and Phi' K Phi = diag(om^2).
37% Checked on the open E, the largest system of the default chart.
38[M, Kmat, ~, ~, meta] = bass_assemble(P, 4, 0);
39Rc = chol(M);
40Am = (Rc') \ (Kmat / Rc);
41Am = (Am + Am') / 2;
42[Uv, Sv, ~] = svd(Am);
43om2 = flipud(diag(Sv));
44Phi = Rc \ fliplr(Uv);
45eOrth = max(max(abs(Phi' * M * Phi - eye(meta.n))));
46KP = Phi' * Kmat * Phi;
47eDiag = max(max(abs(KP - diag(om2)))) / max(om2);
48ok = eOrth < 1e-8 && eDiag < 1e-10;
49fails = fails + (~ok);
50fprintf('eigensolve: orthonormality %.1e, diagonalization %.1e (n = %d) %s\n', ...
51 eOrth, eDiag, meta.n, bv_verdict(ok));
53% ---- 3. a midpoint pluck silences the even partials ----
54% b_k ~ sin(k pi xp/Lv): plucking at the middle of the vibrating length
55% puts every even string mode at a node of the initial shape, so the
56% second partial must all but vanish from the rendered pressure, while an
57% ordinary pluck position leaves it strong. This exercises the projection,
58% the synthesis, and the radiation end to end.
59Q = P;
60Q.normalize = 0;
61Q.tail_s = 3.0;
62Q.pluck_pos = 0.5;
63[y1, b1, d1] = bass_render(Q);
64fs = Q.sample_rate;
65f1 = d1.fundHz;
66f2 = bv_nearest(b1.f, 2 * f1);
67relMid = 20 * log10(bv_level(y1, fs, f2) / bv_level(y1, fs, f1) + 1e-15);
68Q.pluck_pos = 0.30;
69[y2, b2, d2] = bass_render(Q);
70f2b = bv_nearest(b2.f, 2 * d2.fundHz);
71relOrd = 20 * log10(bv_level(y2, fs, f2b) / bv_level(y2, fs, d2.fundHz) + 1e-15);
72ok = relMid < -35 && relOrd > -25;
73fails = fails + (~ok);
74fprintf('midpoint pluck: partial 2 at %.1f dB (ordinary pluck %.1f dB) %s\n', ...
75 relMid, relOrd, bv_verdict(ok));
77% ---- 4. two-way coupling: veering at a string-body coincidence ----
78% Tuning the strongest corpus mode onto the open G fundamental must split
79% it into a mixed pair (this is the wolf-note mechanism), and the splitting
80% must grow as the square root of the bridge mobility, which is what a
81% two-way coupled pair does and a one-way filter does not.
82d05 = bv_split(P, 0.5);
83d20 = bv_split(P, 2.0);
84ratio = d20.split / d05.split;
85ok = abs(ratio - 2) < 0.25 && d20.mix > 0.15;
86fails = fails + (~ok);
87fprintf('veering: split %.2f -> %.2f Hz for 4x mobility (ratio %.2f, mix %.2f) %s\n', ...
88 d05.split, d20.split, ratio, d20.mix, bv_verdict(ok));
90% ---- 5. far-field falloff: -6 dB per doubling of distance ----
91% RMS over the render, with the two listeners on the same ray from the
92% bridge so the dipole directivity is held fixed. At the fundamental the
93% dipole near-field term still contributes a few hundredths of a dB at
94% 4 m, which is physics rather than error; the tolerance covers it.
95Q = P;
96Q.normalize = 0;
97Q.tail_s = 2.5;
98Q.listener_x_m = 1.2; Q.listener_y_m = 0.75; Q.listener_z_m = 3.0;
99[ya, ~, ~] = bass_render(Q);
100Q.listener_x_m = 2.4; Q.listener_y_m = 1.5; Q.listener_z_m = 6.0;
101[yb, ~, ~] = bass_render(Q);
102drop = 10 * log10(sum(yb.^2) / sum(ya.^2));
103ok = abs(drop + 6.02) < 0.35;
104fails = fails + (~ok);
105fprintf('falloff 3.3 m -> 6.6 m: %.2f dB rms (target -6.02) %s\n', drop, bv_verdict(ok));
107% ---- 6. the seeded walking line is exactly reproducible ----
108Q = P;
109Q.note_count = 8;
110[t1, m1, s1, st1, fo1, g1] = bass_notes(Q);
111[t2, m2, s2, st2, fo2, g2] = bass_notes(Q);
112same = max(abs(t1 - t2)) == 0 && max(abs(m1 - m2)) == 0 ...
113 && max(abs(s1 - s2)) == 0 && max(abs(fo1 - fo2)) == 0;
114ok = same ...
115 && g1.count == 8 ...
116 && sum(m1) == 293 ...
117 && abs(g1.meanForce - 14.014948614) < 1e-6 ...
118 && abs(t1(2) - 0.623242473696) < 1e-9;
119fails = fails + (~ok);
120fprintf('walk (seed %d): %d notes, sum midi %d, mean force %.6f N %s\n', ...
121 Q.seed, g1.count, sum(m1), g1.meanForce, bv_verdict(ok));
123% ---- 7. pinned end-to-end diagnostics at the defaults ----
124[y, basis, dinfo] = bass_render(P);
cfafae8Retilt the body radiativity toward the low registerJeremy Magland 125ok = abs(dinfo.peakPa - 0.270351) < 0.001 ...
126 && abs(dinfo.fundHz - 54.4515) < 0.005 ...
128fails = fails + (~ok);
129fprintf('pinned render: peak %.6f Pa (target %.6f), fund %.4f Hz (target %.4f), %d modes %s\n', ...
cfafae8Retilt the body radiativity toward the low registerJeremy Magland 130 dinfo.peakPa, 0.270351, dinfo.fundHz, 54.4515, basis.n, bv_verdict(ok));
132% ---- 8. the note decays ----
133% An open-string pizzicato rings for seconds, so the tail is not silent;
134% but the energy in the last half second must sit well below the opening
135% half second, which carries the attack.
136fs = P.sample_rate;
137n = length(y);
138head = y(1:round(0.5 * fs));
139tail = y(n - round(0.5 * fs):n);
140dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
141ok = dropE < -15;
142fails = fails + (~ok);
143fprintf('note decay: last/first half second %.1f dB %s\n', dropE, bv_verdict(ok));
145if fails > 0
146 fprintf('%d CHECK(S) FAILED\n', fails);
147else
148 fprintf('all checks passed\n');
149end
151function fv = bv_nearest(fvec, target)
152% The entry of fvec nearest target.
153fv = fvec(1);
154best = abs(fvec(1) - target);
155for i = 2:length(fvec)
156 d = abs(fvec(i) - target);
157 if d < best
158 best = d;
159 fv = fvec(i);
160 end
161end
162end
164function a = bv_level(y, fs, f)
165% Amplitude of y at frequency f over the first two seconds (plain
166% quadrature dot products; the window is long against 1/f).
167n = min(length(y), floor(2.0 * fs));
168t = (0:n - 1) / fs;
169seg = y(1:n);
170a = sqrt((seg * (cos(2 * pi * f * t))')^2 + (seg * (sin(2 * pi * f * t))')^2);
171end
173function r = bv_split(P, mobility)
174% Frequency splitting of the coupled pair when the strongest corpus mode
175% is tuned onto the open G fundamental, plus how mixed the pair is.
176Q = P;
177Q.bridge_mobility = mobility;
178Q.body_trim = 97.9989 / 97.0;
179basis = bass_basis(Q, 1, 0);
180f0 = 97.9989;
181% the two coupled modes nearest the coincidence
182d = abs(basis.f - f0);
183[~, order] = sort(d);
184i1 = order(1);
185i2 = order(2);
186r.split = abs(basis.f(i1) - basis.f(i2));
187r.mix = min(min(basis.bodyFrac(i1), 1 - basis.bodyFrac(i1)), ...
188 min(basis.bodyFrac(i2), 1 - basis.bodyFrac(i2)));
189end
191function s = bv_verdict(ok)
192if ok
193 s = 'ok';
194else
195 s = 'MISMATCH';
196end
197end
↑↓ move↵openescclose