Listener quantity: acceleration (default), velocity or displacement; acceleration from forces with boundary-corrected damping; zero-mode threshold 1e-6 of top frequency
10 changed files+118−21
README.mdmodified+4−4View file
@@ -1,6 +1,6 @@
11 # bead-percussion
22
3-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](https://numbl.org), integrated accurately in time, and heard through the displacement of a listener bead.
3+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](https://numbl.org), integrated accurately in time, and heard through the motion (by default the acceleration) of a listener bead.
44
55 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.
66
@@ -22,7 +22,7 @@ Each bead has a mass and a position in 3D. The script defines:
2222 - *Friction* on each bead, force -c0 v, which damps all frequencies alike.
2323 - *Pinned* beads, held fixed.
2424 - *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.
25-- A *listener*: a bead (or the average of several) and a direction; the displacement along that direction is the audio signal.
25+- 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.
2626
2727 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. 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.
2828
@@ -56,10 +56,10 @@ The "perturbed twin" option runs the system again with every hit changed by one
5656
5757 | Example | What it shows |
5858 | --- | --- |
59-| `drum` | A 15 cm drum head (823 beads, fundamental 237 Hz) struck by a soft mallet. Nearly linear at first (2% from the small-strike reference over 10 ms); the slight pitch change at large amplitude accumulates to 25% by 0.3 s. About 7x real time. |
59+| `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. |
6060 | `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. |
6161 | `string` | 120 beads on a string, struck slightly off-axis; fast. |
62-| `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. Mildly nonlinear at the given force (it is "thicker" than a real cymbal: nonlinear once the motion reaches about 4 mm). About 20x real time. |
62+| `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. |
6363 | `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. |
6464 | `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. |
6565
matlab/examples/plate.mmodified+4−3View file
@@ -27,9 +27,10 @@ sys.bends = [trip, kb * ones(size(trip, 1), 1)];
2727 bp_decay(pos, sys.mass, sys.springs, [100 2.5; 5000 1.2], sys.bends);
2828
2929 % A stick on the edge: a sharp 0.3 ms push on the beads within 2 cm of the
30-% rim point. A flat plate is only mildly nonlinear at this force (a few
31-% percent from the small-strike reference); at 4x it differs by about 10%
32-% after 0.3 s. A longer, softer push mostly rocks the plate on its mount.
30+% rim point. The plate is mildly nonlinear at this force (14% from the
31+% small-strike reference by 0.3 s); at 4x it differs by 14% within 10 ms
32+% and 93% by 0.3 s. A longer, softer push mostly rocks the plate on its
33+% mount.
3334 hit.t = 0;
3435 hit.beads = bp_near(pos, [0.18 0 0], 0.02);
3536 hit.force = [0 0 -150]; % peak force per bead (N)
src/App.tsxmodified+4−1View file
@@ -379,7 +379,10 @@ function Report({ p }: { p: RunPayload }) {
379379 const SYS_DOC = `Required
380380 sys.pos n-by-3 bead positions (m)
381381 sys.mass scalar or n-by-1 (kg)
382- sys.listen struct: beads (index or indices, averaged), dir (default [0 0 1])
382+ sys.listen struct: beads (index or indices, averaged), dir (default [0 0 1]),
383+ quantity: 'acceleration' (default), 'velocity' or 'displacement'
384+ (sound in air follows acceleration more closely; displacement
385+ favors low frequencies by 1/f^2)
383386 Elements
384387 sys.springs rows [i j k L0] or [i j k L0 k3]; force k*e + k3*e^3, e = length - L0
385388 sys.bends rows [i j k kb] or [i j k kb cb], j the middle bead: bending,
src/sim/analysis.tsmodified+1−1View file
@@ -91,7 +91,7 @@ export function modeBands(
9191 ): { edges: Float64Array; e: Float32Array; nBands: number } {
9292 const nModes = freqs.length;
9393 // modes at (nearly) 0 Hz, from free rigid motion, go in the lowest band
94- const fLow = freqs.find(f => f > 1e-4 * freqs[nModes - 1]) ?? freqs[nModes - 1];
94+ const fLow = freqs.find(f => f > 1e-6 * freqs[nModes - 1]) ?? freqs[nModes - 1];
9595 const fmin = fLow * 0.97;
9696 const fmax = freqs[nModes - 1] * 1.03;
9797 const edges = new Float64Array(nBands + 1);
src/sim/parse.tsmodified+6−1View file
@@ -2,7 +2,7 @@
22 // The script uses 1-based bead indices; they become 0-based here.
33
44 import type { RuntimeValue } from "numbl";
5-import type { BeadSystem, Hit, Method } from "./types";
5+import type { BeadSystem, Hit, ListenQuantity, Method } from "./types";
66
77 type Mat = { rows: number; cols: number; data: Float64Array };
88
@@ -220,6 +220,10 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
220220 const listenBeads = indices(L.get("beads") ?? L.get("bead"), "sys.listen.beads", n);
221221 if (listenBeads.length === 0) throw new SysError("sys.listen.beads is required");
222222 const dir = vec3(L.get("dir"), "sys.listen.dir", [0, 0, 1]);
223+ const listenQuantity = str(L.get("quantity"), "sys.listen.quantity", "acceleration") as ListenQuantity;
224+ if (!["acceleration", "velocity", "displacement"].includes(listenQuantity)) {
225+ throw new SysError("sys.listen.quantity must be 'acceleration' (default), 'velocity' or 'displacement'");
226+ }
223227 const dn = Math.hypot(...dir);
224228 if (dn === 0) throw new SysError("sys.listen.dir must be nonzero");
225229
@@ -229,6 +233,7 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
229233 const sys: BeadSystem = {
230234 n, pos0, mass, pinned, si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb, di, dj, dc, friction, hits, listenBeads,
231235 listenDir: [dir[0] / dn, dir[1] / dn, dir[2] / dn],
236+ listenQuantity,
232237 duration: scalar(f.get("duration"), "sys.duration", 2),
233238 sampleRate: scalar(f.get("sampleRate"), "sys.sampleRate", 48000),
234239 dt: scalar(f.get("dt"), "sys.dt", 0),
src/sim/solver.tsmodified+44−7View file
@@ -331,14 +331,51 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
331331 const lb = sys.listenBeads;
332332 const [ldx, ldy, ldz] = sys.listenDir;
333333 const pos0 = sys.pos0;
334- const listen = () => {
334+ // Acceleration is (F - damping force) / m on each listener bead, with F the
335+ // spring, bending and pulse forces at the current positions and time (the
336+ // force array is current at every sample boundary). The damping force uses
337+ // the dashpots and bending triplets that touch each listener bead.
338+ const touchD = Array.from(lb, b => Array.from(di.keys()).filter(d => di[d] === b || dj[d] === b));
339+ const touchB = Array.from(lb, b => Array.from(bi.keys()).filter(t => (bcb[t] !== 0 && (bi[t] === b || bj[t] === b || bk[t] === b))));
340+ const dampForce = (v: Float64Array, q: number, b: number, c: number) => {
341+ let f = -friction[b] * v[3 * b + c];
342+ for (const d of touchD[q]) {
343+ const other = di[d] === b ? dj[d] : di[d];
344+ f -= dc[d] * (v[3 * b + c] - v[3 * other + c]);
345+ }
346+ for (const t of touchB[q]) {
347+ const r = v[3 * bi[t] + c] - 2 * v[3 * bj[t] + c] + v[3 * bk[t] + c];
348+ const w = bj[t] === b ? -2 : 1;
349+ f -= bcb[t] * w * r;
350+ }
351+ return f;
352+ };
353+ // At a sample boundary the damping for the time pendingD has not been
354+ // applied yet; vb = v - pendingD A v brings the velocities to the boundary
355+ // (to second order) for the listener's velocity and damping force.
356+ const vb = new Float64Array(N3);
357+ const Av = new Float64Array(N3);
358+ const listen = (pendingD: number) => {
335359 let s = 0;
336- for (const b of lb) {
337- s +=
338- (x[3 * b] - pos0[3 * b]) * ldx +
339- (x[3 * b + 1] - pos0[3 * b + 1]) * ldy +
340- (x[3 * b + 2] - pos0[3 * b + 2]) * ldz;
360+ if (sys.listenQuantity !== "displacement") {
361+ if (pendingD !== 0 && hasDamping) {
362+ applyA(v, Av);
363+ for (let q = 0; q < N3; q++) vb[q] = v[q] - pendingD * Av[q];
364+ } else vb.set(v);
341365 }
366+ lb.forEach((b, q) => {
367+ const dir = [ldx, ldy, ldz];
368+ for (let c = 0; c < 3; c++) {
369+ if (dir[c] === 0) continue;
370+ let val: number;
371+ if (sys.listenQuantity === "displacement") val = x[3 * b + c] - pos0[3 * b + c];
372+ else {
373+ if (sys.listenQuantity === "velocity") val = vb[3 * b + c];
374+ else val = (F[3 * b + c] + dampForce(vb, q, b, c)) * invm[b];
375+ }
376+ s += val * dir[c];
377+ }
378+ });
342379 return s / lb.length;
343380 };
344381 const frameStride = Math.max(1, Math.round(fs / sys.frameRate));
@@ -428,7 +465,7 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
428465 }
429466 }
430467 }
431- if (sample >= 0 && sample < nSamples) audio[sample] = listen();
468+ if (sample >= 0 && sample < nSamples) audio[sample] = listen(pendingD);
432469 if (isFrame) {
433470 frameTimes[frame] = t;
434471 energy[frame] = kinetic() + pe - peRest;
src/sim/types.tsmodified+2−0View file
@@ -2,6 +2,7 @@
22 // script uses 1-based indices and parse.ts converts them.
33
44 export type Method = "verlet" | "yoshida4";
5+export type ListenQuantity = "acceleration" | "velocity" | "displacement";
56
67 export interface Hit {
78 t: number; // start time (s)
@@ -39,6 +40,7 @@ export interface BeadSystem {
3940 hits: Hit[];
4041 listenBeads: Int32Array;
4142 listenDir: [number, number, number]; // unit vector
43+ listenQuantity: ListenQuantity;
4244 duration: number; // s
4345 sampleRate: number; // Hz, audio output
4446 dt: number; // s; 0 = automatic
src/ui/Plots.tsxmodified+5−3View file
@@ -32,7 +32,9 @@ export function Waveform({ p, cursor }: { p: RunPayload; cursor: number | null }
3232 for (const v of audio) peak = Math.max(peak, Math.abs(v));
3333 peak = peak || 1;
3434 const T = audio.length / sampleRate;
35- drawAxes(f, [0, T], [-peak * 1e3, peak * 1e3], { xlabel: "time (s)", ylabel: "disp (mm)" });
35+ const q = p.sys.listenQuantity;
36+ const [unit, label] = q === "displacement" ? [1e3, "disp (mm)"] : q === "velocity" ? [1, "vel (m/s)"] : [1, "accel (m/s²)"];
37+ drawAxes(f, [0, T], [-peak * unit, peak * unit], { xlabel: "time (s)", ylabel: label });
3638 const { ctx, x0, y0, pw, ph } = f;
3739 ctx.fillStyle = f.css("--series");
3840 const per = audio.length / pw;
@@ -311,10 +313,10 @@ export function DecayTimes({ p }: { p: RunPayload }) {
311313 c => {
312314 const f = setupCanvas(c);
313315 const fTop = freqs[freqs.length - 1];
314- const shown = (k: number) => freqs[k] > 1e-4 * fTop; // skip free rigid motion at 0 Hz
316+ const shown = (k: number) => freqs[k] > 1e-6 * fTop; // skip free rigid motion at 0 Hz
315317 const T60 = Array.from(sig, (s, k) => (s > 0 && shown(k) ? Math.log(1000) / s : Infinity));
316318 const finite = T60.filter(Number.isFinite);
317- const fx0 = Math.log10((freqs.find(f => f > 1e-4 * fTop) ?? fTop) * 0.9);
319+ const fx0 = Math.log10((freqs.find(f => f > 1e-6 * fTop) ?? fTop) * 0.9);
318320 const fx1 = Math.log10(freqs[freqs.length - 1] * 1.1);
319321 let ty0 = Math.floor(Math.log10(Math.min(...finite, 1)) * 2) / 2;
320322 let ty1 = Math.ceil(Math.log10(Math.max(...finite, 1)) * 2) / 2;
src/worker.tsmodified+3−1View file
@@ -22,6 +22,7 @@ export interface RunPayload {
2222 sj: Int32Array;
2323 nDashpots: number;
2424 nBends: number;
25+ listenQuantity: string;
2526 listenBeads: Int32Array;
2627 hitBeads: Int32Array;
2728 method: string;
@@ -75,7 +76,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
7576 post({ type: "stage", text: `Computing linear modes (${nFree} free beads)` });
7677 modes = computeModes(sys, numblEig);
7778 const fTop = modes.freqs[modes.freqs.length - 1];
78- const nZero = modes.freqs.filter(f => f < 1e-4 * fTop).length;
79+ const nZero = modes.freqs.filter(f => f < 1e-6 * fTop).length;
7980 if (nZero > 0) {
8081 notes.push(
8182 `${nZero} transverse mode${nZero > 1 ? "s are" : " is"} at (nearly) 0 Hz: the beads can move or tilt freely as a whole. ` +
@@ -155,6 +156,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
155156 sj: sys.sj,
156157 nDashpots: sys.dc.filter(c => c !== 0).length,
157158 nBends: sys.bi.length,
159+ listenQuantity: sys.listenQuantity,
158160 listenBeads: sys.listenBeads,
159161 hitBeads,
160162 method: sys.method,
tests/solver.test.tsmodified+45−0View file
@@ -335,3 +335,48 @@ sys.friction = 0; sys.duration = 0.01;`);
335335 [2.756, 5.404, 8.933].forEach((ratio, q) => expect(Math.abs(f[3 + q] / f1 / ratio - 1)).toBeLessThan(0.01));
336336 });
337337 });
338+
339+describe("listener quantity", () => {
340+ // Stiff, heavily damped string (dashpots and bending damping), listened to
341+ // as displacement, velocity and acceleration at the same bead.
342+ const run = (q: string) => {
343+ const { sys } = build(
344+ STRING(`
345+[pos, E, trip] = bp_chain(N + 2, (N + 1) * a);
346+sys.bends = [trip, 50 * ones(size(trip, 1), 1)];
347+[sys.friction, sys.dashpots, sys.bends] = bp_decay(pos, m, sys.springs, [100 0.05; 1000 0.02], sys.bends);
348+sys.hits.velocity = [0 0 1e-4];
349+sys.listen.quantity = '${q}';`)
350+ );
351+ return simulate({ ...sys, duration: 0.05 });
352+ };
353+
354+ it("defaults to acceleration", () => {
355+ expect(build(STRING("")).sys.listenQuantity).toBe("acceleration");
356+ });
357+
358+ it("gives acceleration and velocity consistent with each other and the displacement", () => {
359+ const x = run("displacement").audio;
360+ const v = run("velocity").audio;
361+ const acc = run("acceleration").audio;
362+ const fs = 48000;
363+ const rel = (p: (i: number) => number, ref: Float32Array) => {
364+ let num = 0;
365+ let den = 0;
366+ for (let i = 1; i < ref.length - 1; i++) {
367+ num += (p(i) - ref[i]) ** 2;
368+ den += ref[i] ** 2;
369+ }
370+ return Math.sqrt(num / den);
371+ };
372+ // central differences of the recorded signals; their own error is second
373+ // order in 1/fs (about 5e-4 for v -> a here, and 1e-5 for x -> v, which is
374+ // limited by the single-precision audio)
375+ expect(rel(i => ((x[i + 1] - x[i - 1]) * fs) / 2, v)).toBeLessThan(1e-4);
376+ expect(rel(i => ((v[i + 1] - v[i - 1]) * fs) / 2, acc)).toBeLessThan(1e-3);
377+ });
378+
379+ it("rejects an unknown quantity", () => {
380+ expect(() => build(STRING("sys.listen.quantity = 'jerk';"))).toThrow(/quantity/);
381+ });
382+});