4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1function [p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, fsHi)
2% CYMBAL_RADIATE Free-field pressure at the listener from the modal envelopes.
3%
4% Two contributions:
5%
6% 1. The mode bank. Each retained mode is a normal-dipole sheet on the plate:
7% a surface element of area dS carries a dipole of volume dS * ldip along
8% the local normal (ldip from cymbal_basis.m), and the element dipoles are
9% summed coherently with their relative retardations at the mode's centre
10% frequency,
11% G = sum_e (i om rho0 / 4 pi) V_e (u_e . rhat_e) (i k + 1/r_e)
12% exp(-i k (r_e - rcen)) / r_e,
13% over n_azimuth azimuths and two sample points per meridian element. The
14% common carrier exp(-i k rcen) is stripped, so relative phase across the
15% plate is kept while the overall propagation delay is not double counted.
16% Unlike the far-field form used for the radiated power, this keeps the
17% 1/r^2 near-field term and the true distances, which matters: a listener
18% half a metre from a 40 cm plate is nowhere near the far field. Only the
19% normal component of the surface motion radiates; the in-plane component
20% of a thin plate displaces almost no air.
21%
22% Because the azimuthal shape is exactly cos(n theta), the sum over
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 23% azimuths carries the multipole cancellation of the ideal mode: on the
24% axis of symmetry every mode with n >= 1 cancels exactly, and off axis a
25% mode radiates only into the cone where its azimuthal trace is supersonic,
26% n < k a sin(theta). That cutoff is why an unhammered plate is heard as a
27% handful of low-order partials; the hybridization applied below is what
28% gives the rest of the modes a way out. Moving the listener re-runs only
29% this stage.
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 30%
31% The synthesis itself is an oscillator bank rather than the closed-form
32% damped exponentials of the sibling repositories: the cascade makes every
33% amplitude time-varying, so each mode is a fixed-frequency oscillator
34% whose envelope is the cascade solution, interpolated linearly across each
35% control frame. Within a frame the two matrix products below are exact,
36% and the per-frame phase advance is a unit-modulus multiply, so nothing
37% drifts.
38%
39% Amplitude comes from the cascade; PHASE comes from a linear impulse
40% response carried alongside it. Each oscillator keeps a unit phasor and the
41% magnitude the mode would have had with no cascade at all; a contact adds
42% its signed modal velocity jump to that complex state and the phasor is
43% renormalized. A mode that was quiet is therefore re-phased by the stroke,
44% which is what makes the attack a coherent click, while a mode already
45% ringing loudly keeps the phase it had, which is also right. What is lost
46% is the intermodulation the real cascade would produce, since the energy
47% it moves into a mode arrives with whatever phase that mode already had.
48%
49% 2. Direct strike clicks. Modes above max_mode_hz are truncated; their
50% collective off-resonance response is the immediate local motion of the
51% plate under the tip. It is reinstated from the driving-point impedance of
52% an infinite plate, which for a plate is a pure resistance
53% Z = 8 sqrt(D rho t): the contact force divided by Z is the local
54% velocity, which radiates as a dipole over about half a flexural
55% wavelength, evaluated at the pulse's own centre frequency om* = pi / tau.
56% Scaled by direct_gain. This is the least principled term of the model,
57% and direct_gain = 0 turns it off.
58%
59% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.
61c = P.speed_of_sound;
62rho0 = P.air_density;
63N = basis.n;
64nsurf = length(basis.surfR);
65nth = max(8, round(P.n_azimuth));
66lx = P.listener_x_m;
67ly = P.listener_y_m;
68lz = P.listener_z_m;
70% ---- surface element geometry, meridian samples x azimuths ----
71th = reshape(((0:nth - 1) + 0.5) / nth * 2 * pi, 1, nth);
72dth = 2 * pi / nth;
73ct = cos(th);
74st = sin(th);
75ex = basis.surfR * ct; % nsurf x nth
76ey = basis.surfR * st;
77ez = basis.surfZ * ones(1, nth);
78nx = (-basis.surfDz) * ct;
79ny = (-basis.surfDz) * st;
80nz = basis.surfDr * ones(1, nth);
81dSa = (basis.surfR .* basis.surfLen) * (dth * ones(1, nth));
83dx = lx - ex;
84dy = ly - ey;
85dz = lz - ez;
86re = sqrt(dx.^2 + dy.^2 + dz.^2);
87listenerInside = min(min(re)) < 0.02;
88re = max(re, 0.01);
89cosg = (nx .* dx + ny .* dy + nz .* dz) ./ re; % nsurf x nth
91cx = 0;
92cy = 0;
93cz = sum(basis.surfR .* basis.surfLen .* basis.surfZ) / sum(basis.surfR .* basis.surfLen);
94rcen = sqrt((lx - cx)^2 + (ly - cy)^2 + (lz - cz)^2);
96% ---- per-mode listener transfer ----
97G = zeros(N, 1) + 0i;
98for m = 1:N
99 om = basis.omega(m);
100 k = om / c;
101 azi = cos(basis.ndia(m) * th); % 1 x nth
102 amp = (basis.ldip(m) * dSa .* (basis.SW(:, m) * ones(1, nth))) .* (ones(nsurf, 1) * azi);
103 ker = amp .* cosg .* (1i * k + 1 ./ re) .* exp(-1i * k * (re - rcen)) ./ re;
104 G(m) = (1i * om * rho0 / (4 * pi)) * sum(sum(ker));
105end
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 106% the plate's real modes are the hammered hybrids, and the transfer is linear
107% in the shape, so the mixing of cymbal_basis.m applies directly here
108G = basis.mix * G;
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 109
110% ---- oscillator bank driven by the cascade envelopes ----
111T = floor(tframe(end) * fsHi);
112p = zeros(1, T);
113nf = length(tframe) - 1;
114L = max(1, round((tframe(2) - tframe(1)) * fsHi));
115omd = basis.omega .* sqrt(max(1e-12, 1 - basis.zeta.^2));
116Eb = exp(1i * omd * ((0:L - 1) / fsHi)); % N x L, hoisted
117stepPh = exp(1i * omd * (L / fsHi)); % per-frame advance
118decLin = exp(-basis.zeta .* basis.omega * (L / fsHi));
119ramp = (0:L - 1) / L;
120sPh = ones(N, 1) + 0i;
121aLin = zeros(N, 1);
122pos = 1;
123for f = 1:nf
124 if pos > T
125 break;
126 end
127 hit = find(kf == f);
128 for q = 1:length(hit)
129 cc = aLin .* sPh + DV(:, hit(q));
130 aLin = abs(cc);
131 sPh = cc ./ max(1e-300, aLin);
132 end
133 len = min(L, T - pos + 1);
134 A0 = env(:, f);
135 A1 = env(:, f + 1);
136 c1 = (G .* A0 .* sPh).';
137 c2 = (G .* (A1 - A0) .* sPh).';
138 seg = real((c1 * Eb(:, 1:len)) + (c2 * Eb(:, 1:len)) .* ramp(1:len));
139 p(pos:pos + len - 1) = seg;
140 sPh = sPh .* stepPh;
141 aLin = aLin .* decLin;
142 pos = pos + L;
143end
145% ---- direct strike clicks ----
146E = P.youngs_modulus;
147nu = P.poisson_ratio;
148rho = P.density;
149cL = sqrt(E / (rho * (1 - nu^2)));
150I = length(tstr);
151ks = min(T, max(1, round(tstr * fsHi) + 1));
152ang = P.strike_angle_deg * pi / 180;
153directPeak = 0;
154for i = 1:I
155 [~, ~, thk, drds, dzds] = cymbal_at(basis, sstr(i));
156 [rs, zs] = cymbal_point(basis, sstr(i));
157 dvx = cos(ang) * dzds + sin(ang) * drds; % stroke direction at theta = 0
158 dvz = -cos(ang) * drds + sin(ang) * dzds;
159 rsx = lx - rs;
160 rsy = ly;
161 rsz = lz - zs;
162 rr = max(0.01, sqrt(rsx^2 + rsy^2 + rsz^2));
163 cg = (rsx * dvx + rsz * dvz) / rr;
164 nP = max(2, round(taustr(i) * fsHi));
165 % the click carries the true path delay relative to the carrier that was
166 % stripped from G, which is negative when the contact is nearer the
167 % listener than the plate centroid
168 k0 = max(1, ks(i) + round((rr - rcen) / c * fsHi));
169 if k0 + nP > T
170 continue;
171 end
172 D = E * thk^3 / (12 * (1 - nu^2));
173 Zp = 8 * sqrt(D * rho * thk);
174 omStar = pi / taustr(i);
175 kap = thk / sqrt(12);
176 kB = sqrt(omStar / (kap * cL));
177 lamB = 2 * pi / kB;
178 Aeff = pi * (lamB / 4)^2; % patch driven by the contact
179 ldir = pi / max(kB, omStar / c);
180 mu = Aeff * ldir / Zp; % dipole volume per unit force
181 amp = P.direct_gain * fstr(i) * mu * rho0 * omStar^2 * cg / (4 * pi * rr * c);
182 pulse = amp * sin(pi * (0:nP) / nP);
183 p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
184 if abs(amp) > directPeak
185 directPeak = abs(amp);
186 end
187end
189radinfo.gainAbs = abs(G);
190radinfo.directPeakPa = directPeak;
191radinfo.listenerInside = listenerInside;
192radinfo.listenerDistM = rcen;
193end