Woodblock physical model — specification#

Companion to woodblock.py. This describes the physics, the equations, the state vector, and the traps. The Python script is the reference implementation; where this document and the script disagree, the script is right.

The model produces free-field acoustic pressure in pascals at a point in space. There is no microphone, no directivity, no gain staging.

1. Signal chain#

mallet (Hunt–Crossley contact)
  → body bending modes (free–free bar, modal reduction)      ⇄ two-way coupled
      → air cavity (Helmholtz resonator)                     ⇄
          → distributed surface dipole sources + slit monopole
              → retarded free-field pressure at listener point
                  → anti-aliased decimation → PCM

Three physical sub-systems (body, cavity, mallet) integrate together as one coupled ODE. Radiation is a post-process: it is linear and non-reacting, so it runs once after the solve, not inside the loop.

2. Coordinate system and geometry#

Block-centred, metres. Needed for the 3D view and for radiation.

• Block occupies x ∈ [−L/2, +L/2], y ∈ [−W/2, +W/2], z ∈ [−H/2, +H/2] where L = length_m, W = width_m, H = block_height_m.
• Radiating faces sit at z = ±H/2.
• Strike point is at x_c = (strike_position − ½)·L on the top face.
• Slit opening is a point at (slit_x_m, slit_y_m, slit_z_m); default is on the side wall at (0, +W/2, 0).
• Listener is at (listener_x_m, listener_y_m, listener_z_m).
• wall_thickness_m (h) is the bending depth of the vibrating top wall. It is NOT the block height. It sets pitch; block_height_m sets the dipole arm.

For the 3D view: draw the box, the slot, a marker at the strike point, a marker at the slit, a marker at the listener, and optionally the sampled surface elements along the top face.

3. Body: free–free Euler–Bernoulli bar#

3.1 Modal frequencies#

Eigenvalues β_n L are the roots of cos(βL)·cosh(βL) = 1:

4.73004074, 7.85320462, 10.99560784, 14.13716549,  then (2n+1)π/2

The asymptotic form is exact to ~1e-9 from n = 5 up. With c_L = √(E/ρ), I = W h³/12, A = W h:

f_n = κ · (β_n L)² / (2π L²) · √(E I / (ρ A))
    = κ · (β_n L)² · c_L · h / (2π √12 · L²)

κ = hollow_softening accounts for stiffness removed by the slot.

Pitch scales linearly with h and as √(E/ρ). Verified: h = 6/10/16 mm gives 602/1003/1606 Hz; E = 6e9 vs 1.6e10 gives 614 vs 1003 Hz.

Optional stretch (γ = inharmonicity), applied to partials above the first:

f_n ← f_1 · (f_n / f_1)^γ

Discard modes above max_mode_hz (20 kHz). Defaults keep 6.

3.2 Mode shapes#

On s = x/L ∈ [0,1], with σ_n = (cosh βL − cos βL)/(sinh βL − sin βL):

φ̃_n(s) = cosh(β_n L s) + cos(β_n L s) − σ_n [ sinh(β_n L s) + sin(β_n L s) ]

Mass-normalise so the modal mass is 1:

φ_n(s) = φ̃_n(s) / ( √(∫₀¹ φ̃_n² ds) · √M_p )

M_p = participating_mass_kg. Units: φ_n is kg^−1/2, modal coordinate q_n is m·√kg, physical displacement w = Σ q_n φ_n is metres.

M_p is a lumped correction, the least principled parameter in the model. The bar formula otherwise locks frequency and modal mass together, and the thin-wall mass (ρWhL ≈ 93 g) makes the block far too compliant — 395 µm of deflection and 119 dB SPL at 42 cm. A real woodblock is 300–500 g with the whole box participating. M_p = 0.35 gives 124 µm and 109.7 dB. It sets absolute level, so treat SPL as plausible, not calibrated.

3.3 Damping#

ζ_n = ½ ( η₀ + η₁ ω_n )

η₀ = loss_factor (frequency-independent wood loss), η₁ = loss_slope (extra HF damping). Bigger = duller and shorter. η₀: 0.010 → T60 185 ms; 0.030 → 66 ms.

4. Mallet contact (the only nonlinearity)#

Surface displacement and velocity at the strike point:

w  = Σ_n q_n  φ_n(s_c)
ẇ  = Σ_n q̇_n φ_n(s_c)
δ  = y − w                     (compression; y = mallet position)

Hunt–Crossley, with λ = contact_damping in s/m:

f_c = K δ^α (1 + λ (ẏ − ẇ))     if δ > 0, clamped to ≥ 0
f_c = 0                          otherwise

m_h ÿ = −f_c

Initial conditions: y = 0, ẏ = v₀ = mallet_velocity, all else zero.

Trap. Do not write the loss term as μ · d(δ^α)/dt. That makes μ a relaxation time; at μ = 0.35 s against a 0.1 ms collision, a 3 µm compression produced 1.9 × 10⁷ N. Every strike became a delta function and mallet stiffness had literally zero effect on the spectrum. The form above keeps the correction dimensionless and O(1).

Expected behaviour, worth asserting in tests: contact duration scales roughly as K^−2/5. K = 3e7 / 3e8 / 3e9 / 3e10 → 0.92 / 0.59 / 0.28 / 0.16 ms and 287 / 615 / 1289 / 3012 N. Because the contact is nonlinear, timbre varies with strike velocity on its own — that is physics, not a scripted velocity layer.

5. Cavity: Helmholtz resonator#

Effective plug length with end corrections on both ends:

L_eff = slit_length_m + 2 · 0.85 · √(S/π)          S = slit_area_m2
ω_H   = c √( S / (V L_eff) )                       V = cavity_volume_m3

Coupling coefficients — net volume each mode injects into the cavity, over the wall span [lo, hi] = cavity_wall_span:

A_n = χ · W L ∫_lo^hi φ_n(s) ds            χ = cavity_coupling

Cavity pressure, from wall motion plus slit-plug displacement ξ:

p_c = (ρ₀ c² / V) ( S ξ + Σ_n A_n q_n )

Plug equation (mass per unit area ρ₀ L_eff, driven by −p_c):

ξ̈ = −2 ζ_H ω_H ξ̇ − ω_H² ξ − (c² / (V L_eff)) Σ_n A_n q_n

A_n appears both in the cavity drive and in the modal load below. That symmetry is the correctness check — if the coupling only goes one way, the energy will not be conserved.

6. Coupled system#

State vector, length 2N + 4:

[ q₁..q_N , q̇₁..q̇_N , ξ , ξ̇ , y , ẏ ]

Derivatives:

q̈_n = −2 ζ_n ω_n q̇_n − ω_n² q_n + f_c φ_n(s_c) − A_n p_c
ξ̈   = −2 ζ_H ω_H ξ̇ − ω_H² ξ − (c²/(V L_eff)) Σ_m A_m q_m
ÿ   = −f_c / m_h

Integrate with RK4 at f_hi = sample_rate × oversample. Record q̈_n(t) for all modes and ξ̈(t) — those are the radiation inputs.

The system is linear apart from f_c. Eigenvalues of the coupled linear part are all strictly decaying (slowest ≈ −34 s⁻¹ at defaults); a zero or positive real part means the coupling was assembled wrongly.

Convergence. Relative to oversample = 32: oversample = 8 differs by 5.2%, 16 by 1.8%. Use 16 for render, 8 is acceptable for interactive preview. The residual sits in the top mode near 18.7 kHz.

7. Radiation: pressure at a point#

Each cross-section of a bending bar translates, so the top and bottom faces inject equal and opposite volume. Every slice is a local dipole, and the net volume displacement of the bar is exactly zero.

Trap. This is why you cannot use a single lumped monopole term with B_n = ∫φ_n ds for the body: elastic modes are orthogonal to rigid translation, so that integral is identically zero and the body contributes nothing. Near-field emphasis and HF directivity must come from path lengths and delays, not from a gain constant.

Discretise the top face into M = radiating_elements strips at s_j = (j − ½)/M, x_j = (s_j − ½)L, each of area patch = W L / M. Every strip contributes a pair of point sources:

r_j⁺ = (x_j, 0, +H/2)     sign +1
r_j⁻ = (x_j, 0, −H/2)     sign −1

with volume acceleration

V̈_j(t) = ± patch · Σ_n φ_n(s_j) q̈_n(t)

Plus one monopole at the slit: V̈_slit(t) = S ξ̈(t).

Sum retarded monopoles over all 2M + 1 sources:

p(r, t) = Σ_j  ρ₀ / (4π d_j) · V̈_j( t − d_j / c )        d_j = |r − r_j|

Integer-sample delays are fine: at f_hi = 768 kHz one sample is 0.45 mm.

Element count. Spacing must satisfy L/M < ¼ · c/f_max; at defaults that needs M ≥ 44, so M = 48. Warn the user if M is too coarse, and if the listener is inside the source region (|r| < L/2) where point-source strips are crude.

Validation — the app should reproduce these#

Falloff is −7.8 dB per doubling close in, relaxing to −6.5 dB by 1.6 m: the dipole-to-inverse-square transition. If your port shows a flat −6 dB everywhere, the two faces are not in antiphase.

8. Sampling and output#

1. Solve at f_hi = 48000 × oversample.
2. Radiate → p(t) in pascals at f_hi.
3. Anti-alias and decimate by oversample. Python uses scipy.resample_poly; in the browser implement a windowed-sinc (Kaiser, β ≈ 8) FIR lowpass at ~0.45·48 kHz with roughly 12 × oversample taps, then take every oversample-th sample. Do not decimate without the filter — the double differentiation puts real energy near f_hi/2.
4. Prepend pre_delay_ms of silence, apply a short squared fade at the tail.
5. Either normalise to peak_dbfs, or divide by full_scale_pa for absolute pressure. Report peak Pa and dB SPL (20 log₁₀(p/20 µPa)), and warn on clip.

Trap. Do not "remove DC" by subtracting the mean of the first N samples. The strike is at t = 0, so the head of the buffer is the loudest part of the file; subtracting it injects a large offset across the whole signal. The model has no DC — p.mean() is ~1e-11 — so no DC removal is needed at all.

9. Parameters#

Defaults yield: partials 1003 / 2766 / 5423 / 8964 / 13390 / 18702 Hz (ratios 1.00, 2.76, 5.40, 8.93 — bar inharmonicity, not octaves), Helmholtz 484 Hz, contact 0.30 ms at 794 N, T60 ≈ 149 ms, 109.7 dB SPL.

10. Browser implementation notes#

• Cost. Solve is 4 × T derivative evaluations with N ≈ 6 (T = 345,600 at defaults) — cheap. Radiation is the bottleneck: 2M+1 sources × N × T ≈ 200 M multiply-adds. Compute Σ_n φ_n(s_j) q̈_n(t) once per strip and reuse it for the +z and −z sources with opposite sign — that halves the work.
• Run in a Web Worker, transfer the result as a Float32Array ArrayBuffer, and copy into an AudioBuffer at 48 kHz. Keep the solver in Float64Array.
• For live parameter dragging, preview at oversample = 8 and radiating_elements = 24, then re-render at full quality on release.
• Nothing needs real-time streaming — one strike is one offline render of well under a second.
• Cache the modal basis; it only depends on geometry, material, strike_position, and radiating_elements, not on mallet or listener parameters. Moving the listener requires only re-running §7.
• Good things to visualise: mode shapes animated as φ̃_n(s) along the bar; partial frequencies as a spectrum overlay; the contact-force pulse f_c(t) (its width sets brightness); source strips and listener in the 3D view.