bead-percussion#

A 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, integrated accurately in time, and heard through the motion (by default the acceleration) of a listener bead.

The 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.

Running#

npm install
npm run dev      # the app, at the printed localhost URL
npm test         # solver tests (Node; the WebGPU path is checked in the browser with src/devcompare.ts)

The model#

Each bead has a mass and a position in 3D. The script defines:

• 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.
• 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).
• 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.
• 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.
• Friction on each bead, force -c0 v, which damps all frequencies alike.
• Pinned beads, held fixed.
• 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.
• 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.

Damping 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.

A 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.

The 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".

Time integration#

The 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.

The 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.

Each run reports:

• 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);
• the energy drift after the last excitation, with the energy removed by damping added back (it is counted exactly);
• 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.

WebGPU#

The "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.

Compared 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.

Linear modes#

Linear 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.

The "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.

Chaos test#

The "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).

Examples#

Timings are for a laptop (Chrome, single thread).

Limitations and next steps#

• Bending triplets need straight or flat rest shapes; hinges handle curved ones.
• 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.
• 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.
• 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.