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 displacement 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 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.
• Dashpots between pairs of beads, with force c (v_j - v_i) on bead i. 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) and a direction; the displacement along that direction is the audio signal.

Damping is easiest to set by its effect. [sys.friction, sys.dashpots] = bp_decay(pos, mass, springs, [f1 T1; f2 T2]) chooses them 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 and each spring a parallel dashpot 2 sigma1 times its tension per length, which makes the damping matrix of the transverse motion 2 sigma0 M + 2 sigma1 K; the modes then stay uncoupled and the law holds exactly at small amplitude. Large motion and the in-plane modes are damped by the same elements, but not by this law. The page plots the small-amplitude decay time of every mode for whatever damping the script sets.

Bending stiffness (for plates and cymbals) is not in yet.

The matlab/helpers/bp_*.m functions build common geometries (bp_hexdisk: a triangular lattice clipped to a disk, with its rim beads; bp_chain: a line of beads), list edges (bp_edges), make springs with a given tension at rest (bp_springs), 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.

Linear modes#

For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, and its stiffness matrix is a graph Laplacian weighted by spring tension over length. The app computes its eigenpairs with numbl's eig (up to 1500 free beads; about 1 s for the 721 free beads of the drum example) and 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. 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#

• No bending stiffness yet, so no plates or cymbals; the plan is hinge (dihedral angle) elements on pairs of triangles, with a domed rest shape.
• 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 are computed only for flat nets.