How the cymbal is computed

An interactive walk through the theory and the implementation behind this simulator. The physical model is due to Toshiya Samejima (Faculty of Design, Kyushu University):{' '} “Nonlinear physical modeling sound synthesis of cymbals involving dynamics of washers and sticks/mallets,” Acoustical Science and Technology{' '} 42(6), 314–325 (2021) , which itself builds on Bilbao's nonlinear shell synthesis (2010). Equation numbers below, written (1), (20), and so on, are the paper's, so this page can be read beside it. Where the implementation departs from the paper in order to run stably or quickly, Section 10 says exactly how.

{/* ---------------- 1 ---------------- */}

1. What a render computes

Nothing in this app is a recording, an oscillator bank, or a filter tuned to sound like a cymbal. Pressing Render integrates the equations of motion of three coupled mechanical systems, sample by sample, several hundred thousand times: a thin bronze shell (the cymbal itself), the felt washer it is mounted on, and the mallet striking it, the mallet being a flexible stick held by a moving hand with an elastic head on its tip. The sound you hear is the vertical acceleration of one point on the shell's surface, read off the simulation like a contact pickup.

The whole computation is written in MATLAB syntax, in the .m files you can open (and edit) with the app's Code button, and runs in your browser. Each stage lives in its own file:

The first three boxes are precomputation: they build the grid and assemble, then invert, the small matrices the time loop applies at every step. The loop itself, in cymbal2_sim.m, is where the minutes go (or the seconds, when the browser exposes WebGPU and the loop runs on the GPU). The sections below follow this pipeline: first the physics being solved (Sections 2 to 5), then how it is discretized and integrated (Sections 6 to 8), then how the motion becomes audio (Section 9).

{/* ---------------- 2 ---------------- */}

2. Why the model must be nonlinear

A linear model of a struck plate is a fixed set of modes: the strike decides how much energy each mode receives, and from then on each rings at its own constant frequency and decays. Whatever spectrum exists at the moment of impact is the spectrum forever after, only quieter. A real cymbal does something a linear model cannot: struck hard, its sound changes shape after the strike. Pitches glide, and the bright, noise-like wash of a crash builds up over several hundred milliseconds, well after the stick has left. Chaigne, Touzé and Thomas called cymbals “nonlinear percussion instruments” for this reason, and Samejima's paper opens with the same observation.

The mechanism is geometric. A flat plate can bend without stretching, but a plate deflected into a curve, or a shell deflected at all, must stretch its own mid-surface, and the membrane tension this creates acts back on the bending motion. The effect grows as the square of the deflection measured in plate thicknesses: negligible for a whisper-quiet tap, dominant for a crash, where the shell moves by tens of thicknesses. Two small experiments show what this coupling does, and both are audible: each play button integrates the system's equations right here in the page and renders the result. The first is the simplest system that has it, a single oscillator whose spring stiffens with displacement:

The second experiment shows the other, more consequential effect: nonlinearity couples modes to each other, so energy placed low in the spectrum migrates upward. Here is a chain of masses with a quadratic geometric term in its springs (the Fermi–Pasta–Ulam system), started with all of its energy in its lowest mode:

Both effects, the amplitude-dependent detuning and the upward cascade, emerge in the shell model of the next section from a single nonlinear term. Neither is put in by hand. This is the sense in which the crash is computed rather than modeled: hit the simulated shell gently and it rings like a bell; hit it hard and the cascade fills the spectrum, because the equations say so.

{/* ---------------- 3 ---------------- */}

3. The shell and the von Kármán equations

A real cymbal has a raised cup and a gently tapered bow. Following Bilbao (2010), the paper replaces both by one shallow spherical cap of radius{' '} a, uniform thickness h, and a small rise at the center. The rise fixes the radius R_s of the sphere the cap sits on, and that single curvature is all the geometry the model carries:

The price of this simplification is real and audible: the linear spectrum is that of an idealized cap, not of a hammered, lathed instrument, and no irregularity appears anywhere in the model. What it buys is a formulation in which the strong nonlinearity is actually computable. The transverse deflection{' '} u(r, \theta, t) obeys the von Kármán shallow-shell equations (the paper's Eqs. (1) and (2)):

{String.raw`\rho h\,\frac{\partial^2 u}{\partial t^2} = -D\nabla^4 u \;+\; L(\phi, u)\;-\;\frac{1}{R_s}\nabla^2\phi \;-\;2\sigma_0\rho h\,\frac{\partial u}{\partial t} \;+\;2\sigma_1 a^2 \rho h\,\nabla^2\frac{\partial u}{\partial t} \;+\;\epsilon f`} {String.raw`\nabla^4\phi = -\frac{Eh}{2}\,L(u, u)\;+\;\frac{Eh}{R_s}\nabla^2 u`}

Reading Eq. (1) term by term: \rho h is the mass per unit area;{' '} -D\nabla^4 u is ordinary linear plate bending, with{' '} D = Eh^3/12(1-\nu^2) the bending stiffness and \nabla^4 the biharmonic operator (the Laplacian applied twice); the two damping terms carry the fitted loss constants \sigma_0 (frequency-independent) and{' '} \sigma_1 (growing with the squared wavenumber, so high frequencies die faster, as they do on real cymbals); and \epsilon f is the mallet's contact force, applied at the strike point.

The remaining two terms are the shell physics. The function \phi is the{' '} Airy stress function: a single scalar field whose derivatives encode the membrane (in-plane) stress everywhere in the shell. Eq. (2) says where that stress comes from: deflecting the shell stretches it, through the curvature (the{' '} \nabla^2 u / R_s term) and through the quadratic operator{' '} L(u,u). Eq. (1) says what the stress does: it feeds back on the bending motion, again through the curvature and through L(\phi, u). The bracket{' '} L is a fixed bilinear form; in polar coordinates,

{String.raw`L(\phi, u) = \frac{\partial^2\phi}{\partial r^2} \left(\frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial^2 u}{\partial\theta^2}\right) + \frac{\partial^2 u}{\partial r^2} \left(\frac{1}{r}\frac{\partial\phi}{\partial r} + \frac{1}{r^2}\frac{\partial^2\phi}{\partial\theta^2}\right) - 2\left(\frac{1}{r}\frac{\partial^2\phi}{\partial r\,\partial\theta} - \frac{1}{r^2}\frac{\partial\phi}{\partial\theta}\right) \left(\frac{1}{r}\frac{\partial^2 u}{\partial r\,\partial\theta} - \frac{1}{r^2}\frac{\partial u}{\partial\theta}\right)`}

This closed loop, deflection creates stress (2), stress deflects (1), is the engine of everything in Section 2. Because L(u,u) is quadratic in the deflection, the feedback is negligible at small amplitude and dominant at large, and because L mixes radial and azimuthal derivatives, it couples modes of different shape, which is what drives the cascade. The paper works with a dimensionless form of these equations (its Eqs. (3) to (5)); the app integrates the dimensional form above, which is the same content in physical units.

{/* ---------------- 4 ---------------- */}

4. The washer and the free edge

Fourth-order equations need two boundary conditions at each edge. At the outer rim{' '} r = a the shell is free, so the bending moment and the (Kelvin–Kirchhoff) shear force both vanish (the paper's Eqs. (6) and (7)):

{String.raw`\left.\frac{\partial^2 u}{\partial r^2} + \nu\Bigl(\frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial^2 u}{\partial\theta^2}\Bigr)\right|_{r=a} = 0, \qquad \left.\frac{\partial^3 u}{\partial r^3} + \frac{1}{r}\frac{\partial^2 u}{\partial r^2} - \frac{1}{r^2}\frac{\partial u}{\partial r} + \frac{\nu - 3}{r^3}\frac{\partial^2 u}{\partial\theta^2} + \frac{2-\nu}{r^2}\frac{\partial^3 u}{\partial r\,\partial\theta^2}\right|_{r=a} = 0`}

The inner boundary is the paper's first main contribution. Earlier work (Bilbao 2010) clamped the shell rigidly at its center, but a mounted cymbal is nothing like clamped: it sits on a felt washer, and visibly bounces and rocks on it when struck. Samejima models the washer as a mass–damper–spring per unit length of the mounting circle r = c, and imposes it as a boundary condition: the bending moment vanishes there (the shell may tilt freely), and the shear force balances the washer's reaction (Eqs. (14) and (15)):

{String.raw`-D\,[\text{shear expression}]_{\,r=c} \;=\; M_\mathrm{W}\frac{\partial^2 u}{\partial t^2} \;+\; R_\mathrm{W}\frac{\partial u}{\partial t} \;+\; K_\mathrm{W}\,u`}

The paper measures the stiffness of an actual washer,{' '} {String.raw`K_\mathrm{W} = 1.65\times 10^5\ \mathrm{N/m}`} per meter, and sets the mass and damping to zero. The consequences are on two scales. The whole shell acquires rigid-body modes on its mount:

And the audible modes see a compliant inner edge instead of a rigid one, which shifts the low-order resonances; the paper's Fig. 8 shows exactly this between its clamped and impedance-supported calculations, and its Fig. 9 shows the mode shapes that appear only when the support may move, such as the (0,0) mode in which the cymbal moves as one piece. The app defaults to the washer support;{' '} support_clamped restores the rigid clamp of the earlier literature (Eqs. (12), (13)) for comparison.

{/* ---------------- 5 ---------------- */}

5. The mallet: head, stick, grip, hand

The paper's second main contribution is taking the striking mechanism seriously. Most synthesis work injects a prescribed force pulse (a raised cosine, say) at the strike point. Here the force is not prescribed at all: it is whatever comes out of a mallet, itself a dynamical system, colliding with the moving shell.

The head (mass {String.raw`M_\mathrm{M}`}, 28 g for the paper's measured mallet) is an elastic body whose compression \xi against the shell produces a Hertz-type force (Eq. (32)):

{String.raw`f_{\mathrm{M}1} = K_\mathrm{M}\,\bigl[\,\xi\,\bigr]_+^{\,\alpha}, \qquad \xi = u - (u_\mathrm{M} - d_\mathrm{M})`}

where [\cdot]_+ keeps only positive values: the head pushes when compressed and lets go when it separates, and the same force acts back on the shell (Eq. (33)). This one-sided power law is the standard model for hammer–string and mallet–membrane contact (Chaigne & Askenfelt 1994; Rhaouti et al. 1999), which the paper adopts for the cymbal:

The stick is an Euler–Bernoulli beam (Eq. (35)), discretized by the finite element method into cubic beam elements (Eq. (38)), with the head attached at its tip. The grip is the boundary condition at the other end, and it is the paper's third experiment: a tight grip is a clamp (displacement and rotation both follow the hand), a loose grip lets the stick pivot against a rotational spring {String.raw`K_\mathrm{R}`}, measured at 140 N·m/rad. Because the beam is compliant, the tip whips: it arrives faster than the hand that drives it (by a factor near 1.3 at the default stroke), flexural waves reflect along the stick during contact, and the contact force develops the multiple peaks of the paper's Fig. 12. All of this is audible in the attack, and the paper connects it to the measured effect of grip tightness on cymbal tone (Takashima et al. 2018).

The hand closes the system: the grip end's displacement{' '} {String.raw`u_\mathrm{H}(t)`} is prescribed. The paper drives it at constant acceleration to a terminal speed; this app's hand is a small state machine that accelerates, coasts, holds briefly on contact, and retracts to rest, so that repeated strokes and rolls are possible (Section 10 has the details).

{/* ---------------- 6 ---------------- */}

6. Discretizing space: the grid

To compute with the equations of Section 3, the shell's annulus (washer circle to rim) is sampled on a polar grid: N_r intervals in radius,{' '} N_\theta points around the circumference (the paper's Fig. 4, and{' '} cymbal2_grid.m here). Spatial derivatives become central differences: for example

{String.raw`\left.\frac{\partial u}{\partial r}\right|_{l}\approx\frac{u_{l+1}-u_{l-1}}{2h_r}, \qquad \left.\frac{\partial^2 u}{\partial r^2}\right|_{l}\approx\frac{u_{l+1}-2u_l+u_{l-1}}{h_r^2}`}

The biharmonic \nabla^4 reaches two neighbors deep, so evaluating it on the first real line inside a boundary needs values on lines that lie outside{' '} the shell. The grid therefore carries two virtual lines beyond each boundary, and the four boundary conditions of Section 4 (two per edge) are exactly the equations that determine the values there. Nothing is wasted and nothing is free: the boundary conditions are enforced by construction, at the cost of carrying the virtual unknowns in every solve.

The grid also sets what the model can resolve. The paper estimates that with its production grid (24 by 60) the computation is reliable up to roughly the (6,5) mode near 2 kHz; the crash extends far above that, but its top octaves ride on grid modes that are dispersion-warped (Section 7 shows why), which is the main respect in which the simulated crash is duller than a real one. Raising both counts sharpens it, at a cost that grows with their product.

{/* ---------------- 7 ---------------- */}

7. One wavenumber at a time

The implicit time step of the next section needs two elliptic solves per sample: one to find the stress function \phi from Eq. (2), one to advance{' '} u. Done directly, each is a linear system over the whole grid, about 1,500 unknowns at the default resolution, several hundred thousand times per render. The paper assembles exactly these full-grid operators (its Eqs. (18), (19) and Appendix, as Kronecker products). The implementation here exploits a symmetry instead, and this is its one structural departure from the paper (Section 10): every linear operator and boundary row in the model is axisymmetric, unchanged by rotating the cymbal. Such operators cannot mix azimuthal wavenumbers, so decomposing everything into \cos m\theta and \sin m\theta patterns splits each big solve into N_\theta independent radial solves of size{' '} N_r + 2, each of which is precomputed and inverted once before the loop starts. The result of each step is identical (the decomposition is exact, and the self-test verifies it); the cost falls by about two orders of magnitude.

One subtlety is worth naming. On the discrete ring, differentiating twice in{' '} \theta multiplies the wavenumber-m pattern by{' '} {String.raw`\lambda_2(m) = (2\cos m h_\theta - 2)/h_\theta^2`} rather than by -m^2. At low order the two agree; toward the grid cutoff{' '} m = N_\theta/2 the discrete value falls short, so the highest nodal-diameter modes sit flat of where physics would put them. The same happens radially. This numerical dispersion is harmless for the low modes that carry the pitch and wobble, and increasingly wrong for the top of the crash, which is why the tutorial keeps saying the top octaves are approximate.

The nonlinear terms do mix wavenumbers, deliberately: they are the cascade. They are therefore evaluated in physical space, on the grid, and transformed to the wavenumber basis only for the solves. Every time step crosses between the two representations twice.

{/* ---------------- 8 ---------------- */}

8. Marching in time

Time advances in steps of {String.raw`h_t = 1/(f_s \cdot \text{oversample})`}, about 10 μs at the default 96 kHz simulation rate. The scheme, the paper's Eq. (20), splits the physics by stiffness. The linear terms, bending, curvature coupling, damping, are integrated implicitly: the update solves

{String.raw`\boldsymbol{A}\,\boldsymbol{u}_{k+1} = \boldsymbol{B}\,\boldsymbol{u}_k + \boldsymbol{C}\,\boldsymbol{u}_{k-1} \;+\; h_t^2\,[\,\text{nonlinear terms} + \text{contact force}\,]_k`}

where \boldsymbol A, \boldsymbol B, \boldsymbol C contain the biharmonic and Laplacian operators weighted by a free parameter g (the paper uses g = 0.2498), and the boundary rows of Sections 4 replace the equation-of-motion rows on the boundary and virtual lines. The washer's own dynamics, Eq. (15), is folded into two of those rows, so the support costs nothing extra. Why implicit? Because the alternative fails:

The nonlinear terms and the contact force are evaluated explicitly, at the current step, because L(\phi, u) changes every sample and cannot be pre-inverted. This is the standard compromise (Bilbao 2009), and it has a consequence you will meet if you push the hand speed slider: an explicit quadratic term is only conditionally stable, and the condition tightens with amplitude. Strokes whose tip arrives above roughly 10 m/s drive deflections past forty plate thicknesses and the integration diverges, at any rate we tried. This is also about where the shallow-shell expansion itself stops being valid, so the model refuses more or less where its physics runs out; the renderer detects the divergence and reports it rather than emitting garbage. Less intuitively, the azimuthal resolution matters more than the time step: the cascade climbs the nodal-diameter ladder, and if n_theta cuts the ladder short, energy piles up at the cutoff and blows up where a full grid rides it out.

With the operators prebuilt, one pass through the loop in{' '} cymbal2_sim.m is:

Contact. Compression \xi from the shell and beam-tip positions; force {String.raw`f = K_\mathrm{M}[\xi]_+^\alpha`} if touching, zero if not.
Mallet. Advance the stick beam one step (its own implicit two-step scheme), driven by the hand at one end and the contact force at the tip.
Hand. Advance the stroke state machine (accelerate, coast, hold, retract).
Stress function. Evaluate L(u,u) and \nabla^2 u on the grid, transform to wavenumbers, solve Eq. (2) for \phi (one small radial solve per wavenumber), transform back.
Nonlinear forcing. Evaluate L(\phi, u) and the membrane term on the grid; inject the contact force at the strike node.
Shell. Transform the total forcing, apply the prebuilt{' '} {String.raw`\boldsymbol A^{-1}\boldsymbol B`},{' '} {String.raw`\boldsymbol A^{-1}\boldsymbol C`},{' '} {String.raw`\boldsymbol A^{-1}`} per wavenumber, transform back:{' '} {String.raw`\boldsymbol u_{k+1}`}.
Pickup. Record the second time difference of u at the pickup node: the acceleration sample.

A full render at the defaults is about 580,000 such steps. In the browser the same loop runs either through the numbl interpreter or, when WebGPU is available, as a sequence of GPU compute passes per step, which is roughly twenty times faster; the diagnostics panel reports which path a render took.

{/* ---------------- 9 ---------------- */}

9. From motion to sound

The rendered signal is the shell's normal acceleration at the pickup point, just as in the paper's numerical experiments (its output point is{' '} \bar r = 0.6 on the meridian opposite the strike, and that is the app's default too). cymbal2_output.m conditions it: a gentle DC blocker removes the quasi-static offset the washer relaxation leaves; a windowed-sinc filter anti-aliases and decimates from the simulation rate down to the output rate; a short fade closes the tail; and the result is normalized to the requested peak level.

It is worth being precise about what this signal is and is not. It is the motion of one point of the metal: a contact pickup, not a microphone. Radiating the sound into air, with the interference and directivity of a real plate source, would require the full surface velocity history and an acoustic model on top, and is not part of the paper's method. The audible consequence is that low-order modes which radiate poorly from a real cymbal are present here at full strength, so renders carry somewhat more low end than a room recording of the same instrument would. With that caveat, the paper's trial calculations, and this implementation against the demonstration recording published by the paper's authors, reproduce the features that define the instrument: the clear low-order peaks, the slow buildup of high-frequency energy after a hard stroke, the difference between a washer-mounted and a clamped cymbal, and the effect of the grip on the attack. The spectrogram panel of the simulator is the place to watch the first of these; render the Paper demo preset and look for the energy arriving above 5 kHz over the first half second.

{/* ---------------- 10 ---------------- */}

10. Where the implementation departs from the paper

The physics above is Samejima's throughout. Turning a published scheme into code that must run stably for arbitrary knob settings forced a small number of departures, and honesty requires naming them; the repository README states each with its full justification.

The head is lumped onto the beam tip. The paper integrates the head as its own one-degree body (Eq. (34)) coupled to the beam through a rigid connection (Eqs. (36), (37)). Coupling that explicit update to the implicit beam turns out to be unstable as a discrete pair. Since the connection is rigid, substituting it into Newton's law puts the head mass on the tip node exactly, which is the same physics and inherits the beam scheme's stability.
The implicit weighting of the washer and beam schemes. The paper prints the split with weight (1-\phi)/2 on the old and new steps; that choice is unstable whenever the mass term does not dominate (the massless washer included). Moving the implicit weight \phi/2 onto the average of the new and old steps makes \phi = 0.7 unconditionally stable, and that is what runs here.
The hand can retract. The paper's hand accelerates to a terminal speed and never reverses, which is fine for one stroke and cannot play a second. The app's hand accelerates at an arm-like rate, coasts, and snaps back after contact, so figures and rolls are possible.
The solves are block-diagonalized over wavenumber (Section 7). This changes the cost, not the result.
Guard rails. The loop detects divergence and stops with a message rather than emitting garbage, and the system builder truncates azimuthal columns whose inner boundary layer the radial grid cannot resolve (none, at the default grid).

{/* ---------------- 11 ---------------- */}

11. Credits and references

The physical model implemented by this app, the washer boundary condition, the stick/mallet system, and the numerical scheme, are the work of{' '} Toshiya Samejima (Faculty of Design, Kyushu University), published as the paper below. The defaults are the paper's measured instrument: an 18-inch medium crash on its measured washer, struck by its measured mallet. This tutorial and the implementation are the app's own; errors in either should be laid at the app's door, not the paper's.

T. Samejima, “Nonlinear physical modeling sound synthesis of cymbals involving dynamics of washers and sticks/mallets,” Acoustical Science and Technology 42(6), 314–325 (2021).{' '} doi:10.1250/ast.42.314
S. Bilbao, “Percussion synthesis based on models of nonlinear shell vibration,” IEEE Transactions on Audio, Speech, and Language Processing 18, 872–880 (2010).
S. Bilbao, Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics (Wiley, 2009).
A. Chaigne, C. Touzé and O. Thomas, “Nonlinear vibrations and chaos in gongs and cymbals,” Acoustical Science and Technology{' '} 26, 403–409 (2005).
A. Chaigne and A. Askenfelt, “Numerical simulations of piano strings. I. A physical model for a struck string using finite difference methods,”{' '} Journal of the Acoustical Society of America 95, 1112–1118 (1994).
L. Rhaouti, A. Chaigne and P. Joly, “Time-domain modeling and numerical simulation of a kettledrum,” Journal of the Acoustical Society of America 105, 3545–3562 (1999).
S. Takashima, N. Wakatsuki and K. Mizutani, “Measurement of tone change of cymbal depending on grip tightness,” Proceedings of Musical Acoustics, Acoustical Society of Japan 37(6), 71–74 (2018, in Japanese).