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. It is a statement about the IDEAL plate, so
82% the hammering is switched off: hybridizing the modes is precisely what gives
83% the high-order ones a way out, and with it on they radiate on the axis like
84% anything else.
85Q = P;
86Q.hammer_mix = 0; Q.hammer_detune = 0;
87Q.n_diameters = 6; Q.max_mode_hz = 1500;
88Q.listener_x_m = 0; Q.listener_y_m = 0; Q.listener_z_m = 0.6;
89b6 = cymbal_basis(Q);
90[tstr, ~, jstr, taustr, fstr, sstr, sinf] = cymbal_strikes(Q, b6);
91[env, tfr, DV, kf] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.05);
92[~, ri] = cymbal_radiate(Q, b6, env, tfr, DV, kf, tstr, taustr, fstr, sstr, Q.sample_rate);
93g0 = max(ri.gainAbs(b6.ndia == 0));
94g1 = max(ri.gainAbs(b6.ndia >= 1));
95ok = g1 / g0 < 1e-10;
96fails = fails + (~ok);
97fprintf('on axis : |G| of n>=1 modes is %.1e of the n=0 modes %s\n', g1 / g0, cv_verdict(ok));
99% ---- 4b. the hammering mixing is exactly orthogonal ----
100% The hybridization is a product of Givens rotations, so it must be orthogonal
101% to round-off. That is what makes it energy conserving: it redistributes the
102% stroke among the modes and the radiation among them again, without inventing
103% or destroying any.
104bh = cymbal_basis(P);
105eMix = max(max(abs(bh.mix * bh.mix.' - eye(bh.n))));
106spread = sum(bh.mix.^2 > 0.01, 2);
107ok = eMix < 1e-12 && median(spread) > 2;
108fails = fails + (~ok);
109fprintf('hammering : orthogonality %.1e, median %d ideal modes per hybrid %s\n', ...
110 eMix, round(median(spread)), cv_verdict(ok));
112% ---- 5. Hertz contact-time scaling tau ~ v^(-1/5) ----
113% With a stick much lighter than the plate mass engaged during the contact the
114% reduced mass is constant and the free-sphere scaling must emerge.
115Q = P;
116Q.strike_count = 1;
117Q.stick_mass_kg = 2e-5;
118[~, ~, ~, tau1] = cymbal_strikes(Q, b6);
119Q.strike_velocity = 2 * P.strike_velocity;
120[~, ~, ~, tau2] = cymbal_strikes(Q, b6);
121ratio = tau2(1) / tau1(1);
122ok = abs(ratio - 2^(-1 / 5)) < 0.03;
123fails = fails + (~ok);
124fprintf('hertz : tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), cv_verdict(ok));
126% ---- 6. the cascade conserves energy, and switches off cleanly ----
127% With the damping removed, the top of the ladder reflecting and the
128% fluctuation off, the transfer scheme is a Markov generator on the modal
129% energies and must conserve their sum to round-off. With cascade_gain = 0
130% nothing is exchanged, so the fluctuation freezes at zero on its own and
131% every envelope must collapse onto the closed-form exp(-zeta om t) of an
132% uncoupled damped oscillator, fluctuation still switched on.
133Q = P;
134Q.n_diameters = 6; Q.max_mode_hz = 1500; Q.cascade_leak = 0; Q.turbulence = 0;
135bz = b6;
136bz.zeta = zeros(bz.n, 1);
137[tstr, ~, jstr, taustr, ~, sstr] = cymbal_strikes(Q, bz);
138[envC, ~, DVc] = cymbal_cascade(Q, bz, tstr, jstr, taustr, sstr, 0.4);
139Etot = sum(envC.^2 / 2, 1);
140eCons = max(abs(Etot / Etot(1) - 1));
141Q.cascade_gain = 0;
142Q.turbulence = 0.9; % explicitly on, so the freeze is actually tested
143[envL, tfrL] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.4);
144pred = envL(:, 1) * exp(-(b6.zeta .* b6.omega).' * 0 + 0);
145pred = (envL(:, 1) * ones(1, length(tfrL))) .* exp(-(b6.zeta .* b6.omega) * tfrL);
146sel = envL(:, 1) > 1e-12;
147eLin = max(max(abs(envL(sel, :) ./ pred(sel, :) - 1)));
148ok = eCons < 1e-10 && eLin < 1e-9;
149fails = fails + (~ok);
150fprintf('cascade : energy drift %.1e, linear limit err %.1e %s\n', eCons, eLin, cv_verdict(ok));
152% ---- 7. far-field falloff: -6 dB per doubling of distance ----
153% RMS over the render, with the two listeners on the same ray from the origin
154% so every element dipole keeps its directivity factor.
155Q = P;
156Q.normalize = 0; Q.n_diameters = 12; Q.max_mode_hz = 3000;
157Q.tail_s = 0.4; Q.oversample = 1; Q.n_azimuth = 96;
158Q.listener_x_m = 0.8; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;
159y1 = cymbal_render(Q);
160Q.listener_x_m = 1.6; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;
161y2 = cymbal_render(Q);
162drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
163ok = abs(drop + 6.02) < 0.25;
164fails = fails + (~ok);
165fprintf('falloff : 2 m -> 4 m %.2f dB rms (target -6.02) %s\n', drop, cv_verdict(ok));
167% ---- 8. the seeded ride figure is exactly reproducible ----
168Q = P;
169Q.strike_count = 8;
170[t1, v1, j1, ~, ~, s1, g1] = cymbal_strikes(Q, b6);
171[t2, v2, ~, ~, ~, s2] = cymbal_strikes(Q, b6);
172same = max(abs(t1 - t2)) == 0 && max(abs(v1 - v2)) == 0 && max(abs(s1 - s2)) == 0;
173ok = same ...
174 && g1.count == 8 ...
175 && abs(g1.meanV - 3.111518886) < 2e-6 ...
176 && abs(sum(j1) - 0.150818144) < 1e-9 ...
177 && abs(t1(2) - 0.159427523) < 1e-9;
178fails = fails + (~ok);
179fprintf('ride (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
180 Q.seed, g1.count, g1.meanV, sum(j1), cv_verdict(ok));
182% ---- 9. pinned end-to-end diagnostics at the defaults ----
183% The rms is pinned alongside the peak deliberately. The peak of a sum of
184% several hundred oscillators is a fragile thing to compare across machines,
185% where the rms is not, so a mismatch on the peak alone says the phases have
186% moved and a mismatch on both says the physics has.
187[y, basis, dinfo] = cymbal_render(P);
188yrms = sqrt(sum(y.^2) / length(y));
189ok = abs(dinfo.peakPa - 13.6155) < 0.02 ...
190 && abs(yrms - 0.04652826) < 2e-6 ...
191 && abs(basis.f(1) - 49.4130) < 0.02 ...
192 && basis.n == 378;
193fails = fails + (~ok);
194fprintf('pinned : peak %.4f Pa (target %.4f), rms %.8f (target %.8f), mode1 %.3f Hz, %d modes %s\n', ...
195 dinfo.peakPa, 13.6155, yrms, 0.04652826, basis.f(1), basis.n, cv_verdict(ok));
197% ---- 10. the wash decays ----
198% A cymbal rings for seconds, so the tail is not silent, but the energy in the
199% last half second must sit well below that of the opening half second.
200fs = P.sample_rate;
201n = length(y);
202head = y(1:round(0.5 * fs));
203tail = y(n - round(0.5 * fs):n);
204dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
205ok = dropE < -20 && dropE > -80;
206fails = fails + (~ok);
207fprintf('decay : last/first half second %.1f dB %s\n', dropE, cv_verdict(ok));
209if fails > 0
210 fprintf('%d CHECK(S) FAILED\n', fails);
211else
212 fprintf('all checks passed\n');
213end
215function f = cv_freqs(G, n, P, Lref)
216[K, M] = cymbal_shell(G, n, P.youngs_modulus, P.poisson_ratio, P.density, Lref);
217Rc = chol(M);
218Am = (Rc') \ (K / Rc);
219Am = (Am + Am') / 2;
220[~, Sv, ~] = svd(Am);
221f = sort(sqrt(diag(Sv))) / (2 * pi);
222end
224function lam = cv_plate_roots(n, nu)
225% Roots of the free-edge frequency equation of a circular plate, bracketed on
226% a coarse scan and bisected.
227lam = [];
228prev = cv_plate_det(n, nu, 0.30);
229for x = 0.30:0.01:26
230 cur = cv_plate_det(n, nu, x + 0.01);
231 if prev * cur < 0
232 lo = x;
233 hi = x + 0.01;
234 for it = 1:60
235 mid = (lo + hi) / 2;
236 if cv_plate_det(n, nu, lo) * cv_plate_det(n, nu, mid) < 0
237 hi = mid;
238 else
239 lo = mid;
240 end
241 end
242 lam = [lam, (lo + hi) / 2];
243 end
244 prev = cur;
245end
246end
248function d = cv_plate_det(n, nu, x)
249% W = A J_n(lam r/a) + B I_n(lam r/a); the free edge sets the radial moment
250% and the Kelvin-Kirchhoff shear to zero, and the 2 x 2 determinant vanishes.
251J = besselj(n, x); I = besseli(n, x);
252Jp = besselj(n - 1, x) - (n / x) * J;
253Ip = besseli(n - 1, x) - (n / x) * I;
254Jpp = -Jp / x + (n^2 / x^2 - 1) * J;
255Ipp = -Ip / x + (n^2 / x^2 + 1) * I;
256a11 = Jpp + nu * (Jp / x - n^2 * J / x^2);
257a12 = Ipp + nu * (Ip / x - n^2 * I / x^2);
258a21 = -Jp - (1 - nu) * (n^2 / x^2) * (Jp - J / x);
259a22 = Ip - (1 - nu) * (n^2 / x^2) * (Ip - I / x);
260d = a11 * a22 - a12 * a21;
261end
263function [pd, pb] = cv_farfield(P, b, m, R, nth, tho, pho)
264% Far-field pressure of one mode, per unit modal velocity, computed by direct
265% summation over a theta-discretized dipole sheet and by the Bessel reduction.
266c = P.speed_of_sound;
267rho0 = P.air_density;
268om = b.omega(m);
269k = om / c;
270n = b.ndia(m);
271ld = b.ldip(m);
272rh = [sin(tho) * cos(pho), sin(tho) * sin(pho), cos(tho)];
273th = ((0:nth - 1) + 0.5) / nth * 2 * pi;
274dth = 2 * pi / nth;
275acc = 0;
276for e = 1:length(b.surfR)
277 x = b.surfR(e) * cos(th);
278 y = b.surfR(e) * sin(th);
279 z = b.surfZ(e) * ones(1, nth);
280 nx = -b.surfDz(e) * cos(th);
281 ny = -b.surfDz(e) * sin(th);
282 nz = b.surfDr(e) * ones(1, nth);
283 dx = R * rh(1) - x;
284 dy = R * rh(2) - y;
285 dz = R * rh(3) - z;
286 rr = sqrt(dx.^2 + dy.^2 + dz.^2);
287 proj = (nx .* dx + ny .* dy + nz .* dz) ./ rr * b.SW(e, m) .* cos(n * th);
288 V = ld * b.surfR(e) * b.surfLen(e) * dth;
289 acc = acc + sum((1i * om * rho0 / (4 * pi)) * V * proj .* (1i * k + 1 ./ rr) .* exp(-1i * k * rr) ./ rr);
290end
291pd = abs(acc);
292mu = cos(tho);
293sth = sin(tho);
294al = k * b.surfR * sth;
295Jn = besselj(n, al);
296Jp1 = besselj(n + 1, al);
297if n == 0
298 Jm1 = -Jp1;
299else
300 Jm1 = besselj(n - 1, al);
301end
302ip = 1i^n;
303Tk = ip * (2 * pi * Jn .* (b.surfDr * mu) - 1i * pi * (Jp1 - Jm1) .* (b.surfDz * sth));
304An = sum(ld * b.surfR .* b.surfLen .* b.SW(:, m) .* Tk .* exp(1i * (k * b.surfZ) * mu));
305pb = abs(-(om^2 * rho0 / (4 * pi * c * R)) * exp(-1i * k * R) * An * cos(n * pho));
306end
308function s = cv_verdict(ok)
309if ok
310 s = 'ok';
311else
312 s = 'MISMATCH';
313end
314end