Cymbal#

A physically-based simulation of a struck cymbal, sampled as acoustic pressure at a point in space. The instrument is a bronze shell of revolution with a raised bell and a tapered bow; a Hertz stick contact drives its bending and stretching modes, computed by a one-dimensional finite element eigensolve per nodal-diameter order, and the geometric nonlinearity that turns a hard stroke into a crash is carried as a cascade of energy up the modal frequency ladder. The hammering that makes a real cymbal irregular is carried too, because without it the model is a bell. The plate radiates to a listener point as a sheet of normal dipoles.

There is no microphone model: no directivity, no capsule response, no gain staging. The output is p(r, t) in pascals.

Contents#

The MATLAB code runs from the command line with cd matlab && numbl run cymbal_demo.m, or cymbal_validate.m to run the checks. The web app is a Vite project; see webapp/README.md.

The model#

The signal chain, with each stage in its own .m file:

meridian profile: bell, bow, thickness taper       cymbal_profile.m
  -> thin-shell elements for one Fourier order     cymbal_shell.m
      -> one eigensolve per nodal-diameter order   cymbal_basis.m
         + hammering: hybridize near neighbours
         + radiated power -> radiation damping
stroke gesture -> Hertz stick contacts             cymbal_strikes.m
  -> modal energies, their nonlinear flow up the
     frequency ladder, and their fluctuation       cymbal_cascade.m
      -> oscillator bank + dipole-sheet transfer   cymbal_radiate.m
          -> free-field pressure at the listener
              -> anti-aliased decimation -> PCM    cymbal_output.m

Four points are worth stating plainly, because they are where the model commits itself.

The modes come from a shell of revolution, one Fourier order at a time. A cymbal is axisymmetric, so its displacement field separates exactly in the azimuth: every mode carries n nodal diameters and varies as cos(n theta), and the two-dimensional eigenproblem breaks into one one-dimensional problem per n. Each is discretized along the meridian with two-node Love-Kirchhoff shell elements, linear in the meridional and circumferential displacements and cubic Hermite in the normal one, integrated by Gauss quadrature against the tabulated profile rather than frozen into conical frusta (cymbal_shell.m). Rotary inertia is kept, which shifts the retained frequencies by well under a percent but is what makes the mass matrix positive definite, so each order reduces through the Cholesky factor of the mass matrix and is diagonalized by svd. The curvature of the bow is what earns this formulation its keep: on a flat disc, bending and stretching are independent, and it is the meridian curvature that couples them and pulls the partials of a cymbal away from the Chladni frequencies of a plate. Shear deformation is neglected, and the mounting is taken as free-free with a point dashpot at the hole for the felt washer, which reaches the n = 0 and n = 1 modes and almost nothing else. Only the cos member of each degenerate pair is retained: that is exact for the model as posed, since every contact sits on the meridian theta = 0 where the sin member has a node.

A perfectly axisymmetric plate is a bad cymbal, and the model has to say so. This is the point at which the idealization fails audibly rather than quietly. On an exactly round plate the modes are exact cos(n theta) shapes, and such a mode radiates only into the cone where its azimuthal trace speed is supersonic, n < k a sin(theta). Below a few kilohertz that leaves only the handful of n = 0, 1, 2 modes reaching any listening position at all: of the several hundred modes the eigensolve retains, six carry half the radiated energy and three carry half of the tail, and the result is unmistakably a small bell. Real cymbals are hammered and lathed, and that is not decoration. The irregularity destroys the exact azimuthal separation, so every mode of the real plate is a hybrid of the ideal ones lying near it in frequency, and a hybrid carries low-order content and therefore radiates. The model applies that as first-order degenerate perturbation theory: a seeded orthogonal mixing of modes within hammer_bw of each other in frequency, built from Givens rotations on neighbouring pairs so that it conserves energy exactly, plus a seeded relative detune of hammer_detune which is what makes near neighbours beat. Because the window is relative, the sparse low modes stay pure and only the dense region above a kilohertz hybridizes, which is where the wash lives. hammer_mix = 0 recovers the ideal plate, and is worth listening to once. What this does not do is model the hammering pattern itself; the mixing is statistical, so a specific instrument cannot be reproduced from a photograph of its bow.

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.

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.

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.

Validation#

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.

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.

Defaults#

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.

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.

Setting strike_count to 8 gives a repeated figure alternating between the edge and the bow; moving strike_r to 0.6 and dropping the velocity gives a ride, and to 0.1 the bell. Dropping strike_velocity to a few tenths of a metre per second turns the same instrument back into a ping, because the cascade has nothing left to work with.