physmod

A pseudo-physical model of resonance in an n-dimensional space, built as a finite-difference scheme for sound synthesis.

A grid of points, each displaced into one extra dimension and pulled back by its neighbours. Sample the displacement at a tap once per step and you have an audio signal.

use physmod::{Axis, Medium};

let sr = 44_100.0;
let mut m = Medium::tuned(vec![Axis::fixed(128)], &[1], 110.0, sr)?;

for k in 1..=24 {
    m.excite_mode(&[k], 1.0 / k as f64)?;   // roughly a plucked string
}

let tap = m.tap_at(&[19]).unwrap();
let audio: Vec<f64> = (0..88_200).map(|_| { m.step(); m.read(tap) }).collect();

What governs the thing

The time stepping is Stormer-Verlet / leapfrog, written in its two-variable form. Eliminating velocity recovers the classical three-level scheme x[n+1] - 2x[n] + x[n-1] = c^2 L x[n], so displacements are second-order accurate and the whole classical wave-equation literature applies directly.

One coefficient, not two. The wave speed, in grid cells per step. Older formulations of this model carried separate momentum and pull coefficients, but only their product is physical -- rescaling one and inversely rescaling the other leaves the trajectory bit-identical. What survives is c = sqrt(momentum * pull), and low modes obey omega = c |k|.

The wave speed is capped. Stability requires c <= 1/sqrt(rank), the CFL condition. Medium::new enforces it, so violations surface as an Error at construction rather than as inf in your output.

Pitch has a ceiling, and long grids play low. One period is a round trip across the medium, so pitch rises with wave speed -- which is capped. For a fixed axis of N cells the ceiling is about sample_rate * c_max / 2N. Ask a long grid for a high pitch and tuning fails; shorten the axes. pitch::highest_hz reports the limit.

Nothing damps. Below the CFL limit the scheme is symplectic and conserves Medium::energy exactly -- relative drift of order 1e-15 over 88,200 steps. A plucked medium rings forever. A simulation that seems to die out is almost certainly tuned so low that its fundamental is a fraction of a hertz.

The conserved quantity is

E = c^2 <x, -Lx> + c^2 <x, -Lv> + <v, v>

The cross term is what makes it exact; kinetic-plus-potential merely oscillates. It is positive definite precisely when the CFL condition holds, so it doubles as the stability proof and as a runtime monitor with no false positives.

Bending pitch reintroduces the danger. The CFL bound assumes a constant wave speed. Varying it can pump a mode through parametric resonance, diverging while every instantaneous value stays legal. Modulating at twice a mode's frequency by 1% grows it by 4e5 over 120k steps. Glide::ConstantEnergy and Glide::ConstantAction rescale the state across the change and hold growth to 1.00x; Glide::Free is fine only for bends slow relative to the modes carrying energy.

Boundaries

Boundary::Fixed, Periodic, and Free are Dirichlet, periodic and Neumann. Set per axis end, so Axis::fixed gives the full harmonic series while Axis::stopped (fixed one end, free the other) gives odd harmonics only, an octave below -- the stopped-pipe series.

Instrument files

Geometry, tuning and excitation are what you change while hunting for a sound, so they live in TOML rather than in code. The config feature adds a render binary that turns a definition into a WAV:

cargo run --release --features config --bin render -- instruments/string.toml

sample_rate = 44100
duration    = 2.0
output      = "string.wav"

[[axis]]
len = 128
lower = "fixed"
upper = "fixed"

[tune]
mode = [1]
hz   = 110.0

[[excite]]
mode      = [1]     # or: point = [21], to strike instead of exciting a mode
amplitude = 1.0

[[tap]]
at = [19]           # one tap renders mono, two render stereo

Add hz_end and glide to [tune] for a pitch bend; the sweep is exponential in pitch, so it is linear in cents. Unreachable pitches are rejected before the render starts rather than after it.

instruments/ has five to start from: string (full harmonic series), clarinet (stopped, odd harmonics only), bend, drum (2-D, stereo) and four_d.

Note one TOML rule: a bare key belongs to the table header above it, so top-level keys must come before the first [[axis]].

Examples

The examples demonstrate the API rather than the workflow; for making sounds, prefer the instrument files above.

cargo run --release --example pluck     # a string, writes pluck.wav
cargo run --release --example bend      # glide modes and parametric resonance
cargo run --release --example four_d    # a 4-D medium, writes four_d.wav
cargo run --release --example bench     # throughput at several grid sizes

Performance

Roughly 2e8 point-updates per second single-threaded. The optional parallel feature spreads each step across threads with rayon:

cargo run --release --features parallel --example bench

On 16 cores that is about 4x on large grids -- the scheme is memory-bandwidth bound, not compute bound -- and a large loss below roughly 10k points, where per-step scheduling overhead dominates. Enable it only for big grids.

Memory is 16 bytes per point, two f64 fields and no neighbour cache; neighbour offsets are computed arithmetically from strides.

History

This started as a C program, kept on the old branch. The rewrite was prompted by a coefficient that seemed impossible to tune -- vibrations either died out at once or blew up -- which turned out to be two separate problems. The 4-D configuration sat about 4x above the CFL limit, so it could not have been stable at any setting; and the alternative it was compared against was tuned so low its fundamental was 0.22 Hz, which over a 64k-sample run is a third of one cycle and reads as silence. The usable range between those was wide and never sampled.

git show old:physmod.c for the original.