Blame matlab/cymbal_basis.m at 048e093d0a5deb889f21b4f702a07378fb46f994

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);

f61a91b Canonicalize the eigenvector signsJeremy Magland 1 month ago143 % svd fixes each singular vector only up to sign, and different linear

144 % algebra libraries make different choices. On its own that would not

145 % matter, since flipping one mode flips both its excitation and its

146 % radiation and the two cancel; but the hammering below forms SUMS of

147 % modes, where they do not, and the render would then depend on which

148 % machine it was computed on. Canonicalize instead: the largest component

149 % of every shape is positive, with ties broken by index so that the rule

150 % itself is deterministic.

151 for k = 1:nk

152 av = abs(Phi(:, k));

153 hit = find(av >= max(av) * (1 - 1e-9));

154 if Phi(hit(1), k) < 0

155 Phi(:, k) = -Phi(:, k);

156 end

157 end

4ddcc1b Cymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1 month ago158 Uz = Phi(1:4:4 * nn, :); % meridional, nn x nk

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

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

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

162 SW = zeros(nsurf, nk);

163 for idx = 1:nsurf

164 e = surfEl(idx);

165 xi = surfXi(idx);

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

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

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

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

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

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

172 end

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

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

175 NUall = [NUall, Uz];

176 NWall = [NWall, Wz];

177 NTall = [NTall, Tz];

178 SWall = [SWall, SW];

179end

180N = length(fs);

181if N == 0

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

183end

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

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

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

187NUall = NUall(:, order);

188NWall = NWall(:, order);

189NTall = NTall(:, order);

190SWall = 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 ago191

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

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

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

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

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

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

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

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

200NUall = NUall(:, ord2);

201NWall = NWall(:, ord2);

202NTall = NTall(:, ord2);

203SWall = SWall(:, ord2);

204

205HPASS = 4;

206mix = eye(N);

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

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

209 q = 0;

210 for pass = 1:HPASS

211 for i = 1:N - 1

212 q = q + 1;

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

214 continue;

215 end

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

217 ca = cos(th);

218 sa = sin(th);

219 ri = mix(i, :);

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

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

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

223 end

224 end

225end

226mix2 = mix.^2;

227

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

229

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

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

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

233kap = tMean / sqrt(12);

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

235kAc = omega / c;

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

237

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

239nMu = 24;

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

241dmu = 2 / nMu;

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

243zetaRad = zeros(N, 1);

244sigmaEff = zeros(N, 1);

245for m = 1:N

246 om = omega(m);

247 k = om / c;

248 nn_d = nd(m);

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

250 Jn = besselj(nn_d, alpha);

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

252 if nn_d == 0

253 Jm1 = -Jp1;

254 else

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

256 end

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

258 ipow = 1i^nn_d;

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

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

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

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

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

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

265 if nn_d == 0

266 piN = 2 * pi;

267 else

268 piN = pi;

269 end

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

271 % radiation efficiency, for the diagnostics only

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

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

274end

275

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

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

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

279

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

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

283zetaRad = mix2 * zetaRad;

284zetaMount = mix2 * zetaMount;

285sigmaEff = mix2 * sigmaEff;

286

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

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

289 + P.radiation_loss * zetaRad + zetaMount;

290

291basis.n = N;

292basis.f = fs;

293basis.omega = omega;

294basis.zeta = zeta;

295basis.zetaRad = zetaRad;

296basis.zetaMount = zetaMount;

297basis.sigma = sigmaEff;

298basis.ndia = nd;

299basis.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 ago300basis.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 ago301basis.NU = NUall;

302basis.NW = NWall;

303basis.NT = NTall;

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

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

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

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

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

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

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

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

312basis.surfR = surfR;

313basis.surfZ = surfZ;

314basis.surfDr = surfDr;

315basis.surfDz = surfDz;

316basis.surfLen = surfLen;

317basis.SW = SWall;

318basis.S = G.S;

319basis.area = G.area;

320basis.mass = G.mass;

321basis.thkMean = tMean;

322end