1% Validate the cymbal model. There is no external reference implementation, so
2% the checks are (a) closed-form identities the physics must satisfy, (b)
3% structural identities of the shell eigensolve and of the radiation integral,
4% (c) exact reproducibility of the seeded strike process, and (d) pinned
5% end-to-end diagnostics from the default configuration, recorded when the
6% model was built and treated as a regression reference since.
7% Run with: numbl run cymbal_validate.m
9fails = 0;
11% ---- 1. flat free circular plate against the closed form ----
12% The same shell routine that builds the cymbal, run on a flat meridian of
13% constant thickness with a pinhole at the centre, must reproduce the
14% completely free circular plate: the roots of the 2 x 2 determinant built
15% from J_n and I_n under vanishing radial moment and Kelvin-Kirchhoff shear,
16% and exactly two rigid-body modes at n = 0 (axial translation and rotation
17% about the axis), two at n = 1 (translation and tilt) and none above.
18P = cymbal_defaults();
19Q = P;
20Q.bow_height_m = 0; Q.bell_height_m = 0;
21Q.thickness_bell_m = 0.001; Q.thickness_rim_m = 0.001;
22Q.hole_radius_m = 0.002; Q.radius_m = 0.2; Q.n_elements = 40;
23Gf = cymbal_profile(Q);
24Dp = Q.youngs_modulus * 0.001^3 / (12 * (1 - Q.poisson_ratio^2));
25scal = sqrt(Dp / (Q.density * 0.001)) / Q.radius_m^2 / (2 * pi);
26Lr = Gf.S / (Gf.nn - 1);
27errP = 0;
28rigid = zeros(1, 5);
29for n = 0:4
30 f = cv_freqs(Gf, n, Q, Lr);
31 rigid(n + 1) = sum(f < 3);
32 fe = f(f > 3);
33 lam = cv_plate_roots(n, Q.poisson_ratio);
34 fc = lam.^2 * scal;
35 nc = min(4, min(length(fe), length(fc)));
36 errP = max(errP, max(abs(fe(1:nc).' - fc(1:nc)) ./ fc(1:nc)));
37end
38ok = errP < 3e-3 && rigid(1) == 2 && rigid(2) == 2 && sum(rigid(3:5)) == 0;
39fails = fails + (~ok);
40fprintf('flat plate : max freq err %.1e, rigid modes %d %d %d %d %d %s\n', ...
41 errP, rigid(1), rigid(2), rigid(3), rigid(4), rigid(5), cv_verdict(ok));
43% ---- 2. eigensolve identities on the cymbal itself ----
44% The Cholesky reduction must hand back mass-orthonormal shapes that
45% diagonalize the stiffness: phi' M phi = I and phi' K phi = diag(om^2).
46G = cymbal_profile(P);
47Lref = G.S / (G.nn - 1);
48[K, M] = cymbal_shell(G, 5, P.youngs_modulus, P.poisson_ratio, P.density, Lref);
49Rc = chol(M);
50Am = (Rc') \ (K / Rc);
51Am = (Am + Am') / 2;
52[Uv, Sv, ~] = svd(Am);
53om2 = flipud(diag(Sv));
54Phi = Rc \ fliplr(Uv);
55eOrth = max(max(abs(Phi' * M * Phi - eye(length(om2)))));
56eDiag = max(max(abs(Phi' * K * Phi - diag(om2)))) / max(om2);
57ok = eOrth < 1e-8 && eDiag < 1e-10;
58fails = fails + (~ok);
59fprintf('eigensolve : orthonormality %.1e, diagonalization %.1e %s\n', eOrth, eDiag, cv_verdict(ok));
61% ---- 3. the Bessel far field against a direct sum over the dipole sheet ----
62% cymbal_basis reduces the azimuthal integral of the far-field kernel to
63% Bessel functions in order to get the radiated power; that reduction must
64% agree with brute force summation over a finely discretized dipole sheet.
65Q = P;
66Q.n_diameters = 4; Q.max_mode_hz = 1200;
67b4 = cymbal_basis(Q);
68errF = 0;
69for m = [1 3 5 9]
70 [pd, pb] = cv_farfield(Q, b4, m, 200, 720, 0.6, 0.4);
71 errF = max(errF, abs(pb / pd - 1));
72end
73ok = errF < 0.01;
74fails = fails + (~ok);
75fprintf('far field : Bessel vs direct dipole sum, max err %.2e %s\n', errF, cv_verdict(ok));
77% ---- 4. on the axis, every mode with a nodal diameter cancels exactly ----
78% A cos(n theta) source on a surface of revolution has no monopole, dipole or
79% any lower multipole than 2^n, and on the symmetry axis all of them vanish.
80% This is the sharpest available check that the azimuthal sum in the listener
81% transfer is assembled correctly.
82Q = P;
83Q.n_diameters = 6; Q.max_mode_hz = 1500;
84Q.listener_x_m = 0; Q.listener_y_m = 0; Q.listener_z_m = 0.6;
85b6 = cymbal_basis(Q);
86[tstr, ~, jstr, taustr, fstr, sstr, sinf] = cymbal_strikes(Q, b6);
87[env, tfr, DV, kf] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.05);
88[~, ri] = cymbal_radiate(Q, b6, env, tfr, DV, kf, tstr, taustr, fstr, sstr, Q.sample_rate);
89g0 = max(ri.gainAbs(b6.ndia == 0));
90g1 = max(ri.gainAbs(b6.ndia >= 1));
91ok = g1 / g0 < 1e-10;
92fails = fails + (~ok);
93fprintf('on axis : |G| of n>=1 modes is %.1e of the n=0 modes %s\n', g1 / g0, cv_verdict(ok));
95% ---- 5. Hertz contact-time scaling tau ~ v^(-1/5) ----
96% With a stick much lighter than the plate mass engaged during the contact the
97% reduced mass is constant and the free-sphere scaling must emerge.
98Q = P;
99Q.strike_count = 1;
100Q.stick_mass_kg = 2e-5;
101[~, ~, ~, tau1] = cymbal_strikes(Q, b6);
102Q.strike_velocity = 2 * P.strike_velocity;
103[~, ~, ~, tau2] = cymbal_strikes(Q, b6);
104ratio = tau2(1) / tau1(1);
105ok = abs(ratio - 2^(-1 / 5)) < 0.03;
106fails = fails + (~ok);
107fprintf('hertz : tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), cv_verdict(ok));
109% ---- 6. the cascade conserves energy, and switches off cleanly ----
110% With the damping removed and the top of the ladder reflecting, the transfer
111% scheme is a Markov generator on the modal energies and must conserve their
112% sum to round-off. With cascade_gain = 0 every envelope must instead collapse
113% onto the closed-form exp(-zeta om t) of an uncoupled damped oscillator.
114Q = P;
115Q.n_diameters = 6; Q.max_mode_hz = 1500; Q.cascade_leak = 0;
116bz = b6;
117bz.zeta = zeros(bz.n, 1);
118[tstr, ~, jstr, taustr, ~, sstr] = cymbal_strikes(Q, bz);
119[envC, ~, DVc] = cymbal_cascade(Q, bz, tstr, jstr, taustr, sstr, 0.4);
120Etot = sum(envC.^2 / 2, 1);
121eCons = max(abs(Etot / Etot(1) - 1));
122Q.cascade_gain = 0;
123[envL, tfrL] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.4);
124pred = envL(:, 1) * exp(-(b6.zeta .* b6.omega).' * 0 + 0);
125pred = (envL(:, 1) * ones(1, length(tfrL))) .* exp(-(b6.zeta .* b6.omega) * tfrL);
126sel = envL(:, 1) > 1e-12;
127eLin = max(max(abs(envL(sel, :) ./ pred(sel, :) - 1)));
128ok = eCons < 1e-10 && eLin < 1e-9;
129fails = fails + (~ok);
130fprintf('cascade : energy drift %.1e, linear limit err %.1e %s\n', eCons, eLin, cv_verdict(ok));
132% ---- 7. far-field falloff: -6 dB per doubling of distance ----
133% RMS over the render, with the two listeners on the same ray from the origin
134% so every element dipole keeps its directivity factor.
135Q = P;
136Q.normalize = 0; Q.n_diameters = 12; Q.max_mode_hz = 3000;
137Q.tail_s = 0.4; Q.oversample = 1;
138Q.listener_x_m = 0.8; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;
139y1 = cymbal_render(Q);
140Q.listener_x_m = 1.6; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;
141y2 = cymbal_render(Q);
142drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
143ok = abs(drop + 6.02) < 0.25;
144fails = fails + (~ok);
145fprintf('falloff : 2 m -> 4 m %.2f dB rms (target -6.02) %s\n', drop, cv_verdict(ok));
147% ---- 8. the seeded ride figure is exactly reproducible ----
148Q = P;
149Q.strike_count = 8;
150[t1, v1, j1, ~, ~, s1, g1] = cymbal_strikes(Q, b6);
151[t2, v2, ~, ~, ~, s2] = cymbal_strikes(Q, b6);
152same = max(abs(t1 - t2)) == 0 && max(abs(v1 - v2)) == 0 && max(abs(s1 - s2)) == 0;
153ok = same ...
154 && g1.count == 8 ...
155 && abs(g1.meanV - 3.111518886) < 1e-6 ...
156 && abs(sum(j1) - 0.150818144) < 1e-9 ...
157 && abs(t1(2) - 0.159427523) < 1e-9;
158fails = fails + (~ok);
159fprintf('ride (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
160 Q.seed, g1.count, g1.meanV, sum(j1), cv_verdict(ok));
162% ---- 9. pinned end-to-end diagnostics at the defaults ----
163[y, basis, dinfo] = cymbal_render(P);
164ok = abs(dinfo.peakPa - 13.4297) < 0.02 ...
165 && abs(basis.f(1) - 49.0808) < 0.02 ...
166 && basis.n == 303;
167fails = fails + (~ok);
168fprintf('pinned : peak %.4f Pa (target %.4f), mode1 %.3f Hz (target %.3f), %d modes %s\n', ...
169 dinfo.peakPa, 13.4297, basis.f(1), 49.0808, basis.n, cv_verdict(ok));
171% ---- 10. the wash decays ----
172% A cymbal rings for seconds, so the tail is not silent, but the energy in the
173% last half second must sit well below that of the opening half second.
174fs = P.sample_rate;
175n = length(y);
176head = y(1:round(0.5 * fs));
177tail = y(n - round(0.5 * fs):n);
178dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
179ok = dropE < -20 && dropE > -80;
180fails = fails + (~ok);
181fprintf('decay : last/first half second %.1f dB %s\n', dropE, cv_verdict(ok));
183if fails > 0
184 fprintf('%d CHECK(S) FAILED\n', fails);
185else
186 fprintf('all checks passed\n');
187end
189function f = cv_freqs(G, n, P, Lref)
190[K, M] = cymbal_shell(G, n, P.youngs_modulus, P.poisson_ratio, P.density, Lref);
191Rc = chol(M);
192Am = (Rc') \ (K / Rc);
193Am = (Am + Am') / 2;
194[~, Sv, ~] = svd(Am);
195f = sort(sqrt(diag(Sv))) / (2 * pi);
196end
198function lam = cv_plate_roots(n, nu)
199% Roots of the free-edge frequency equation of a circular plate, bracketed on
200% a coarse scan and bisected.
201lam = [];
202prev = cv_plate_det(n, nu, 0.30);
203for x = 0.30:0.01:26
204 cur = cv_plate_det(n, nu, x + 0.01);
205 if prev * cur < 0
206 lo = x;
207 hi = x + 0.01;
208 for it = 1:60
209 mid = (lo + hi) / 2;
210 if cv_plate_det(n, nu, lo) * cv_plate_det(n, nu, mid) < 0
211 hi = mid;
212 else
213 lo = mid;
214 end
215 end
216 lam = [lam, (lo + hi) / 2];
217 end
218 prev = cur;
219end
220end
222function d = cv_plate_det(n, nu, x)
223% W = A J_n(lam r/a) + B I_n(lam r/a); the free edge sets the radial moment
224% and the Kelvin-Kirchhoff shear to zero, and the 2 x 2 determinant vanishes.
225J = besselj(n, x); I = besseli(n, x);
226Jp = besselj(n - 1, x) - (n / x) * J;
227Ip = besseli(n - 1, x) - (n / x) * I;
228Jpp = -Jp / x + (n^2 / x^2 - 1) * J;
229Ipp = -Ip / x + (n^2 / x^2 + 1) * I;
230a11 = Jpp + nu * (Jp / x - n^2 * J / x^2);
231a12 = Ipp + nu * (Ip / x - n^2 * I / x^2);
232a21 = -Jp - (1 - nu) * (n^2 / x^2) * (Jp - J / x);
233a22 = Ip - (1 - nu) * (n^2 / x^2) * (Ip - I / x);
234d = a11 * a22 - a12 * a21;
235end
237function [pd, pb] = cv_farfield(P, b, m, R, nth, tho, pho)
238% Far-field pressure of one mode, per unit modal velocity, computed by direct
239% summation over a theta-discretized dipole sheet and by the Bessel reduction.
240c = P.speed_of_sound;
241rho0 = P.air_density;
242om = b.omega(m);
243k = om / c;
244n = b.ndia(m);
245ld = b.ldip(m);
246rh = [sin(tho) * cos(pho), sin(tho) * sin(pho), cos(tho)];
247th = ((0:nth - 1) + 0.5) / nth * 2 * pi;
248dth = 2 * pi / nth;
249acc = 0;
250for e = 1:length(b.surfR)
251 x = b.surfR(e) * cos(th);
252 y = b.surfR(e) * sin(th);
253 z = b.surfZ(e) * ones(1, nth);
254 nx = -b.surfDz(e) * cos(th);
255 ny = -b.surfDz(e) * sin(th);
256 nz = b.surfDr(e) * ones(1, nth);
257 dx = R * rh(1) - x;
258 dy = R * rh(2) - y;
259 dz = R * rh(3) - z;
260 rr = sqrt(dx.^2 + dy.^2 + dz.^2);
261 proj = (nx .* dx + ny .* dy + nz .* dz) ./ rr * b.SW(e, m) .* cos(n * th);
262 V = ld * b.surfR(e) * b.surfLen(e) * dth;
263 acc = acc + sum((1i * om * rho0 / (4 * pi)) * V * proj .* (1i * k + 1 ./ rr) .* exp(-1i * k * rr) ./ rr);
264end
265pd = abs(acc);
266mu = cos(tho);
267sth = sin(tho);
268al = k * b.surfR * sth;
269Jn = besselj(n, al);
270Jp1 = besselj(n + 1, al);
271if n == 0
272 Jm1 = -Jp1;
273else
274 Jm1 = besselj(n - 1, al);
275end
276ip = 1i^n;
277Tk = ip * (2 * pi * Jn .* (b.surfDr * mu) - 1i * pi * (Jp1 - Jm1) .* (b.surfDz * sth));
278An = sum(ld * b.surfR .* b.surfLen .* b.SW(:, m) .* Tk .* exp(1i * (k * b.surfZ) * mu));
279pb = abs(-(om^2 * rho0 / (4 * pi * c * R)) * exp(-1i * k * R) * An * cos(n * pho));
280end
282function s = cv_verdict(ok)
283if ok
284 s = 'ok';
285else
286 s = 'MISMATCH';
287end
288end