/ simulated_instruments / cymbal
Sign in
simulated_instruments / cymbal
7da0671
Branches
main
cymbal / matlab / cymbal_cascade.m
206 lines · 9.0 KBCodeBlameHistory
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.
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.
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.
7da0671Hammer the cymbal: hybridize near-degenerate modes, add the turbulent fluctuation about the mean cascade, and flatten the decay across the bandJeremy Magland 45% FLUCTUATION. What the equations above track is the MEAN energy of each mode,
46% and the mean is not what a listener hears. A mode sitting in a turbulent
47% field does not hold a smooth amplitude; it is continually being filled and
48% drained by its neighbours, and its amplitude wanders about that mean on the
49% timescale of the exchange. This matters more than it sounds: a few hundred
50% modes with smooth envelopes are a few hundred audible partials, spaced tens
51% of hertz apart with bandwidths of a couple of hertz, which is a chime. Give
52% each one a wandering amplitude with a correlation time of a few milliseconds
53% and each partial broadens by a couple of hundred hertz, the modal overlap
54% passes one, and the comb closes into the continuum a wash actually is.
56% So each mode carries a log-amplitude that follows an Ornstein-Uhlenbeck
57% process discretized at the control rate,
58% rho_m = exp(-(energy exchanged this step) / E_m)
59% g_m <- rho_m g_m + sqrt(1 - rho_m^2) turbulence z_m
60% with the correlation rate taken as the mode's own fractional energy turnover
61% rather than as a fitted time: a mode being fed while nearly empty
62% decorrelates completely every step, a mode merely ringing down does not
63% decorrelate at all. A running variance keeps exp(g) unbiased, so the
64% fluctuation redistributes energy in time without creating or destroying any,
65% and the process is seeded, so a render is reproducible. Where the cascade is
66% off the turnover is zero, rho is one, g stays where it started, and the
67% envelopes are exactly the closed-form exponentials of a linear plate.
68% turbulence = 0 removes the fluctuation and leaves the mean-energy solve.
4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 70% What this buys is a cymbal that responds to dynamics the way the instrument
71% does, at the cost of one control-rate loop rather than a full nonlinear
72% integration. What it gives up is phase: energy moves between modes but the
73% intermodulation products of a real cascade, and the split-pair beating of an
74% out-of-round cymbal, are not reproduced. cascade_gain = 0 turns the whole
75% mechanism off and leaves an exactly linear plate.
77% Each contact adds ½ (dv_m)^2 to the energy of every mode, where
78% dv_m = J S(om_m tau) (phi_m(x_c) . dhat)
79% with S(x) = cos(x/2) / (1 - (x/pi)^2) the spectrum of the half-sine contact
80% pulse and dhat the stroke direction (strike_angle_deg away from the surface
82% hammering mixing of cymbal_basis.m, since the modes the stick actually
83% 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 84% for a plate already ringing and exact for the first stroke into a silent one.
86% Note that only the MAGNITUDE of the excitation survives into the energy
87% state. The signed jumps DV and the frame each contact lands in are returned
88% alongside, because cymbal_radiate.m needs them for the oscillator phases:
89% amplitude comes from the cascade, phase from the linear impulse response.
91% Returns env (N x nf+1, modal VELOCITY amplitude at each frame boundary),
92% tframe (1 x nf+1), DV (N x number of contacts), kf (the frame index of each
93% contact) and cinfo diagnostics.
95N = basis.n;
96omega = basis.omega;
97zeta = basis.zeta;
98h = max(1e-4, P.cascade_step_ms * 1e-3);
99nf = max(1, ceil(duration / h));
100tframe = (0:nf) * h;
102% ---- modal velocity jump of each contact ----
103I = length(tstr);
104DV = zeros(N, I);
105ang = P.strike_angle_deg * pi / 180;
106for i = 1:I
107 [uh, wh] = cymbal_at(basis, sstr(i));
108 proj = -cos(ang) * wh + sin(ang) * uh; % 1 x N, along the stroke
111 den = 1 - (x / pi).^2;
112 S = cos(x / 2) ./ den;
113 sing = abs(den) < 1e-6;
114 S(sing) = pi / 4; % limit at x -> pi
115 DV(:, i) = jstr(i) * S .* reshape(proj, N, 1);
116end
117kf = min(nf, max(1, floor(tstr / h) + 1)); % frame each contact lands in
119decay = exp(-2 * zeta .* omega * h);
120gmax = 0.5 / (h * (1 + P.cascade_back));
121Mtot = basis.mass;
122teff = basis.thkMean;
124% ---- the log-frequency ladder ----
125% upIdx(m) is the lowest mode at least cascade_ratio times f(m); modes with no
126% such partner point at N + 1, the sink standing in for the truncated band.
127% dnIdx(m) is the highest mode at most f(m) / cascade_ratio, or m itself at the
128% bottom, which reflects.
129ratio = max(1.02, P.cascade_ratio);
130f = basis.f;
131upIdx = zeros(N, 1);
132dnIdx = zeros(N, 1);
133for m = 1:N
134 ju = find(f > f(m) * ratio);
135 if length(ju) == 0
136 upIdx(m) = N + 1;
137 else
138 upIdx(m) = ju(1);
139 end
140 jd = find(f < f(m) / ratio);
141 if length(jd) == 0
142 dnIdx(m) = m;
143 else
144 dnIdx(m) = jd(length(jd));
145 end
146end
148env = zeros(N, nf + 1);
149Ee = zeros(N, 1);
150peakStrength = 0;
151peakWrms = 0;
152lost = 0;
154% one independent LCG stream per mode, for the amplitude fluctuation
155sig = P.turbulence;
156st1 = mod(floor(abs(P.seed)) * 2654435761 + (1:N).' * 40503, 4294967296);
157st2 = mod(st1 * 69069 + 12345, 4294967296);
158gl = zeros(N, 1);
159gv = zeros(N, 1);
162 hit = find(kf == f);
163 for q = 1:length(hit)
164 Ee = Ee + 0.5 * DV(:, hit(q)).^2;
165 end
168 wrms2 = sum(Ee ./ omega.^2) / Mtot;
169 strength = wrms2 / teff^2;
170 if strength > peakStrength
171 peakStrength = strength;
172 peakWrms = sqrt(wrms2);
173 end
174 gam = min(gmax, P.cascade_gain * Ee ./ (omega * Mtot * teff^2));
175 up = gam .* Ee * h;
176 dn = P.cascade_back * gam .* Ee * h;
177 arrU = accumarray(upIdx, up, [N + 1, 1]);
178 arrD = accumarray(dnIdx, dn, [N, 1]);
179 leak = arrU(N + 1);
180 En = Ee - up - dn + arrU(1:N) + arrD;
181 En(N) = En(N) + (1 - P.cascade_leak) * leak;
182 lost = lost + P.cascade_leak * leak;
184 if sig > 0
185 % fractional energy turnover this step sets the decorrelation
186 moved = up + dn + arrU(1:N) + arrD;
187 rho = exp(-min(30, moved ./ max(Ee, 1e-300)));
188 [u1, st1] = cymbal_rngv(st1);
189 [u2, st2] = cymbal_rngv(st2);
190 zg = sqrt(-2 * log(u1)) .* cos(2 * pi * u2);
191 w = 1 - rho.^2;
192 gl = rho .* gl + sqrt(w) * sig .* zg;
193 gv = rho.^2 .* gv + w * sig^2;
194 end
197end
200cinfo.frames = nf;
201cinfo.stepS = h;
202cinfo.peakStrength = peakStrength;
203cinfo.peakWrmsUm = peakWrms * 1e6;
204cinfo.cascadeLostJ = lost;
205cinfo.injectedJ = sum(sum(0.5 * DV.^2));
206end
↑↓ move↵openescclose