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%
45% What this buys is a cymbal that responds to dynamics the way the instrument
46% does, at the cost of one control-rate loop rather than a full nonlinear
47% integration. What it gives up is phase: energy moves between modes but the
48% intermodulation products of a real cascade, and the split-pair beating of an
49% out-of-round cymbal, are not reproduced. cascade_gain = 0 turns the whole
50% mechanism off and leaves an exactly linear plate.
51%
52% Each contact adds ½ (dv_m)^2 to the energy of every mode, where
53% dv_m = J S(om_m tau) (phi_m(x_c) . dhat)
54% with S(x) = cos(x/2) / (1 - (x/pi)^2) the spectrum of the half-sine contact
55% pulse and dhat the stroke direction (strike_angle_deg away from the surface
56% normal, in the meridian plane). Energies add incoherently, which is right
57% for a plate already ringing and exact for the first stroke into a silent one.
58%
59% Note that only the MAGNITUDE of the excitation survives into the energy
60% state. The signed jumps DV and the frame each contact lands in are returned
61% alongside, because cymbal_radiate.m needs them for the oscillator phases:
62% amplitude comes from the cascade, phase from the linear impulse response.
63%
64% Returns env (N x nf+1, modal VELOCITY amplitude at each frame boundary),
65% tframe (1 x nf+1), DV (N x number of contacts), kf (the frame index of each
66% contact) and cinfo diagnostics.
68N = basis.n;
69omega = basis.omega;
70zeta = basis.zeta;
71h = max(1e-4, P.cascade_step_ms * 1e-3);
72nf = max(1, ceil(duration / h));
73tframe = (0:nf) * h;
75% ---- modal velocity jump of each contact ----
76I = length(tstr);
77DV = zeros(N, I);
78ang = P.strike_angle_deg * pi / 180;
79for i = 1:I
80 [uh, wh] = cymbal_at(basis, sstr(i));
81 proj = -cos(ang) * wh + sin(ang) * uh; % 1 x N, along the stroke
82 x = omega * taustr(i);
83 den = 1 - (x / pi).^2;
84 S = cos(x / 2) ./ den;
85 sing = abs(den) < 1e-6;
86 S(sing) = pi / 4; % limit at x -> pi
87 DV(:, i) = jstr(i) * S .* reshape(proj, N, 1);
88end
89kf = min(nf, max(1, floor(tstr / h) + 1)); % frame each contact lands in
91decay = exp(-2 * zeta .* omega * h);
92gmax = 0.5 / (h * (1 + P.cascade_back));
93Mtot = basis.mass;
94teff = basis.thkMean;
96% ---- the log-frequency ladder ----
97% upIdx(m) is the lowest mode at least cascade_ratio times f(m); modes with no
98% such partner point at N + 1, the sink standing in for the truncated band.
99% dnIdx(m) is the highest mode at most f(m) / cascade_ratio, or m itself at the
100% bottom, which reflects.
101ratio = max(1.02, P.cascade_ratio);
102f = basis.f;
103upIdx = zeros(N, 1);
104dnIdx = zeros(N, 1);
105for m = 1:N
106 ju = find(f > f(m) * ratio);
107 if length(ju) == 0
108 upIdx(m) = N + 1;
109 else
110 upIdx(m) = ju(1);
111 end
112 jd = find(f < f(m) / ratio);
113 if length(jd) == 0
114 dnIdx(m) = m;
115 else
116 dnIdx(m) = jd(length(jd));
117 end
118end
120env = zeros(N, nf + 1);
121Ee = zeros(N, 1);
122peakStrength = 0;
123peakWrms = 0;
124lost = 0;
125for f = 1:nf
126 hit = find(kf == f);
127 for q = 1:length(hit)
128 Ee = Ee + 0.5 * DV(:, hit(q)).^2;
129 end
130 env(:, f) = sqrt(2 * Ee);
132 wrms2 = sum(Ee ./ omega.^2) / Mtot;
133 strength = wrms2 / teff^2;
134 if strength > peakStrength
135 peakStrength = strength;
136 peakWrms = sqrt(wrms2);
137 end
138 gam = min(gmax, P.cascade_gain * Ee ./ (omega * Mtot * teff^2));
139 up = gam .* Ee * h;
140 dn = P.cascade_back * gam .* Ee * h;
141 arrU = accumarray(upIdx, up, [N + 1, 1]);
142 arrD = accumarray(dnIdx, dn, [N, 1]);
143 leak = arrU(N + 1);
144 En = Ee - up - dn + arrU(1:N) + arrD;
145 En(N) = En(N) + (1 - P.cascade_leak) * leak;
146 lost = lost + P.cascade_leak * leak;
147 Ee = En .* decay;
148end
149env(:, nf + 1) = sqrt(2 * Ee);
151cinfo.frames = nf;
152cinfo.stepS = h;
153cinfo.peakStrength = peakStrength;
154cinfo.peakWrmsUm = peakWrms * 1e6;
155cinfo.cascadeLostJ = lost;
156cinfo.injectedJ = sum(sum(0.5 * DV.^2));
157end