4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1function [env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, taustr, sstr, duration)
2% CYMBAL_CASCADE Modal excitation and the nonlinear flow of energy up the
3% frequency ladder.
4%
5% This is the stage that makes the instrument a cymbal rather than a gong or a
6% bell. A cymbal struck hard is not a linear resonator: the transverse
7% displacement is comparable to the plate thickness, so the von Karman
8% stretching of the mid-surface couples the modes, and energy injected in the
9% low modes by the stick migrates upward over tens of milliseconds. That is
10% what the ear hears as a crash swelling into shimmer after its attack, and it
11% is why the same cymbal struck softly just pings: the coupling strength goes
12% as the square of the amplitude, so a quiet stroke stays linear.
13%
14% Integrating the von Karman equations is out of reach here (the cubic
15% coupling is a dense N^3 tensor, and the sibling repositories render in
16% seconds). What is integrated instead is the ENERGY of each mode. Each mode
17% is given an up-target, the lowest retained mode at least cascade_ratio times
18% its frequency, and a down-target the same ratio below, and energy performs a
19% biased random walk between them:
20% dE_m/dt = -(1 + beta) Gam_m E_m + (arrivals from below and above)
21% - 2 zeta_m om_m E_m
22% Gam_m = cascade_gain om_m (w_m / t_mean)^2 = cascade_gain E_m / (om_m M t^2)
23% w_m^2 = E_m / (om_m^2 M)
24% with beta = cascade_back < 1 the downward-to-upward ratio, so the net flux is
25% toward high frequency. Steps are taken in log frequency rather than in mode
26% index on purpose: four-wave interactions move energy by a fraction of an
27% octave at a time, not by one mode, and a ladder defined on the index would
28% make the swell depend on how densely the band happens to be resolved. The
29% rate is proportional to om_m because the higher a mode sits the faster its
30% interactions run, and to the square of that mode's own
31% amplitude-to-thickness ratio because that is the order of the geometric
32% nonlinearity. Taking the amplitude locally rather than as a plate-wide rms
33% is what keeps the cascade honest at the top of the band: the flux out of a
34% mode goes as the SQUARE of its energy, so a nearly empty mode passes nothing
35% on, the front of the cascade advances rather than the whole band draining at
36% once, and the swell stops on its own as the energy thins out.
37%
38% The scheme conserves energy exactly when the damping is off and the top of
39% the ladder reflects (cascade_leak = 0), which the validation checks; with
40% cascade_leak = 1 the flux out of the top is lost, standing in for the
41% truncated modes above max_mode_hz. cascade_gain is the one fitted number in
42% the model: it is a dimensionless coupling constant, set so that a default
43% crash takes about the 40 ms a real one takes to open up.
44%
3a32a66Default the turbulent fluctuation offJeremy Magland 45% FLUCTUATION, off by default. What the equations above track is the MEAN
46% energy of each mode, and the mean is not quite what a listener hears: a mode
47% sitting in a turbulent field is continually being filled and drained by its
48% neighbours, and its amplitude wanders about that mean on the timescale of
49% the exchange. Restoring that wander broadens each partial and raises the
50% modal overlap, which on paper is what turns a few hundred audible partials
51% into a continuum. In practice it does not sound like a cymbal. The
52% decorrelation rate this model assigns runs up to the control rate itself, so
53% what it adds is closer to modulation grit than to a wash, and turbulence is
54% therefore 0 at the defaults. It is kept because the physics is right even
55% where this discretization of it is not, and because it is the honest place
56% to start if the graininess is ever worth chasing; small values, a few tenths,
57% are the gentle end.
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 58%
3a32a66Default the turbulent fluctuation offJeremy Magland 59% Each mode carries a log-amplitude that follows an Ornstein-Uhlenbeck
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 60% process discretized at the control rate,
61% rho_m = exp(-(energy exchanged this step) / E_m)
62% g_m <- rho_m g_m + sqrt(1 - rho_m^2) turbulence z_m
63% with the correlation rate taken as the mode's own fractional energy turnover
64% rather than as a fitted time: a mode being fed while nearly empty
65% decorrelates completely every step, a mode merely ringing down does not
66% decorrelate at all. A running variance keeps exp(g) unbiased, so the
67% fluctuation redistributes energy in time without creating or destroying any,
68% and the process is seeded, so a render is reproducible. Where the cascade is
69% off the turnover is zero, rho is one, g stays where it started, and the
70% envelopes are exactly the closed-form exponentials of a linear plate.
71%
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 72% What this buys is a cymbal that responds to dynamics the way the instrument
73% does, at the cost of one control-rate loop rather than a full nonlinear
74% integration. What it gives up is phase: energy moves between modes but the
75% intermodulation products of a real cascade, and the split-pair beating of an
76% out-of-round cymbal, are not reproduced. cascade_gain = 0 turns the whole
77% mechanism off and leaves an exactly linear plate.
78%
79% Each contact adds ½ (dv_m)^2 to the energy of every mode, where
80% dv_m = J S(om_m tau) (phi_m(x_c) . dhat)
81% with S(x) = cos(x/2) / (1 - (x/pi)^2) the spectrum of the half-sine contact
82% pulse and dhat the stroke direction (strike_angle_deg away from the surface
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 83% normal, in the meridian plane). The projection is taken through the
84% hammering mixing of cymbal_basis.m, since the modes the stick actually
85% drives are the hybrids, not the ideal cos(n theta) shapes. Energies add incoherently, which is right
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 86% for a plate already ringing and exact for the first stroke into a silent one.
87%
88% Note that only the MAGNITUDE of the excitation survives into the energy
89% state. The signed jumps DV and the frame each contact lands in are returned
90% alongside, because cymbal_radiate.m needs them for the oscillator phases:
91% amplitude comes from the cascade, phase from the linear impulse response.
92%
93% Returns env (N x nf+1, modal VELOCITY amplitude at each frame boundary),
94% tframe (1 x nf+1), DV (N x number of contacts), kf (the frame index of each
95% contact) and cinfo diagnostics.
97N = basis.n;
98omega = basis.omega;
99zeta = basis.zeta;
100h = max(1e-4, P.cascade_step_ms * 1e-3);
101nf = max(1, ceil(duration / h));
102tframe = (0:nf) * h;
104% ---- modal velocity jump of each contact ----
105I = length(tstr);
106DV = zeros(N, I);
107ang = P.strike_angle_deg * pi / 180;
108for i = 1:I
109 [uh, wh] = cymbal_at(basis, sstr(i));
110 proj = -cos(ang) * wh + sin(ang) * uh; % 1 x N, along the stroke
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 111 proj = (basis.mix * reshape(proj, N, 1)).'; % hammered modes are hybrids
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 112 x = omega * taustr(i);
113 den = 1 - (x / pi).^2;
114 S = cos(x / 2) ./ den;
115 sing = abs(den) < 1e-6;
116 S(sing) = pi / 4; % limit at x -> pi
117 DV(:, i) = jstr(i) * S .* reshape(proj, N, 1);
118end
119kf = min(nf, max(1, floor(tstr / h) + 1)); % frame each contact lands in
121decay = exp(-2 * zeta .* omega * h);
122gmax = 0.5 / (h * (1 + P.cascade_back));
123Mtot = basis.mass;
124teff = basis.thkMean;
126% ---- the log-frequency ladder ----
127% upIdx(m) is the lowest mode at least cascade_ratio times f(m); modes with no
128% such partner point at N + 1, the sink standing in for the truncated band.
129% dnIdx(m) is the highest mode at most f(m) / cascade_ratio, or m itself at the
130% bottom, which reflects.
131ratio = max(1.02, P.cascade_ratio);
132f = basis.f;
133upIdx = zeros(N, 1);
134dnIdx = zeros(N, 1);
135for m = 1:N
136 ju = find(f > f(m) * ratio);
137 if length(ju) == 0
138 upIdx(m) = N + 1;
139 else
140 upIdx(m) = ju(1);
141 end
142 jd = find(f < f(m) / ratio);
143 if length(jd) == 0
144 dnIdx(m) = m;
145 else
146 dnIdx(m) = jd(length(jd));
147 end
148end
150env = zeros(N, nf + 1);
151Ee = zeros(N, 1);
152peakStrength = 0;
153peakWrms = 0;
154lost = 0;
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 155
156% one independent LCG stream per mode, for the amplitude fluctuation
157sig = P.turbulence;
158st1 = mod(floor(abs(P.seed)) * 2654435761 + (1:N).' * 40503, 4294967296);
159st2 = mod(st1 * 69069 + 12345, 4294967296);
160gl = zeros(N, 1);
161gv = zeros(N, 1);
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 163for f = 1:nf
164 hit = find(kf == f);
165 for q = 1:length(hit)
166 Ee = Ee + 0.5 * DV(:, hit(q)).^2;
167 end
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 168 env(:, f) = sqrt(2 * Ee) .* exp(gl - gv / 2);
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 169
170 wrms2 = sum(Ee ./ omega.^2) / Mtot;
171 strength = wrms2 / teff^2;
172 if strength > peakStrength
173 peakStrength = strength;
174 peakWrms = sqrt(wrms2);
175 end
176 gam = min(gmax, P.cascade_gain * Ee ./ (omega * Mtot * teff^2));
177 up = gam .* Ee * h;
178 dn = P.cascade_back * gam .* Ee * h;
179 arrU = accumarray(upIdx, up, [N + 1, 1]);
180 arrD = accumarray(dnIdx, dn, [N, 1]);
181 leak = arrU(N + 1);
182 En = Ee - up - dn + arrU(1:N) + arrD;
183 En(N) = En(N) + (1 - P.cascade_leak) * leak;
184 lost = lost + P.cascade_leak * leak;
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 185
186 if sig > 0
187 % fractional energy turnover this step sets the decorrelation
188 moved = up + dn + arrU(1:N) + arrD;
189 rho = exp(-min(30, moved ./ max(Ee, 1e-300)));
190 [u1, st1] = cymbal_rngv(st1);
191 [u2, st2] = cymbal_rngv(st2);
192 zg = sqrt(-2 * log(u1)) .* cos(2 * pi * u2);
193 w = 1 - rho.^2;
194 gl = rho .* gl + sqrt(w) * sig .* zg;
195 gv = rho.^2 .* gv + w * sig^2;
196 end
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 198 Ee = En .* decay;
199end
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 200env(:, nf + 1) = sqrt(2 * Ee) .* exp(gl - gv / 2);
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 201
202cinfo.frames = nf;
203cinfo.stepS = h;
204cinfo.peakStrength = peakStrength;
205cinfo.peakWrmsUm = peakWrms * 1e6;
206cinfo.cascadeLostJ = lost;
207cinfo.injectedJ = sum(sum(0.5 * DV.^2));
208end