Blame matlab/cymbal_basis.m at a136eac459c3f77b8dcb0e902940a6e3677b246f

4ddcc1b Cymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1 month ago1function basis = cymbal_basis(P)

2% CYMBAL_BASIS Mode set of the cymbal: one shell eigensolve per nodal-diameter

3% order, plus damping and the radiation data the later stages need.

5% Because the plate is a shell of revolution, the modes separate exactly in

6% the azimuth (see cymbal_shell.m), so the two-dimensional eigenproblem breaks

7% into n_diameters + 1 one-dimensional ones. Each is reduced to symmetric

8% standard form through the Cholesky factor of the consistent mass matrix and

9% diagonalized by svd rather than eig: the reduced matrix is symmetric

10% positive semidefinite, so its singular vectors are its eigenvectors, and svd

11% keeps them orthonormal for the closely spaced pairs this geometry produces.

12% The eigenvectors come back mass-normalized, phi' M phi = I. Rigid-body modes

13% (two at n = 0, two at n = 1, none above) sit at zero frequency and are

14% dropped with everything below min_mode_hz; modes above max_mode_hz are

15% truncated, and the broadband part they carry is reinstated by the direct

16% term in cymbal_radiate.m.

17%

18% Only the cos(n theta) member of each degenerate pair is retained. That is

19% exact for the model as posed: every contact sits on the meridian theta = 0,

7da0671 Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 1 month ago20% where the sin member has a node.

21%

22% HAMMERING. A perfectly axisymmetric plate is a bad cymbal, and the model has

23% to say so. Its modes are exact cos(n theta) shapes, and a mode with n nodal

24% diameters radiates only into the cone where its azimuthal trace speed is

25% supersonic, n < k a sin(theta). Below a few kilohertz that leaves only the

26% handful of n = 0, 1, 2 modes audible at any listening position, so an ideal

27% cymbal sounds like a small bell. Real cymbals are hammered and lathed, and

28% that is not decoration: the irregularity it leaves destroys the exact

29% separation, and every mode of the real plate is a hybrid of the ideal ones

30% that lie near it in frequency. A hybrid carries low-order content, so it

31% radiates.

32%

33% The hybridization is first-order degenerate perturbation theory, applied as

34% a seeded orthogonal mixing of modes lying within hammer_bw of each other in

35% frequency (a sequence of Givens rotations on neighbouring pairs, which is

36% orthogonal by construction and therefore conserves energy exactly), plus a

37% seeded relative detune of hammer_detune, which is what makes near neighbours

38% beat. The mixing matrix is stored as basis.mix and applied downstream

39% wherever the shapes enter linearly: to the stroke's modal projection in

40% cymbal_cascade.m and to the radiation transfer in cymbal_radiate.m. Because

41% the window is relative, the sparse low modes stay pure and only the dense

42% region above a kilohertz hybridizes, which is where a real cymbal's wash

43% lives. hammer_mix = 0 recovers the ideal axisymmetric plate.

4ddcc1b Cymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1 month ago44%

45% Damping has three parts.

46% internal damp_rate / om + 0.5 (loss_factor + loss_slope om). The first

47% term is a frequency-independent decay RATE, which is what a

48% resistive path out of the plate gives and what keeps the low

49% modes of a cymbal ringing for tens of seconds; the other two are

50% the usual constant and frequency-proportional loss factors, and

51% they are what makes the shimmer die well before the hum does;

52% mount a point dashpot of mount_loss N s/m on the axial motion at the

53% inner edge, standing in for the felt washer and the stand. With

54% mass-normalized shapes its modal damping is c phi_axial^2 / 2 om,

55% which falls off steeply with n because the shapes of high

56% nodal-diameter modes vanish at the hole. This is what kills the

57% low n = 0 gong modes a mounted cymbal does not sustain;

58% radiation computed from the actual radiated power, below.

59%

60% Radiated power is evaluated in the far field, where it is defined, rather

61% than by summing element powers as the sibling triangle does (an incoherent

62% element sum is mesh dependent for a continuous surface: halving the element

63% size halves the estimate). The plate is treated as a sheet of normal dipoles

64% of volume dS * ldip, and for a cos(n theta) source on a surface of revolution

65% the azimuthal integral of the far-field kernel is a Bessel function, so the

66% far field is exactly

67% p(th_o, ph_o) = -(om^2 rho0 / (4 pi c R)) e^{-ikR} A_n(th_o) cos(n ph_o)

68% A_n = ldip sum_e r_e L_e W_e e^{i k z_e cos th_o}

69% [ 2 pi i^n J_n(al) r'_e cos th_o

70% - pi (i^{n+1} J_{n+1}(al) + i^{n-1} J_{n-1}(al)) z'_e sin th_o ]

71% al = k r_e sin th_o

72% and the power integral over the sphere collapses to one quadrature in

73% cos th_o. The cos(n ph_o) factor is the model's multipole structure made

74% explicit: a mode with n nodal diameters is a 2^n-pole, so high-n modes

75% radiate very little at low frequency, and none of them radiate on the axis.

76%

77% ldip is the one phenomenological length in the radiation model: the distance

78% over which the air short-circuits from one face of the plate to the other.

79% It is taken as half the shorter of the flexural and acoustic wavelengths,

80% capped at the plate radius, and scaled by dipole_scale. Below the plate's

81% coincidence frequency (about 18 kHz for 1 mm bronze) it is the flexural

82% wavelength that sets it, which is why a cymbal is such a poor radiator in

83% the low register and gets steadily better as the cascade fills the top.

85G = cymbal_profile(P);

86E = P.youngs_modulus;

87nu = P.poisson_ratio;

88rho = P.density;

89c = P.speed_of_sound;

90rho0 = P.air_density;

91nn = G.nn;

92nel = nn - 1;

93Lref = G.S / nel;

94nmax = max(0, round(P.n_diameters));

96% ---- radiation sample points: two per element, Hermite-evaluated ----

97nsub = 2;

98nsurf = nel * nsub;

99surfS = zeros(nsurf, 1);

100surfLen = zeros(nsurf, 1);

101surfEl = zeros(nsurf, 1);

102surfXi = zeros(nsurf, 1);

103for e = 1:nel

104 L = G.snode(e + 1) - G.snode(e);

105 for j = 1:nsub

106 idx = (e - 1) * nsub + j;

107 xi = (j - 0.5) / nsub;

108 surfS(idx) = G.snode(e) + xi * L;

109 surfLen(idx) = L / nsub;

110 surfEl(idx) = e;

111 surfXi(idx) = xi;

112 end

113end

114surfR = reshape(interp1(G.s, G.r, surfS), nsurf, 1);

115surfZ = reshape(interp1(G.s, G.z, surfS), nsurf, 1);

116surfDr = reshape(interp1(G.s, G.drds, surfS), nsurf, 1);

117surfDz = reshape(interp1(G.s, G.dzds, surfS), nsurf, 1);

118

119% ---- eigensolve per circumferential order ----

120fs = [];

121nd = [];

122NUall = [];

123NWall = [];

124NTall = [];

125SWall = [];

126omLo = 2 * pi * P.min_mode_hz;

127omHi = 2 * pi * P.max_mode_hz;

128for n = 0:nmax

129 [K, M] = cymbal_shell(G, n, E, nu, rho, Lref);

130 Rc = chol(M);

131 Am = (Rc') \ (K / Rc);

132 Am = (Am + Am') / 2;

133 [Uv, Sv, ~] = svd(Am);

134 om2 = flipud(diag(Sv));

135 Vv = fliplr(Uv);

136 om = sqrt(max(0, om2)) * P.stiffness_trim;

137 keep = find(om > omLo & om < omHi);

138 if length(keep) == 0

139 continue;

140 end

141 Phi = Rc \ Vv(:, keep);

142 nk = length(keep);

143 Uz = Phi(1:4:4 * nn, :); % meridional, nn x nk

144 Wz = Phi(3:4:4 * nn, :); % normal

145 Tz = Phi(4:4:4 * nn, :) * Lref; % dW/ds, unscaled back to physical

146 % normal displacement at the radiation sample points (Hermite)

147 SW = zeros(nsurf, nk);

148 for idx = 1:nsurf

149 e = surfEl(idx);

150 xi = surfXi(idx);

151 L = G.snode(e + 1) - G.snode(e);

152 h1 = 1 - 3 * xi^2 + 2 * xi^3;

153 h2 = L * (xi - 2 * xi^2 + xi^3);

154 h3 = 3 * xi^2 - 2 * xi^3;

155 h4 = L * (-xi^2 + xi^3);

156 SW(idx, :) = h1 * Wz(e, :) + h2 * Tz(e, :) + h3 * Wz(e + 1, :) + h4 * Tz(e + 1, :);

157 end

158 fs = [fs; om(keep) / (2 * pi)];

159 nd = [nd; n * ones(nk, 1)];

160 NUall = [NUall, Uz];

161 NWall = [NWall, Wz];

162 NTall = [NTall, Tz];

163 SWall = [SWall, SW];

164end

165N = length(fs);

166if N == 0

167 error('No cymbal modes between min_mode_hz and max_mode_hz -- check geometry/material.');

168end

169[fs, order] = sort(fs);

170fs = reshape(fs, N, 1);

171nd = reshape(nd(order), N, 1);

172NUall = NUall(:, order);

173NWall = NWall(:, order);

174NTall = NTall(:, order);

175SWall = SWall(:, order);

7da0671 Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 1 month ago176

177% ---- hammering: detune, then hybridize near-degenerate modes ----

178hstate = mod(floor(abs(P.hammer_seed)), 4294967296);

179[hu, hstate] = cymbal_rng(hstate, 2 * N);

180hz = sqrt(-2 * log(hu(1:N))) .* cos(2 * pi * hu(N + 1:2 * N));

181fs = fs .* (1 + P.hammer_detune * hz);

182[fs, ord2] = sort(fs);

183fs = reshape(fs, N, 1);

184nd = reshape(nd(ord2), N, 1);

185NUall = NUall(:, ord2);

186NWall = NWall(:, ord2);

187NTall = NTall(:, ord2);

188SWall = SWall(:, ord2);

189

190HPASS = 4;

191mix = eye(N);

192if P.hammer_mix > 0 && N > 1

193 [hu2, hstate] = cymbal_rng(hstate, HPASS * (N - 1));

194 q = 0;

195 for pass = 1:HPASS

196 for i = 1:N - 1

197 q = q + 1;

198 if (fs(i + 1) - fs(i)) / fs(i) > P.hammer_bw

199 continue;

200 end

201 th = P.hammer_mix * (hu2(q) - 0.5) * pi;

202 ca = cos(th);

203 sa = sin(th);

204 ri = mix(i, :);

205 rj = mix(i + 1, :);

206 mix(i, :) = ca * ri - sa * rj;

207 mix(i + 1, :) = sa * ri + ca * rj;

208 end

209 end

210end

211mix2 = mix.^2;

212

4ddcc1b Cymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1 month ago213omega = 2 * pi * fs;

214

215% ---- the acoustic short-circuit length, per mode ----

216tMean = 2 * pi * trapz(G.s, G.r .* G.thk) / G.area;

217cL = sqrt(E / (rho * (1 - nu^2)));

218kap = tMean / sqrt(12);

219kFlex = sqrt(omega / (kap * cL));

220kAc = omega / c;

221ldip = P.dipole_scale * min(P.radius_m, pi ./ max(kFlex, kAc));

222

223% ---- radiated power, and the radiation damping that follows from it ----

224nMu = 24;

225mu = reshape(((1:nMu) - 0.5) / nMu * 2 - 1, 1, nMu); % cos(theta_o) midpoints

226dmu = 2 / nMu;

227sth = sqrt(max(0, 1 - mu.^2));

228zetaRad = zeros(N, 1);

229sigmaEff = zeros(N, 1);

230for m = 1:N

231 om = omega(m);

232 k = om / c;

233 nn_d = nd(m);

234 alpha = (k * surfR) * sth; % nsurf x nMu

235 Jn = besselj(nn_d, alpha);

236 Jp1 = besselj(nn_d + 1, alpha);

237 if nn_d == 0

238 Jm1 = -Jp1;

239 else

240 Jm1 = besselj(nn_d - 1, alpha);

241 end

242 % i^(n+1) J_{n+1} + i^(n-1) J_{n-1} = i^n * i * (J_{n+1} - J_{n-1})

243 ipow = 1i^nn_d;

244 Tk = ipow * (2 * pi * Jn .* (surfDr * mu) ...

245 - 1i * pi * (Jp1 - Jm1) .* (surfDz * sth));

246 carrier = exp(1i * (k * surfZ) * mu);

247 wgt = ldip(m) * surfR .* surfLen .* SWall(:, m);

248 An = sum((wgt * ones(1, nMu)) .* Tk .* carrier, 1); % 1 x nMu

249 Kn = sum(abs(An).^2) * dmu;

250 if nn_d == 0

251 piN = 2 * pi;

252 else

253 piN = pi;

254 end

255 zetaRad(m) = om^3 * rho0 * piN * Kn / (32 * pi^2 * c^3);

256 % radiation efficiency, for the diagnostics only

257 vv = sum(surfR .* surfLen .* SWall(:, m).^2) * piN;

258 sigmaEff(m) = 2 * zetaRad(m) * om / max(1e-30, rho0 * c * vv);

259end

260

261% ---- mount: a point dashpot on the axial motion at the inner edge ----

262axHole = NUall(1, :).' * G.dzdsnode(1) + NWall(1, :).' * G.drdsnode(1);

263zetaMount = P.mount_loss * axHole.^2 ./ (2 * omega);

264

266% over to the hybridized ones as the diagonal of the transformed damping

267% matrix, which to first order is the mixing weights squared.

268zetaRad = mix2 * zetaRad;

269zetaMount = mix2 * zetaMount;

270sigmaEff = mix2 * sigmaEff;

271

4ddcc1b Cymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1 month ago272zeta = P.damp_rate ./ omega ...

273 + 0.5 * (P.loss_factor + P.loss_slope * omega) ...

274 + P.radiation_loss * zetaRad + zetaMount;

275

276basis.n = N;

277basis.f = fs;

278basis.omega = omega;

279basis.zeta = zeta;

280basis.zetaRad = zetaRad;

281basis.zetaMount = zetaMount;

282basis.sigma = sigmaEff;

283basis.ndia = nd;

284basis.ldip = ldip;

7da0671 Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 1 month ago285basis.mix = mix;

4ddcc1b Cymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1 month ago286basis.NU = NUall;

287basis.NW = NWall;

288basis.NT = NTall;

289basis.rnode = reshape(G.rnode, nn, 1);

290basis.znode = reshape(G.znode, nn, 1);

291basis.snode = reshape(G.snode, nn, 1);

292basis.drdsnode = reshape(G.drdsnode, nn, 1);

293basis.dzdsnode = reshape(G.dzdsnode, nn, 1);

294basis.thknode = reshape(G.thknode, nn, 1);

295basis.csnode = reshape(interp1(G.s, G.cs, G.snode), nn, 1);

296basis.ctnode = reshape(interp1(G.s, G.ct, G.snode), nn, 1);

297basis.surfR = surfR;

298basis.surfZ = surfZ;

299basis.surfDr = surfDr;

300basis.surfDz = surfDz;

301basis.surfLen = surfLen;

302basis.SW = SWall;

303basis.S = G.S;

304basis.area = G.area;

305basis.mass = G.mass;

306basis.thkMean = tMean;

307end