/ simulated_instruments / cymbal
Sign in
simulated_instruments / cymbal
Default the turbulent fluctuation off
It broadens the partials as intended but its decorrelation rate runs up to the control rate, so what it adds is modulation grit rather than a wash. The parameter and the code stay, since the physics behind it is not what is wrong; the validation now switches it on explicitly where it checks that it freezes in the linear limit.
Jeremy Magland <jmagland@flatironinstitute.org> committed commit 3a32a66f0dcc parent 7da0671 Browse files
7 changed files+46−23
README.mdmodified+3−3View file
@@ -39,19 +39,19 @@ Four points are worth stating plainly, because they are where the model commits
3939
4040 **The nonlinearity is carried as energy, not as a coupled integration.** A cymbal struck hard is not a linear resonator. Its transverse displacement is comparable to the plate thickness, the von Karman stretching of the mid-surface couples the modes, and energy injected low by the stick migrates upward over tens of milliseconds, which is what the ear hears as a crash opening into shimmer after its attack. Integrating that directly means a dense cubic coupling tensor over hundreds of modes, so what is integrated instead is the energy of each mode (`cymbal_cascade.m`): each mode has an up-target a fixed frequency ratio above it and a down-target the same ratio below, and energy performs a biased random walk between them at a rate `cascade_gain * om * (w/t)^2`, proportional to the mode's own squared amplitude-to-thickness ratio because that is the order of the geometric nonlinearity. Steps are taken in log frequency rather than in mode index on purpose, so the swell does not depend on how densely the band happens to be resolved. Because the flux out of a mode goes as the square of its energy, a nearly empty mode passes nothing on, the front of the cascade advances rather than the whole band draining at once, and the swell stops on its own. What this buys is the instrument's central behaviour: the same cymbal pings when tapped and crashes when hit, with no velocity layers anywhere in the model. What it gives up is phase. Amplitude comes from the cascade and phase from a linear impulse response carried alongside it, so a stroke re-phases the modes it finds quiet and a coherent attack survives, but the intermodulation products of a real cascade do not. `cascade_gain` is the one fitted number in the model; setting it to zero leaves an exactly linear plate, which sounds like a small gong.
4141
42-What those equations track is the *mean* energy of each mode, and the mean is not what a listener hears. A mode sitting in a turbulent field is continually being filled and drained by its neighbours, and its amplitude wanders about that mean on the timescale of the exchange. This matters more than it sounds. A few hundred modes with smooth envelopes are a few hundred audible partials, spaced tens of hertz apart with bandwidths of a couple of hertz, and that is a chime however the energy is distributed among them; give each one a wandering amplitude with a correlation time of a few milliseconds and each partial broadens by a couple of hundred hertz, the modal overlap passes one, and the comb closes into the continuum a wash actually is. Each mode therefore carries a log-amplitude following an Ornstein-Uhlenbeck process at the control rate, with the correlation rate taken as that mode's own fractional energy turnover rather than as a fitted time, so a mode being fed while nearly empty decorrelates every step and a mode merely ringing down does not decorrelate at all. A running variance keeps the exponential unbiased, so the fluctuation redistributes energy in time without creating or destroying any, and it switches itself off exactly where the cascade does. `turbulence = 0` removes it and leaves the mean-energy solve, which chimes.
42+What those equations track is the *mean* energy of each mode, and the mean is not quite what a listener hears: a mode sitting in a turbulent field is continually being filled and drained by its neighbours, and its amplitude wanders about that mean on the timescale of the exchange. `cymbal_cascade.m` carries that wander as an Ornstein-Uhlenbeck process on each mode's log-amplitude, with the correlation rate taken as that mode's own fractional energy turnover rather than as a fitted time, and a running variance keeping the exponential unbiased so that energy is redistributed in time without being created or destroyed. It is switched off at the defaults, and that is worth being straightforward about. On paper it does what it should, broadening each partial past modal overlap and closing the comb into a continuum; in practice the turnover rate it derives runs up to the control rate itself, so what reaches the ear is closer to modulation grit than to a wash. The parameter is `turbulence`, it is kept because the physics behind it is not the thing that is wrong, and a few tenths is the gentle end of it. With it at zero the wash is carried by the hammering and by modal density alone, and both of those could stand to be larger.
4343
4444 **Radiation is a dipole sheet, computed two ways because two regimes are needed.** Each surface element carries a dipole along the local normal of volume `dS * ldip`, where `ldip` is half the shorter of the flexural and acoustic wavelengths, capped at the plate radius: the distance over which the air short-circuits from one face to the other. That length is the one phenomenological quantity in the radiation model, and below the plate's coincidence frequency, about 18 kHz for 1 mm bronze, it is the flexural wavelength that sets it, which is why a cymbal is such a poor radiator in the low register and improves steadily as the cascade fills the top. For the listener, who at half a metre from a 40 cm plate is nowhere near the far field, the element dipoles are summed coherently over the azimuth with their true distances and the near-field term kept (`cymbal_radiate.m`). For the radiated power, which sets the radiation damping and is defined only in the far field, the azimuthal integral of the same kernel collapses onto Bessel functions and the power integral reduces to one quadrature in `cos(theta)` (`cymbal_basis.m`); the two agree to a few parts in a thousand, which the validation checks. Summing element powers incoherently, as the sibling `triangle` does for a rod, would not do here: for a continuous surface that estimate halves when the elements do. The multipole structure comes out of the azimuthal sum on its own. A mode with `n` nodal diameters is a `2^n`-pole, so the high-order modes are weak until their wavelength is short, and on the symmetry axis every one of them cancels exactly.
4545
4646 ## Validation
4747
48-`matlab/cymbal_validate.m` (also `npm run selftest` in `webapp/`) checks the shell elements in their flat limit against the closed-form free circular plate, whose frequency equation is the vanishing of a two-by-two determinant in `J_n` and `I_n` under zero radial moment and zero Kelvin-Kirchhoff shear, together with the rigid-body count each Fourier order must have (two at `n = 0`, two at `n = 1`, none above); the mass-orthonormality and stiffness-diagonalization identities of the eigensolve; the exact orthogonality of the hammering mixing, which is what makes it energy conserving; the Bessel far-field reduction against brute-force summation over a finely discretized dipole sheet; the exact cancellation of every mode with a nodal diameter on the symmetry axis, with the hammering switched off, since hybridizing the modes is precisely what breaks that cancellation; the Hertz `v^(-1/5)` contact-time scaling in the light-stick limit; exact energy conservation of the cascade scheme with the damping off, the top of the ladder reflecting and the fluctuation suppressed, and its collapse onto closed-form damped exponentials when `cascade_gain` is zero, with the fluctuation left switched on so that its freezing in the linear limit is checked too; the far-field `-6 dB` per doubling of distance; exact reproducibility of the seeded stroke pattern; and pinned end-to-end diagnostics.
48+`matlab/cymbal_validate.m` (also `npm run selftest` in `webapp/`) checks the shell elements in their flat limit against the closed-form free circular plate, whose frequency equation is the vanishing of a two-by-two determinant in `J_n` and `I_n` under zero radial moment and zero Kelvin-Kirchhoff shear, together with the rigid-body count each Fourier order must have (two at `n = 0`, two at `n = 1`, none above); the mass-orthonormality and stiffness-diagonalization identities of the eigensolve; the exact orthogonality of the hammering mixing, which is what makes it energy conserving; the Bessel far-field reduction against brute-force summation over a finely discretized dipole sheet; the exact cancellation of every mode with a nodal diameter on the symmetry axis, with the hammering switched off, since hybridizing the modes is precisely what breaks that cancellation; the Hertz `v^(-1/5)` contact-time scaling in the light-stick limit; exact energy conservation of the cascade scheme with the damping off, the top of the ladder reflecting and the fluctuation suppressed, and its collapse onto closed-form damped exponentials when `cascade_gain` is zero, with the fluctuation switched explicitly on so that its freezing in the linear limit is checked too; the far-field `-6 dB` per doubling of distance; exact reproducibility of the seeded stroke pattern; and pinned end-to-end diagnostics.
4949
5050 As in the sibling `maraca`, `triangle`, `recorder` and `bass` repositories, there is no separate Python reference implementation: the MATLAB code is the reference, and the pinned numbers were recorded from it when the model was built.
5151
5252 ## Defaults
5353
54-A 16-inch B20 bronze crash-ride: 20.3 cm radius, a 4.5 cm bell rising 1.7 cm above a bow that itself rises 2 cm from the rim, and a thickness tapering from 1.15 mm at the bell to 0.75 mm at the edge, which comes to 1.05 kg. Solving to 56 nodal diameters retains 378 modes between 49 Hz and 10.0 kHz, each a hybrid of about five ideal ones. The default gesture is a single crash on the edge at 3 m/s with a 30 g nylon-tipped stick held 40 degrees off the surface normal: an 88 microsecond Hertz contact, a 291 N peak force, a peak plate displacement of 149 micrometres or 0.16 thicknesses, and a peak of 20.5 Pa (120 dB SPL) at the 56 cm listener before normalization. About 6% of the struck energy leaves the top of the retained band, which is the shimmer above 10 kHz the model does not carry; raising `max_mode_hz` and `n_diameters` together recovers it at the cost of the eigensolves. A full render takes about six seconds, of which the eigensolves are three.
54+A 16-inch B20 bronze crash-ride: 20.3 cm radius, a 4.5 cm bell rising 1.7 cm above a bow that itself rises 2 cm from the rim, and a thickness tapering from 1.15 mm at the bell to 0.75 mm at the edge, which comes to 1.05 kg. Solving to 56 nodal diameters retains 378 modes between 49 Hz and 10.0 kHz, each a hybrid of about five ideal ones. The default gesture is a single crash on the edge at 3 m/s with a 30 g nylon-tipped stick held 40 degrees off the surface normal: an 88 microsecond Hertz contact, a 291 N peak force, a peak plate displacement of 149 micrometres or 0.16 thicknesses, and a peak of 12.9 Pa (116 dB SPL) at the 56 cm listener before normalization. About 6% of the struck energy leaves the top of the retained band, which is the shimmer above 10 kHz the model does not carry; raising `max_mode_hz` and `n_diameters` together recovers it at the cost of the eigensolves. A full render takes about six seconds, of which the eigensolves are three.
5555
5656 The damping is deliberately flat across the band, about 4 seconds of `T60` at 300 Hz falling to 2 seconds at 6 kHz. That is a choice, and an audible one: a steeper law, with the low modes ringing for the fifteen seconds a frequency-independent loss factor would give them, leaves a handful of low partials sounding alone once the wash has gone, which is again a bell. `damp_rate` is the term that holds the whole band together, and `loss_factor` and `loss_slope` are what tilt it.
5757
matlab/cymbal_cascade.mmodified+14−12View file
@@ -42,18 +42,21 @@ function [env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, tau
4242 % the model: it is a dimensionless coupling constant, set so that a default
4343 % crash takes about the 40 ms a real one takes to open up.
4444 %
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.
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.
5558 %
56-% So each mode carries a log-amplitude that follows an Ornstein-Uhlenbeck
59+% Each mode carries a log-amplitude that follows an Ornstein-Uhlenbeck
5760 % process discretized at the control rate,
5861 % rho_m = exp(-(energy exchanged this step) / E_m)
5962 % g_m <- rho_m g_m + sqrt(1 - rho_m^2) turbulence z_m
@@ -65,7 +68,6 @@ function [env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, tau
6568 % and the process is seeded, so a render is reproducible. Where the cascade is
6669 % off the turnover is zero, rho is one, g stays where it started, and the
6770 % envelopes are exactly the closed-form exponentials of a linear plate.
68-% turbulence = 0 removes the fluctuation and leaves the mean-energy solve.
6971 %
7072 % What this buys is a cymbal that responds to dynamics the way the instrument
7173 % does, at the cost of one control-rate loop rather than a full nonlinear
matlab/cymbal_defaults.mmodified+1−1View file
@@ -60,7 +60,7 @@ P.cascade_back = 0.45; % downward / upward transfer ratio
6060 P.cascade_ratio = 1.19; % frequency ratio of one transfer step (2^0.25)
6161 P.cascade_leak = 1.0; % 1 = energy leaving the top mode is lost
6262 P.cascade_step_ms = 1.0;
63-P.turbulence = 0.9; % log-amplitude spread of the fluctuation
63+P.turbulence = 0.0; % log-amplitude spread of the fluctuation (off)
6464
6565 % Radiation / listener
6666 P.direct_gain = 1.0;
matlab/cymbal_validate.mmodified+3−3View file
@@ -139,7 +139,7 @@ bz.zeta = zeros(bz.n, 1);
139139 Etot = sum(envC.^2 / 2, 1);
140140 eCons = max(abs(Etot / Etot(1) - 1));
141141 Q.cascade_gain = 0;
142-Q.turbulence = P.turbulence;
142+Q.turbulence = 0.9; % explicitly on, so the freeze is actually tested
143143 [envL, tfrL] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.4);
144144 pred = envL(:, 1) * exp(-(b6.zeta .* b6.omega).' * 0 + 0);
145145 pred = (envL(:, 1) * ones(1, length(tfrL))) .* exp(-(b6.zeta .* b6.omega) * tfrL);
@@ -181,12 +181,12 @@ fprintf('ride (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
181181
182182 % ---- 9. pinned end-to-end diagnostics at the defaults ----
183183 [y, basis, dinfo] = cymbal_render(P);
184-ok = abs(dinfo.peakPa - 20.5278) < 0.02 ...
184+ok = abs(dinfo.peakPa - 12.9002) < 0.02 ...
185185 && abs(basis.f(1) - 49.4130) < 0.02 ...
186186 && basis.n == 378;
187187 fails = fails + (~ok);
188188 fprintf('pinned : peak %.4f Pa (target %.4f), mode1 %.3f Hz (target %.3f), %d modes %s\n', ...
189- dinfo.peakPa, 20.5278, basis.f(1), 49.4130, basis.n, cv_verdict(ok));
189+ dinfo.peakPa, 12.9002, basis.f(1), 49.4130, basis.n, cv_verdict(ok));
190190
191191 % ---- 10. the wash decays ----
192192 % A cymbal rings for seconds, so the tail is not silent, but the energy in the
matlab/d_now.madded+21−0View file
@@ -0,0 +1,21 @@
1+P = cymbal_defaults(); Q = P; Q.normalize = 0; Q.full_scale_pa = 1;
2+fs = P.sample_rate; win = 4096;
3+for turb = [0 0.9]
4+ R = Q; R.turbulence = turb;
5+ [y, b, d] = cymbal_render(R);
6+ fprintf('turbulence %.1f peak %6.2f Pa flat@ ', turb, d.peakPa);
7+ for tt = [0.01 0.10 0.40 1.00 2.00]
8+ i0 = round(tt*fs)+1;
9+ x = y(i0:i0+win-1) .* (0.5*(1-cos(2*pi*(0:win-1)/(win-1))));
10+ Y = abs(fft(x)).^2; Y = Y(1:win/2); fq = (0:win/2-1)*fs/win;
11+ Ys = Y(fq>=1000 & fq<8000) + 1e-24;
12+ fprintf('%.3f ', exp(mean(log(Ys)))/mean(Ys));
13+ end
14+ i0 = round(1.0*fs)+1;
15+ x = y(i0:i0+win-1) .* (0.5*(1-cos(2*pi*(0:win-1)/(win-1))));
16+ Y = abs(fft(x)).^2; Y = Y(1:win/2); fq = (0:win/2-1)*fs/win;
17+ fprintf(' | 1s low-vs-high %.0f dB\n', 10*log10(sum(Y(fq<600))+1e-30) - 10*log10(sum(Y(fq>=1500&fq<8000))+1e-30));
18+end
19+P2 = cymbal_defaults();
20+[y2, b2, d2] = cymbal_render(P2);
21+fprintf('\npinned: peakPa %.6f f1 %.6f n %d\n', d2.peakPa, b2.f(1), b2.n);
webapp/README.mdmodified+2−2View file
@@ -2,7 +2,7 @@
22
33 Interactive browser implementation of the physical model described in `../README.md`. The algorithm itself is written in MATLAB syntax, in `../matlab/*.m`, and executed client-side by [numbl](https://numbl.org) inside a Web Worker; the TypeScript here is only the UI, the plots, the WAV encoder, and the marshalling of parameters and tensors in and out of the numbl runtime.
44
5-The signal chain: a Hertz stick contact (or a seeded repeated figure, with `Strokes` above 1), the modes of a bronze shell of revolution from one finite element eigensolve per nodal-diameter order, hybridized to stand for the hammering that makes a real cymbal irregular, a nonlinear cascade of modal energy up the frequency ladder with the turbulent fluctuation about it, an oscillator bank radiating through a coherent sum of surface dipoles to the listener point, and Kaiser-windowed anti-aliased decimation to 48 kHz.
5+The signal chain: a Hertz stick contact (or a seeded repeated figure, with `Strokes` above 1), the modes of a bronze shell of revolution from one finite element eigensolve per nodal-diameter order, hybridized to stand for the hammering that makes a real cymbal irregular, a nonlinear cascade of modal energy up the frequency ladder, an oscillator bank radiating through a coherent sum of surface dipoles to the listener point, and Kaiser-windowed anti-aliased decimation to 48 kHz.
66
77 ## Usage
88
@@ -47,7 +47,7 @@ The plot the sibling repositories do not have. The cascade solve tracks the ener
4747
4848 ## Three sliders worth pulling first
4949
50-`Hybridization` at 0 is the ideal, perfectly round plate. Its high nodal-diameter modes cannot radiate at all, so a few low-order partials carry everything and the instrument becomes a small bell; sliding it back up is hammering the cymbal. `Fluctuation` at 0 leaves every mode with a smooth envelope, which is a few hundred clean partials, that is to say a chime; up at 0.9 each partial broadens into a band and they close into a wash. `Stroke velocity` sweeps the whole dynamic behaviour through one physical quantity, because the cascade strength goes as the square of the amplitude: at 0.7 m/s the same cymbal pings, at 6 m/s it opens right up.
50+`Hybridization` at 0 is the ideal, perfectly round plate. Its high nodal-diameter modes cannot radiate at all, so a few low-order partials carry everything and the instrument becomes a small bell; sliding it back up is hammering the cymbal. `Fluctuation` is off at the defaults: it broadens each partial by letting its amplitude wander about the mean the cascade solves for, which on paper closes the comb into a wash and in practice adds grit, so a few tenths is as far as it is worth taking. `Stroke velocity` sweeps the whole dynamic behaviour through one physical quantity, because the cascade strength goes as the square of the amplitude: at 0.7 m/s the same cymbal pings, at 6 m/s it opens right up.
5151
5252 ## A note on the seeded generator
5353
webapp/src/params.tsmodified+2−2View file
@@ -124,7 +124,7 @@ export const DEFAULTS: CymbalParams = {
124124 cascade_ratio: 1.19,
125125 cascade_leak: 1.0,
126126 cascade_step_ms: 1.0,
127- turbulence: 0.9,
127+ turbulence: 0.0,
128128 direct_gain: 1.0,
129129 dipole_scale: 1.0,
130130 n_azimuth: 224,
@@ -206,7 +206,7 @@ export const PARAM_GROUPS: ParamGroup[] = [
206206 { key: 'cascade_gain', label: 'Cascade gain', min: 0, max: 4000, step: 10, help: 'Coupling constant of the energy cascade. 0 leaves an exactly linear plate, which sounds like a gong; raising it opens the crash up. This is the one fitted number in the model.' },
207207 { key: 'cascade_ratio', label: 'Step ratio', min: 1.03, max: 2.0, step: 0.01, help: 'Frequency ratio of one transfer step. Smaller means a slower, smoother climb; 1.19 is a quarter of an octave.' },
208208 { key: 'cascade_back', label: 'Backward share', min: 0, max: 0.95, step: 0.01, help: 'Downward transfers per upward one. At 1 the flux stops and energy only diffuses.' },
209- { key: 'turbulence', label: 'Fluctuation', min: 0, max: 2.0, step: 0.05, help: 'Log-amplitude spread of the turbulent wander about the mean modal energy. This is what broadens each partial into a band and turns a few hundred audible tones into a wash; at 0 the modes hold smooth envelopes and the cymbal chimes.' },
209+ { key: 'turbulence', label: 'Fluctuation', min: 0, max: 2.0, step: 0.05, help: 'Log-amplitude wander about the mean modal energy, off by default. On paper it broadens each partial and closes the comb into a wash; in practice it decorrelates too fast to sound like one, and adds grit. A few tenths is the gentle end.' },
210210 { key: 'cascade_leak', label: 'Top-of-band leak', min: 0, max: 1, step: 0.05, help: 'Fraction of the flux leaving the highest retained mode that is lost, standing in for the truncated band above. 0 reflects it back.' },
211211 ],
212212 },
↑↓ move↵openescclose