# 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. | Axis | Direction | |---|---| | +x | along the bar's long axis | | +y | across the width | | +z | up, out of the struck face | - 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 | Listener (x, y, z) m | Peak | Note | |---|---|---| | 0.15, 0.25, 0.30 | 109.7 dB | default, oblique | | 0.00, 0.00, 0.40 | 100.9 dB | above centre; 2766 Hz partial −27 dB (symmetry null) | | 0.00, 0.40, 0.00 | 90.9 dB | **edge-on**: faces equidistant, body radiation cancels — 2766 Hz −61 dB, 8964 Hz −90 dB, leaving near-pure cavity monopole | | 0.00, 0.05, 0.02 | 125.4 dB | close to slit, cavity-dominated | 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 | Parameter | Default | Effect | |---|---|---| | `length_m` | 0.200 | pitch ∝ 1/L² | | `width_m` | 0.055 | source area / level | | `wall_thickness_m` | 0.010 | **pitch ∝ h** | | `block_height_m` | 0.050 | dipole arm; bigger = more low end | | `density` | 850 | pitch ∝ 1/√ρ | | `youngs_modulus` | 1.6e10 | pitch ∝ √E | | `hollow_softening` | 0.90 | overall pitch trim | | `participating_mass_kg` | 0.35 | level ∝ 1/M (lumped; see §3.2) | | `inharmonicity` | 1.00 | >1 stretches partials (2766 → 3061 Hz at 1.10) | | `loss_factor` | 0.010 | T60 (185 ms → 66 ms at 0.030) | | `loss_slope` | 1.2e-7 | HF damping; main brightness control | | `strike_position` | 0.42 | partial balance. **Avoid 0.5** — a node of every antisymmetric mode | | `mallet_mass` | 0.025 | contact time, level | | `mallet_velocity` | 2.5 | dynamics *and* timbre (nonlinear contact) | | `mallet_stiffness` | 3.0e8 | hardness; centroid 1953 → 5550 Hz at 3e10 | | `contact_exponent` | 1.5 | 1.5 Hertzian, 2–3 felt | | `contact_damping` | 0.8 | Hunt–Crossley λ, **s/m not seconds** | | `cavity_volume_m3` | 2.6e-4 | Helmholtz pitch | | `slit_area_m2` | 8.0e-4 | Helmholtz pitch, cavity output | | `slit_length_m` | 0.012 | Helmholtz pitch | | `cavity_damping` | 0.06 | hollow-resonance Q | | `cavity_coupling` | 1.0 | 0 = dry bar, >1 = hollower | | `listener_*_m` | .15/.25/.30 | see §7 | | `radiating_elements` | 48 | accuracy, not timbre (see §7) | | `oversample` | 16 | accuracy vs speed | 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.