1# Woodblock physical model — specification
3Companion to `woodblock.py`. This describes the physics, the equations, the state
4vector, and the traps. The Python script is the reference implementation; where
5this document and the script disagree, the script is right.
7The model produces **free-field acoustic pressure in pascals at a point in
8space**. There is no microphone, no directivity, no gain staging.
10---
12## 1. Signal chain
14```
15mallet (Hunt–Crossley contact)
16 → body bending modes (free–free bar, modal reduction) ⇄ two-way coupled
17 → air cavity (Helmholtz resonator) ⇄
18 → distributed surface dipole sources + slit monopole
19 → retarded free-field pressure at listener point
20 → anti-aliased decimation → PCM
21```
23Three physical sub-systems (body, cavity, mallet) integrate together as one
24coupled ODE. Radiation is a post-process: it is linear and non-reacting, so it
25runs once after the solve, not inside the loop.
27---
29## 2. Coordinate system and geometry
31Block-centred, metres. Needed for the 3D view and for radiation.
33| Axis | Direction |
34|---|---|
35| +x | along the bar's long axis |
36| +y | across the width |
37| +z | up, out of the struck face |
39- Block occupies `x ∈ [−L/2, +L/2]`, `y ∈ [−W/2, +W/2]`, `z ∈ [−H/2, +H/2]`
40 where `L = length_m`, `W = width_m`, `H = block_height_m`.
41- **Radiating faces** sit at `z = ±H/2`.
42- **Strike point** is at `x_c = (strike_position − ½)·L` on the top face.
43- **Slit opening** is a point at `(slit_x_m, slit_y_m, slit_z_m)`; default is on
44 the side wall at `(0, +W/2, 0)`.
45- **Listener** is at `(listener_x_m, listener_y_m, listener_z_m)`.
46- `wall_thickness_m` (`h`) is the *bending depth* of the vibrating top wall. It
47 is NOT the block height. It sets pitch; `block_height_m` sets the dipole arm.
49For the 3D view: draw the box, the slot, a marker at the strike point, a marker
50at the slit, a marker at the listener, and optionally the sampled surface
51elements along the top face.
53---
55## 3. Body: free–free Euler–Bernoulli bar
57### 3.1 Modal frequencies
59Eigenvalues `β_n L` are the roots of `cos(βL)·cosh(βL) = 1`:
61```
624.73004074, 7.85320462, 10.99560784, 14.13716549, then (2n+1)π/2
63```
65The asymptotic form is exact to ~1e-9 from n = 5 up. With `c_L = √(E/ρ)`,
66`I = W h³/12`, `A = W h`:
68```
69f_n = κ · (β_n L)² / (2π L²) · √(E I / (ρ A))
70 = κ · (β_n L)² · c_L · h / (2π √12 · L²)
71```
73`κ = hollow_softening` accounts for stiffness removed by the slot.
75Pitch scales **linearly with h** and as **√(E/ρ)**. Verified: h = 6/10/16 mm
76gives 602/1003/1606 Hz; E = 6e9 vs 1.6e10 gives 614 vs 1003 Hz.
78Optional stretch (`γ = inharmonicity`), applied to partials above the first:
80```
81f_n ← f_1 · (f_n / f_1)^γ
82```
84Discard modes above `max_mode_hz` (20 kHz). Defaults keep 6.
86### 3.2 Mode shapes
88On `s = x/L ∈ [0,1]`, with `σ_n = (cosh βL − cos βL)/(sinh βL − sin βL)`:
90```
91φ̃_n(s) = cosh(β_n L s) + cos(β_n L s) − σ_n [ sinh(β_n L s) + sin(β_n L s) ]
92```
94Mass-normalise so the modal mass is 1:
96```
97φ_n(s) = φ̃_n(s) / ( √(∫₀¹ φ̃_n² ds) · √M_p )
98```
100`M_p = participating_mass_kg`. Units: `φ_n` is `kg^−1/2`, modal coordinate `q_n`
101is `m·√kg`, physical displacement `w = Σ q_n φ_n` is metres.
103> `M_p` is a **lumped correction, the least principled parameter in the model.**
104> The bar formula otherwise locks frequency and modal mass together, and the
105> thin-wall mass (`ρWhL` ≈ 93 g) makes the block far too compliant — 395 µm of
106> deflection and 119 dB SPL at 42 cm. A real woodblock is 300–500 g with the
107> whole box participating. `M_p = 0.35` gives 124 µm and 109.7 dB. It sets
108> absolute level, so treat SPL as plausible, not calibrated.
110### 3.3 Damping
112```
113ζ_n = ½ ( η₀ + η₁ ω_n )
114```
116`η₀ = loss_factor` (frequency-independent wood loss), `η₁ = loss_slope` (extra
117HF damping). Bigger = duller and shorter. `η₀`: 0.010 → T60 185 ms;
1180.030 → 66 ms.
120---
122## 4. Mallet contact (the only nonlinearity)
124Surface displacement and velocity at the strike point:
126```
127w = Σ_n q_n φ_n(s_c)
128ẇ = Σ_n q̇_n φ_n(s_c)
129δ = y − w (compression; y = mallet position)
130```
132**Hunt–Crossley**, with `λ = contact_damping` in **s/m**:
134```
135f_c = K δ^α (1 + λ (ẏ − ẇ)) if δ > 0, clamped to ≥ 0
136f_c = 0 otherwise
138m_h ÿ = −f_c
139```
141Initial conditions: `y = 0`, `ẏ = v₀ = mallet_velocity`, all else zero.
143> **Trap.** Do not write the loss term as `μ · d(δ^α)/dt`. That makes `μ` a
144> relaxation *time*; at μ = 0.35 s against a 0.1 ms collision, a 3 µm
145> compression produced **1.9 × 10⁷ N**. Every strike became a delta function and
146> mallet stiffness had literally zero effect on the spectrum. The form above
147> keeps the correction dimensionless and O(1).
149Expected behaviour, worth asserting in tests: contact duration scales roughly as
150`K^−2/5`. `K` = 3e7 / 3e8 / 3e9 / 3e10 → 0.92 / 0.59 / 0.28 / 0.16 ms and
151287 / 615 / 1289 / 3012 N. Because the contact is nonlinear, timbre varies with
152strike velocity on its own — that is physics, not a scripted velocity layer.
154---
156## 5. Cavity: Helmholtz resonator
158Effective plug length with end corrections on both ends:
160```
161L_eff = slit_length_m + 2 · 0.85 · √(S/π) S = slit_area_m2
162ω_H = c √( S / (V L_eff) ) V = cavity_volume_m3
163```
165Coupling coefficients — net volume each mode injects into the cavity, over the
166wall span `[lo, hi] = cavity_wall_span`:
168```
169A_n = χ · W L ∫_lo^hi φ_n(s) ds χ = cavity_coupling
170```
172Cavity pressure, from wall motion plus slit-plug displacement `ξ`:
174```
175p_c = (ρ₀ c² / V) ( S ξ + Σ_n A_n q_n )
176```
178Plug equation (mass per unit area `ρ₀ L_eff`, driven by `−p_c`):
180```
181ξ̈ = −2 ζ_H ω_H ξ̇ − ω_H² ξ − (c² / (V L_eff)) Σ_n A_n q_n
182```
184`A_n` appears both in the cavity drive and in the modal load below. **That
185symmetry is the correctness check** — if the coupling only goes one way, the
186energy will not be conserved.
188---
190## 6. Coupled system
192State vector, length `2N + 4`:
194```
195[ q₁..q_N , q̇₁..q̇_N , ξ , ξ̇ , y , ẏ ]
196```
198Derivatives:
200```
201q̈_n = −2 ζ_n ω_n q̇_n − ω_n² q_n + f_c φ_n(s_c) − A_n p_c
202ξ̈ = −2 ζ_H ω_H ξ̇ − ω_H² ξ − (c²/(V L_eff)) Σ_m A_m q_m
203ÿ = −f_c / m_h
204```
206Integrate with **RK4** at `f_hi = sample_rate × oversample`. Record `q̈_n(t)` for
207all modes and `ξ̈(t)` — those are the radiation inputs.
209The system is linear apart from `f_c`. Eigenvalues of the coupled linear part are
210all strictly decaying (slowest ≈ −34 s⁻¹ at defaults); a zero or positive real
211part means the coupling was assembled wrongly.
213**Convergence.** Relative to `oversample = 32`: `oversample = 8` differs by 5.2%,
214`16` by 1.8%. Use 16 for render, 8 is acceptable for interactive preview. The
215residual sits in the top mode near 18.7 kHz.
217---
219## 7. Radiation: pressure at a point
221Each cross-section of a bending bar **translates**, so the top and bottom faces
222inject equal and opposite volume. Every slice is a local dipole, and the net
223volume displacement of the bar is exactly zero.
225> **Trap.** This is why you cannot use a single lumped monopole term with
226> `B_n = ∫φ_n ds` for the body: elastic modes are orthogonal to rigid
227> translation, so that integral is **identically zero** and the body contributes
228> nothing. Near-field emphasis and HF directivity must come from path lengths and
229> delays, not from a gain constant.
231Discretise the top face into `M = radiating_elements` strips at
232`s_j = (j − ½)/M`, `x_j = (s_j − ½)L`, each of area `patch = W L / M`. Every
233strip contributes a pair of point sources:
235```
236r_j⁺ = (x_j, 0, +H/2) sign +1
237r_j⁻ = (x_j, 0, −H/2) sign −1
238```
240with volume acceleration
242```
243V̈_j(t) = ± patch · Σ_n φ_n(s_j) q̈_n(t)
244```
246Plus one monopole at the slit: `V̈_slit(t) = S ξ̈(t)`.
248Sum retarded monopoles over all `2M + 1` sources:
250```
251p(r, t) = Σ_j ρ₀ / (4π d_j) · V̈_j( t − d_j / c ) d_j = |r − r_j|
252```
254Integer-sample delays are fine: at `f_hi` = 768 kHz one sample is 0.45 mm.
256**Element count.** Spacing must satisfy `L/M < ¼ · c/f_max`; at defaults that
257needs M ≥ 44, so M = 48. Warn the user if `M` is too coarse, and if the listener
258is inside the source region (`|r| < L/2`) where point-source strips are crude.
260### Validation — the app should reproduce these
262| Listener (x, y, z) m | Peak | Note |
263|---|---|---|
264| 0.15, 0.25, 0.30 | 109.7 dB | default, oblique |
265| 0.00, 0.00, 0.40 | 100.9 dB | above centre; 2766 Hz partial −27 dB (symmetry null) |
266| 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 |
267| 0.00, 0.05, 0.02 | 125.4 dB | close to slit, cavity-dominated |
269Falloff is −7.8 dB per doubling close in, relaxing to −6.5 dB by 1.6 m: the
270dipole-to-inverse-square transition. If your port shows a flat −6 dB everywhere,
271the two faces are not in antiphase.
273---
275## 8. Sampling and output
2771. Solve at `f_hi = 48000 × oversample`.
2782. Radiate → `p(t)` in pascals at `f_hi`.
2793. **Anti-alias and decimate** by `oversample`. Python uses `scipy.resample_poly`;
280 in the browser implement a windowed-sinc (Kaiser, β ≈ 8) FIR lowpass at
281 ~0.45·48 kHz with roughly `12 × oversample` taps, then take every
282 `oversample`-th sample. Do not decimate without the filter — the double
283 differentiation puts real energy near `f_hi/2`.
2844. Prepend `pre_delay_ms` of silence, apply a short squared fade at the tail.
2855. Either normalise to `peak_dbfs`, or divide by `full_scale_pa` for absolute
286 pressure. Report peak Pa and dB SPL (`20 log₁₀(p/20 µPa)`), and warn on clip.
288> **Trap.** Do not "remove DC" by subtracting the mean of the first N samples.
289> The strike is at t = 0, so the head of the buffer is the loudest part of the
290> file; subtracting it injects a large offset across the whole signal. The model
291> has no DC — `p.mean()` is ~1e-11 — so no DC removal is needed at all.
293---
295## 9. Parameters
297| Parameter | Default | Effect |
298|---|---|---|
299| `length_m` | 0.200 | pitch ∝ 1/L² |
300| `width_m` | 0.055 | source area / level |
301| `wall_thickness_m` | 0.010 | **pitch ∝ h** |
302| `block_height_m` | 0.050 | dipole arm; bigger = more low end |
303| `density` | 850 | pitch ∝ 1/√ρ |
304| `youngs_modulus` | 1.6e10 | pitch ∝ √E |
305| `hollow_softening` | 0.90 | overall pitch trim |
306| `participating_mass_kg` | 0.35 | level ∝ 1/M (lumped; see §3.2) |
307| `inharmonicity` | 1.00 | >1 stretches partials (2766 → 3061 Hz at 1.10) |
308| `loss_factor` | 0.010 | T60 (185 ms → 66 ms at 0.030) |
309| `loss_slope` | 1.2e-7 | HF damping; main brightness control |
310| `strike_position` | 0.42 | partial balance. **Avoid 0.5** — a node of every antisymmetric mode |
311| `mallet_mass` | 0.025 | contact time, level |
312| `mallet_velocity` | 2.5 | dynamics *and* timbre (nonlinear contact) |
313| `mallet_stiffness` | 3.0e8 | hardness; centroid 1953 → 5550 Hz at 3e10 |
314| `contact_exponent` | 1.5 | 1.5 Hertzian, 2–3 felt |
315| `contact_damping` | 0.8 | Hunt–Crossley λ, **s/m not seconds** |
316| `cavity_volume_m3` | 2.6e-4 | Helmholtz pitch |
317| `slit_area_m2` | 8.0e-4 | Helmholtz pitch, cavity output |
318| `slit_length_m` | 0.012 | Helmholtz pitch |
319| `cavity_damping` | 0.06 | hollow-resonance Q |
320| `cavity_coupling` | 1.0 | 0 = dry bar, >1 = hollower |
321| `listener_*_m` | .15/.25/.30 | see §7 |
322| `radiating_elements` | 48 | accuracy, not timbre (see §7) |
323| `oversample` | 16 | accuracy vs speed |
325Defaults yield: partials 1003 / 2766 / 5423 / 8964 / 13390 / 18702 Hz (ratios
3261.00, 2.76, 5.40, 8.93 — bar inharmonicity, **not** octaves), Helmholtz 484 Hz,
327contact 0.30 ms at 794 N, T60 ≈ 149 ms, 109.7 dB SPL.
329---
331## 10. Browser implementation notes
333- **Cost.** Solve is `4 × T` derivative evaluations with N ≈ 6 (T = 345,600 at
334 defaults) — cheap. Radiation is the bottleneck: `2M+1` sources × N × T ≈ 200 M
335 multiply-adds. Compute `Σ_n φ_n(s_j) q̈_n(t)` **once per strip** and reuse it
336 for the `+z` and `−z` sources with opposite sign — that halves the work.
337- Run in a **Web Worker**, transfer the result as a `Float32Array`
338 `ArrayBuffer`, and copy into an `AudioBuffer` at 48 kHz. Keep the solver in
339 `Float64Array`.
340- For live parameter dragging, preview at `oversample = 8` and
341 `radiating_elements = 24`, then re-render at full quality on release.
342- Nothing needs real-time streaming — one strike is one offline render of well
343 under a second.
344- Cache the modal basis; it only depends on geometry, material,
345 `strike_position`, and `radiating_elements`, not on mallet or listener
346 parameters. Moving the listener requires only re-running §7.
347- Good things to visualise: mode shapes animated as `φ̃_n(s)` along the bar;
348 partial frequencies as a spectrum overlay; the contact-force pulse `f_c(t)`
349 (its width sets brightness); source strips and listener in the 3D view.