1# bead-percussion
3A browser app for idealized percussion: a network of beads (point masses) joined by springs, set up by a MATLAB-syntax script that runs in the browser with [numbl](https://numbl.org), integrated accurately in time, and heard through the motion (by default the acceleration) of a listener bead.
5The aim is to reproduce the phenomena of drums and cymbals (resonant modes, nonlinear coupling between them, transfer of energy to higher frequencies, amplitude-dependent detuning) in a system that is physical but much simpler than the continuum equations of a membrane or a shell. Nothing here discretizes a PDE. The bead network is the system.
7## Running
9```bash
10npm install
11npm run dev # the app, at the printed localhost URL
12npm test # solver tests (Node; the WebGPU path is checked in the browser with src/devcompare.ts)
13```
15## The model
17Each bead has a mass and a position in 3D. The script defines:
19- *Springs* between pairs of beads, with stiffness k, rest length L0 and optional cubic stiffness k3: the force along the spring is k e + k3 e^3, with e the current length minus L0. A spring whose rest length is shorter than the distance between its beads pulls with a force at rest, and this pre-tension gives a flat net its transverse stiffness, as tension does for a drum head. Because the springs are real springs in 3D, the geometric nonlinearity (a transverse displacement stretches the springs) is exact rather than a truncated expansion.
20- *Bending triplets* on three beads in a row, i, j, k, with energy kb/2 |x_i - 2 x_j + x_k|^2: the middle bead resists leaving the straight line through its neighbours. The force is linear and the energy is unchanged by translations and rotations. This gives plates and bars their stiffness without tension, so their edges can be free. The rest shape must be straight and evenly spaced along each triplet (flat plates and straight bars; a domed shape will need a different element).
21- *Hinges* across the shared edge j-k of two triangles (i, j, k) and (l, k, j), with energy kh/2 (theta - theta0)^2, where theta is the angle between the triangles and theta0 its value at the rest positions. Any rest shape is therefore free of stress, which is what a dome or a bell needs; `bp_hinges(tri)` lists the hinges of a triangle mesh. The gradient of theta is that of Bridson, Marino and Fedkiw (2003). At rest the stiffness is exactly kh grad(theta) grad(theta)^T, so the time step and the linear modes remain analytic.
22- *Dashpots* between pairs of beads, with force c (v_j - v_i) on bead i, plus an optional axial part ca (d . (v_j - v_i)) d along the current direction d between the beads. Short-wavelength motion has larger relative velocities between neighbours, so this damping grows with frequency.
23- *Friction* on each bead, force -c0 v, which damps all frequencies alike.
24- *Pinned* beads, held fixed.
25- *Hits*: an instant velocity given to a set of beads, or a raised-cosine force pulse of given peak and duration. Any number, at any times.
26- A *listener*: a bead (or the average of several), a direction, and a quantity: the acceleration (the default), velocity or displacement along that direction is the audio signal. Acceleration is the closest of the three to what a small vibrating surface radiates into the air; displacement weights each frequency by 1/f^2 relative to it, so it sounds dull and is dominated by the lowest modes. The acceleration is computed from the forces on the listener bead at each sample (springs, bending, force pulses and damping), not by differencing positions.
28Damping is easiest to set by its effect. `sys = bp_decay(sys, [f1 T1; f2 T2])`, called after the elements are set, chooses the friction, the dashpots and the damping of the bending triplets and hinges so that every transverse mode of a flat net falls 60 dB in time T(f) following the law sigma(f) = sigma0 + sigma1 (2 pi f)^2 for the amplitude decay rate, through the two given points (one point gives the same decay time at every frequency). It does this by giving each bead friction 2 sigma0 m; each spring a parallel dashpot whose isotropic part is 2 sigma1 times the spring's tension per length g and whose axial part is 2 sigma1 (keff - g) (a spring's stiffness at rest is (keff - g) d d^T + g I); and each triplet and hinge a damping 2 sigma1 times its stiffness. The damping matrix at rest is then 2 sigma0 M + 2 sigma1 K for all the motion, so the modes stay uncoupled and the law holds exactly at small amplitude, for flat and curved nets and for in-plane as well as transverse modes. At large amplitude the same elements damp the motion, but not exactly by this law. The page plots the small-amplitude decay time of every mode for whatever damping the script sets. A script that sets neither `sys.friction` nor `sys.dashpots` gets a default under which every mode falls 60 dB in 2.5 s; `sys.friction = 0` turns damping off.
30A free object needs support somewhere, or it can drift and spin: pin a small cluster of beads (a single pinned bead lets a plate tilt about it), or hold beads softly with `bp_anchor`, which ties each to four fixed anchors in its own plane with pre-tensioned springs of a chosen out-of-plane stiffness, like a cord or a felt. The run warns when the linear analysis finds modes at 0 Hz.
32The `matlab/helpers/bp_*.m` functions build common geometries (`bp_hexdisk`: a triangular lattice clipped to a disk, with its rim beads and bending triplets; `bp_chain`: a line of beads, with its triplets), add soft supports (`bp_anchor`), list edges (`bp_edges`), make springs with a given tension at rest (`bp_springs`), list hinges of a mesh (`bp_hinges`), set damping from decay times (`bp_decay`) and find beads near a point (`bp_near`). They are ordinary scripts; the app lists them, and a script can ignore them and build `sys` directly. The full list of `sys` fields is in the app under "sys fields".
34## Time integration
36The conservative part (springs and force pulses) is integrated with velocity Verlet, or with Yoshida's fourth-order composition of it (`sys.method = 'yoshida4'`). Damping is linear in the velocities, dv/dt = -A v, and is applied as a separate step v <- exp(-tau A) v, computed to round-off with a scaled Taylor series, once per audio sample (Strang splitting). The splitting error scales with the damping rate times the sample interval, which is small for the damping of interest here.
38The time step is chosen automatically. We estimate the highest frequency of the system linearized about rest by power iteration, take `sys.safety` (default 0.5) of the method's stability limit, and round so that a whole number of steps fits in each audio sample (48 kHz by default). The in-plane (stretching) modes of a net are much stiffer than its transverse modes, so they, not the audible frequencies, usually set the step; a smaller pre-stretch makes this worse.
40Each run reports:
42- the integrator's frequency error at 1 to 20 kHz, in cents (exact for a linear oscillator; with the default step for the drum example, Verlet is about 3.4 cents sharp at 5 kHz and 0.14 cents at 1 kHz);
43- the energy drift after the last excitation, with the energy removed by damping added back (it is counted exactly);
44- optionally, the relative difference of the listener signal from a run at dt/2 over [0, t] for several t. In a chaotic run this difference grows with t no matter how small dt is, so the short windows are the meaningful ones.
46## WebGPU
48The "WebGPU" toggle runs the time integration on the GPU (`src/sim/gpu.ts`, `gpu.wgsl`): the same scheme as on the CPU, in single precision, since WebGPU has no double precision. To keep single precision accurate, the GPU stores each bead's displacement from rest rather than its position, every element works from its rest geometry (computed in double precision on the CPU) plus displacement differences, and a spring's stretch is computed as (2 r0 . du + |du|^2)/(L + |r0|), without subtracting nearly equal lengths. Each phase of a step (element forces, bead update, damping terms, recording) is a dispatch over beads or elements; the phases for a chunk of audio samples are encoded into one command buffer, and the CPU computes energies and mode energies from the positions and velocities read back at each frame.
50Compared with the CPU (`src/devcompare.ts`, in Chrome on a laptop with an Intel Iris Xe GPU), the listener signals agree to single precision: relative differences of 2e-5 (drum) to 3e-3 (plate, dome) over the first 10 to 30 ms, with matching energy drift and dissipated energy; in a chaotic run the difference then grows like any other perturbation. The GPU has a floor of about 11 to 15x real time, the cost of about 15 dispatches per audio sample, so it is slower than the CPU for small nets (the 122-bead string, 14x against 3x) and faster for large ones: 367-bead dome 17x against 39x, 1015-bead plate 20x against 50x, a 6481-bead drum 97x against 304x. The perturbed twin needs double precision and is skipped when WebGPU is on.
52## Linear modes
54Linear modes are computed only when the "linear modes" toggle is on (off by default), or when the script sets `sys.modes = true`; the mode-energy and decay-time plots need them. The stiffness at rest is assembled exactly from the elements. For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, so its modes need only the z coordinates; the app computes them with numbl's `eig` (up to 1500 free beads; about 1 s for the 721 free beads of the drum example). For a curved net (a dome) the modes need all three coordinates of every free bead (up to 1200 coordinates, about 12 s). The app records the energy of every mode over time. The "energy by linear mode" plot shows the share of the quadratic energy in log-spaced frequency bands: in a linear system each band only decays, so any band that gains energy is fed by the nonlinearity.
56The "small-strike reference" option runs the same system with every hit scaled down by 1000 and the output scaled back up, which is what the system would do if it were linear. The app plays it for comparison and reports the difference from the actual run.
58## Chaos test
60The "perturbed twin" option runs the system again with every hit changed by one part in 10^9 and plots, over time, the largest difference in any bead coordinate between the two runs, relative to the largest displacement of the original run over the preceding 50 ms (so the measure stays meaningful as the motion decays). In chaotic motion this separation grows exponentially (a straight line on the log plot) until the runs are unrelated; in regular motion it grows slowly, roughly in proportion to time. The app fits the exponential rate over the stretch where the separation rises from 100 times the perturbation to 10^-2 and reports the e-folding time. Positions are compared at full precision, which needs the original run's positions at every frame in memory (limit 400 MB).
62## Examples
64| Example | What it shows |
65| --- | --- |
66| `drum` | A 15 cm drum head (823 beads, fundamental 237 Hz) struck by a soft mallet. Nearly linear at first (the listener's acceleration is 5% from the small-strike reference over 10 ms); the slight pitch change at large amplitude accumulates to 55% by 0.3 s. About 7x real time. |
67| `cascade` | The same net with a 0.5% pre-stretch, weaker dashpots and a hard strike. The share of energy in modes above 1 kHz grows from 4% (small strike) to about 25% within 30 ms. About 10x real time. |
68| `string` | 120 beads on a string, struck slightly off-axis; fast. |
69| `plate` | A flat plate of 1015 beads, 40 cm across, with a free rim and its center clamped (seven pinned beads), struck sharply near the rim: a flat cymbal or gong. Stiffness from bending triplets only; lowest modes 39 and 58 Hz, highest about 22 kHz. Bright (spectral centroid of the first 50 ms about 3.2 kHz). Mildly nonlinear at the given force, 1% from the small-strike reference over 10 ms and 14% by 0.3 s; at 4x the force, 14% and 93%. It is "thicker" than a real cymbal: nonlinear once the motion reaches about 4 mm. About 20x real time. |
70| `dome` | A domed cymbal: 367 beads, 40 cm across, with a 2 cm bow and a 2 cm bell, free at the rim and clamped at the top of the bell, its bending stiffness from hinges. Modes from 63 Hz to about 20 kHz (full 3D analysis, about 12 s). Curvature couples bending to stretching at small amplitude: at the given force the listener signal is 2.2% from the small-strike reference within 10 ms, against 0.3% for the same plate made flat; at 4x the force, the share of energy in modes above 2 kHz rises from 9% to 16% within 10 ms, while the flat plate's stays at 8%. About 20x real time. |
71| `bar` | A free bar of 41 beads on two soft supports at the nodes of its lowest mode, tuned to 880 Hz (it comes out at 838 Hz). Its overtones are at 2.75, 5.39 and 8.88 times the fundamental, against 2.76, 5.40 and 8.93 for an ideal free bar. |
72| `chaos` | The minimal chaotic example: two beads on a short string, free to move in 3D (6 degrees of freedom), struck hard and off-axis, with every mode decaying 60 dB in 2.5 s. With the perturbed twin, the separation e-folds every 2 ms and the runs are unrelated within 0.15 s; at 0.025 of the strike the motion is regular and the separation stays below 1e-8 for 3 s. With half the strike, the damping pulls the motion out of the chaotic range before the runs fully separate (transient chaos). Runs in a fraction of a second. |
74Timings are for a laptop (Chrome, single thread).
76## Limitations and next steps
78- Bending triplets need straight or flat rest shapes; hinges handle curved ones.
79- The solver is plain TypeScript in one Web Worker. Cost scales with the number of springs times the number of steps; nets of several thousand beads will need WASM or WebGPU.
80- A bead net has an upper frequency set by its spacing (about 4 kHz for the 1 cm drum example), so a dense high-frequency spectrum needs many beads.
81- Linear modes of curved nets need all three coordinates of each free bead and a dense eigensolver, so they are limited to about 400 free beads (1200 coordinates, about 12 s); flat nets go to 1500 free beads.