982850fWoodblock strike simulator: Python reference model and interactive web appJeremy Magland 1#!/usr/bin/env python3
2"""
3woodblock.py -- physically-based simulation of a single woodblock strike,
4sampled as acoustic pressure at a point in space.
6Model chain:
8 mallet (Hunt-Crossley contact)
9 -> body modes (free-free bar, modal reduction)
10 <-> air cavity (Helmholtz resonator, two-way coupled)
11 -> distributed surface sources + slit monopole
12 -> free-field pressure at the listener point
13 -> anti-aliased decimation -> .wav
15There is no microphone model: no directivity, no capsule response, no
16gain staging. The output is p(r, t) in pascals, sampled at one point,
17optionally normalized for convenience when writing the file.
19Everything interesting is in WoodblockConfig. Any field can be overridden
20from the command line, e.g.
22 python woodblock.py --listener_z_m 0.05 --strike_position 0.30
23 python woodblock.py --normalize 0 --full_scale_pa 50 -o absolute.wav
25Requires: numpy, scipy
26"""
28from __future__ import annotations
30import argparse
31import dataclasses
32import math
33import wave
34from dataclasses import dataclass
36import numpy as np
37from scipy.signal import resample_poly
40# ----------------------------------------------------------------------------
41# Configuration
42# ----------------------------------------------------------------------------
45class WoodblockConfig:
46 # --- Body geometry -------------------------------------------------------
47 length_m: float = 0.200 # long axis of the block
48 width_m: float = 0.055 # across the grain
49 wall_thickness_m: float = 0.010 # the hollowing means the *top wall* is the
50 # vibrating element, not the whole block.
51 # Pitch scales linearly with this.
52 block_height_m: float = 0.050 # outer height. Sets the spacing between the
53 # antiphase top/bottom faces, i.e. the
54 # dipole arm -- bigger = more low end.
56 # --- Wood material -------------------------------------------------------
57 density: float = 850.0 # kg/m^3 (rosewood ~850, maple ~700)
58 youngs_modulus: float = 1.6e10 # Pa, along the grain
59 # Frequencies scale as sqrt(E/rho), so stiffer/lighter = brighter + higher.
61 hollow_softening: float = 0.90 # 0.5-1.0; the cut slot removes stiffness
62 participating_mass_kg: float = 0.35 # modal mass used for normalization.
63 # The bar formula sets the *frequencies*,
64 # but a real box vibrates with far more
65 # mass than an isolated thin wall, and the
66 # two are otherwise locked together. This
67 # is a lumped stand-in for that; set 0 to
68 # fall back to the thin-wall beam mass.
69 # Scales amplitude as 1/mass.
70 inharmonicity: float = 1.00 # exponent on partial ratios. >1 stretches
71 # the upper partials (more "knock")
72 n_modes: int = 26
73 max_mode_hz: float = 20000.0 # discard modes above this
75 # --- Damping: zeta_n = 0.5 * (eta0 + eta1 * omega_n) ---------------------
76 loss_factor: float = 0.010 # eta0, frequency-independent wood loss
77 loss_slope: float = 1.2e-7 # eta1, extra damping for high partials
78 # (bigger = duller, faster decay)
80 # --- Mallet / strike -----------------------------------------------------
81 strike_position: float = 0.42 # 0..1 along the bar. Avoid exactly 0.5:
82 # dead centre is a node of every
83 # antisymmetric mode and sounds too pure
84 mallet_mass: float = 0.025 # kg
85 mallet_velocity: float = 2.5 # m/s at first contact -- dynamics
86 mallet_stiffness: float = 3.0e8 # K in f = K*delta^alpha. Bigger = harder,
87 # shorter contact, brighter click.
88 # ~1e7 soft rubber ... ~1e10 wood on wood
89 contact_exponent: float = 1.5 # 1.5 Hertzian (wood/plastic), 2-3 for felt
90 contact_damping: float = 0.8 # Hunt-Crossley lambda, in s/m (NOT a
91 # time). 0 = lossless contact; 0.5-2 is
92 # typical. Damps rebound, softens attack
94 # --- Air cavity (Helmholtz) ---------------------------------------------
95 cavity_volume_m3: float = 2.6e-4 # enclosed air volume
96 slit_area_m2: float = 8.0e-4 # opening cross-section
97 slit_length_m: float = 0.012 # wall thickness at the slot
98 cavity_damping: float = 0.06 # zeta_H; viscous + radiation losses
99 cavity_coupling: float = 1.0 # scales body <-> cavity coupling.
100 # 0 = dry bar, 1 = nominal, >1 = hollower
101 cavity_wall_span: tuple = (0.30, 0.70) # fraction of bar facing the cavity
103 # --- Acoustic field ------------------------------------------------------
104 # Block-centred coordinates, in metres: +x along the bar's long axis,
105 # +y across the width, +z up out of the struck face. The block spans
106 # x in [-length/2, +length/2] and its faces sit at z = +-block_height/2.
107 listener_x_m: float = 0.15
108 listener_y_m: float = 0.25
109 listener_z_m: float = 0.30
110 # The slit opening, treated as a monopole at this point.
111 slit_x_m: float = 0.0
112 slit_y_m: float = 0.0275
113 slit_z_m: float = 0.0
114 radiating_elements: int = 48 # surface discretisation along the bar.
115 # Must be fine enough that element
116 # spacing << wavelength (see check_setup)
117 speed_of_sound: float = 343.0
118 air_density: float = 1.204
120 # --- Rendering -----------------------------------------------------------
121 sample_rate: int = 48000
122 oversample: int = 16 # solver runs at oversample * sample_rate
123 duration_s: float = 0.45
124 pre_delay_ms: float = 3.0 # silence before the attack
125 normalize: bool = True # scale so the peak hits peak_dbfs
126 peak_dbfs: float = -3.0 # used only when normalize is True
127 full_scale_pa: float = 20.0 # when normalize is False, this many
128 # pascals maps to digital full scale
129 bit_depth: int = 24 # 16 or 24
131 def __post_init__(self):
132 assert 0.0 < self.strike_position < 1.0
133 assert self.bit_depth in (16, 24)
136 def listener(self) -> np.ndarray:
137 return np.array([self.listener_x_m, self.listener_y_m, self.listener_z_m])
140 def slit_point(self) -> np.ndarray:
141 return np.array([self.slit_x_m, self.slit_y_m, self.slit_z_m])
144# ----------------------------------------------------------------------------
145# Modal basis: free-free Euler-Bernoulli bar
146# ----------------------------------------------------------------------------
148# Roots of cos(bL)cosh(bL) = 1 for the free-free beam.
149_BETA_L = np.array([4.73004074, 7.85320462, 10.99560784, 14.13716549])
152def _beta_l(n: int) -> float:
153 """n is 1-based. Above the 4th root the asymptotic form is exact to 1e-9."""
154 if n <= len(_BETA_L):
155 return float(_BETA_L[n - 1])
156 return (2 * n + 1) * math.pi / 2.0
159def _mode_shape(bl: float, s: np.ndarray) -> np.ndarray:
160 """Free-free mode shape on s in [0,1], unnormalized."""
161 sigma = (math.cosh(bl) - math.cos(bl)) / (math.sinh(bl) - math.sin(bl))
162 x = bl * s
163 return (np.cosh(x) + np.cos(x)) - sigma * (np.sinh(x) + np.sin(x))
166class ModalBasis:
167 """Mass-normalized bending modes and their couplings.
169 phi_c[n] : mode shape at the strike point (1/sqrt(kg))
170 A[n] : volume displacement into the cavity (m^2/sqrt(kg))
171 phi_elem[n,j] : mode shape at surface element j (1/sqrt(kg))
172 """
174 def __init__(self, cfg: WoodblockConfig):
175 L, W, H = cfg.length_m, cfg.width_m, cfg.wall_thickness_m
176 area = W * H
177 mass = cfg.participating_mass_kg or (cfg.density * area * L)
178 c_long = math.sqrt(cfg.youngs_modulus / cfg.density)
180 # f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)), I = W H^3/12
181 prefactor = c_long * H / (2.0 * math.pi * math.sqrt(12.0) * L * L)
183 s = np.linspace(0.0, 1.0, 4096)
184 M = cfg.radiating_elements
185 s_elem = (np.arange(M) + 0.5) / M # element centres along the bar
187 freqs, phi_c, A, phi_elem = [], [], [], []
188 lo, hi = cfg.cavity_wall_span
189 cav_mask = (s >= lo) & (s <= hi)
191 for n in range(1, cfg.n_modes + 1):
192 bl = _beta_l(n)
193 f = prefactor * bl * bl * cfg.hollow_softening
194 if n > 1 and cfg.inharmonicity != 1.0:
195 f0 = prefactor * _beta_l(1) ** 2 * cfg.hollow_softening
196 f = f0 * (f / f0) ** cfg.inharmonicity
197 if f > cfg.max_mode_hz:
198 break
200 shape = _mode_shape(bl, s)
201 # Normalize so that (1/L) * integral(phi^2 dx) = 1, then divide by
202 # sqrt(total mass) to get the mass-normalized mode.
203 rms = math.sqrt(np.trapezoid(shape ** 2, s))
204 shape = shape / (rms * math.sqrt(mass))
206 freqs.append(f)
207 phi_c.append(float(np.interp(cfg.strike_position, s, shape)))
208 A.append(W * L * float(np.trapezoid(np.where(cav_mask, shape, 0.0), s)))
209 phi_elem.append(np.interp(s_elem, s, shape))
211 self.f = np.array(freqs)
212 self.omega = 2.0 * math.pi * self.f
213 self.phi_c = np.array(phi_c)
214 self.A = np.array(A) * cfg.cavity_coupling
215 self.phi_elem = np.array(phi_elem) # (N, M)
216 self.s_elem = s_elem
217 self.zeta = 0.5 * (cfg.loss_factor + cfg.loss_slope * self.omega)
218 self.n = len(freqs)
219 self.mass = mass
222# ----------------------------------------------------------------------------
223# Time-domain solver
224# ----------------------------------------------------------------------------
226_LAST_FC = [0.0] # scratch slot for contact-force diagnostics
229def simulate(cfg: WoodblockConfig):
230 """Integrate the coupled system.
232 Returns (modal accelerations (N,T), slit-plug acceleration (T,), basis, diag).
233 """
234 basis = ModalBasis(cfg)
235 N = basis.n
236 if N == 0:
237 raise ValueError("No modes below max_mode_hz -- check geometry/material.")
239 fs_hi = cfg.sample_rate * cfg.oversample
240 dt = 1.0 / fs_hi
241 n_steps = int(cfg.duration_s * fs_hi)
243 om2 = basis.omega ** 2
244 two_zeta_om = 2.0 * basis.zeta * basis.omega
245 phi_c, A = basis.phi_c, basis.A
247 rho0, c = cfg.air_density, cfg.speed_of_sound
248 V, S, Lslit = cfg.cavity_volume_m3, cfg.slit_area_m2, cfg.slit_length_m
249 # End correction on the slit: ~0.85 * equivalent radius per open end.
250 Leff = Lslit + 2.0 * 0.85 * math.sqrt(S / math.pi)
251 omega_H = c * math.sqrt(S / (V * Leff))
252 two_zeta_H_om = 2.0 * cfg.cavity_damping * omega_H
253 p_scale = rho0 * c ** 2 / V
254 cav_drive = c ** 2 / (V * Leff)
256 K = cfg.mallet_stiffness
257 alpha = cfg.contact_exponent
258 mu = cfg.contact_damping
259 inv_mh = 1.0 / cfg.mallet_mass
261 # State: [q(N), qdot(N), xi, xidot, y, ydot]
262 IX, IXD, IY, IYD = 2 * N, 2 * N + 1, 2 * N + 2, 2 * N + 3
263 state = np.zeros(2 * N + 4)
264 state[IYD] = cfg.mallet_velocity
266 def deriv(s):
267 q, qd = s[:N], s[N:2 * N]
268 xi, xid, y, yd = s[IX], s[IXD], s[IY], s[IYD]
270 # Mallet contact
271 w = float(phi_c @ q)
272 wd = float(phi_c @ qd)
273 delta = y - w
274 if delta > 0.0:
275 # Hunt-Crossley: f = K d^a (1 + lambda * ddot). lambda has units of
276 # s/m, so the correction is dimensionless and stays O(1) over a
277 # sub-millisecond collision. (Writing the loss as mu*d(d^a)/dt
278 # instead makes mu a relaxation *time*; any plausible-looking value
279 # then dwarfs the elastic term and fires one huge impulse.)
280 fc = K * delta ** alpha * (1.0 + mu * (yd - wd))
281 if fc < 0.0:
282 fc = 0.0
283 else:
284 fc = 0.0
285 _LAST_FC[0] = fc
287 # Cavity pressure from combined volume change of walls and slit plug
288 vol = S * xi + float(A @ q)
289 pc = p_scale * vol
291 qdd = -two_zeta_om * qd - om2 * q + fc * phi_c - A * pc
292 xidd = -two_zeta_H_om * xid - omega_H ** 2 * xi - cav_drive * float(A @ q)
293 ydd = -fc * inv_mh
295 d = np.empty_like(s)
296 d[:N] = qd
297 d[N:2 * N] = qdd
298 d[IX] = xid
299 d[IXD] = xidd
300 d[IY] = yd
301 d[IYD] = ydd
302 return d, qdd, xidd
304 Qdd = np.empty((N, n_steps))
305 Xdd = np.empty(n_steps)
306 force = np.zeros(n_steps)
308 for k in range(n_steps):
309 k1, qdd, xidd = deriv(state)
310 Qdd[:, k] = qdd
311 Xdd[k] = xidd
312 force[k] = _LAST_FC[0]
314 k2, _, _ = deriv(state + 0.5 * dt * k1)
315 k3, _, _ = deriv(state + 0.5 * dt * k2)
316 k4, _, _ = deriv(state + dt * k3)
317 state += (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
319 contact = force > 0.0
320 diag = {
321 "contact_force": force,
322 "peak_force_N": float(force.max()),
323 "contact_ms": float(contact.sum() * dt * 1e3),
324 "helmholtz_hz": omega_H / (2.0 * math.pi),
325 }
326 return Qdd, Xdd, basis, diag
329# ----------------------------------------------------------------------------
330# Radiation: free-field pressure at a point
331# ----------------------------------------------------------------------------
333def radiate(cfg: WoodblockConfig, Qdd: np.ndarray, Xdd: np.ndarray,
334 basis: ModalBasis, fs_hi: float):
335 """Sum retarded monopole contributions from every source to the listener.
337 p(r,t) = sum_j rho0 / (4 pi d_j) * Qddot_j(t - d_j/c)
339 The struck face and the opposite face move together (each cross-section
340 translates), so they inject equal and opposite volume -- each slice of the
341 bar is a local dipole. Summed over the bar the net volume displacement is
342 exactly zero, which is the correct statement that a bending bar has no
343 monopole component. The audible near-field asymmetry comes entirely from
344 the different path lengths and delays, not from any fudge factor.
345 """
346 rho0, c = cfg.air_density, cfg.speed_of_sound
347 T = Qdd.shape[1]
348 L, W = cfg.length_m, cfg.width_m
349 half_h = 0.5 * cfg.block_height_m
350 M = cfg.radiating_elements
351 patch = W * L / M # area represented by one element
353 # Surface element positions along the bar, on both faces
354 x_elem = (basis.s_elem - 0.5) * L
355 src_pos, src_sign, src_mode = [], [], []
356 for j in range(M):
357 src_pos.append([x_elem[j], 0.0, +half_h]); src_sign.append(+1.0); src_mode.append(j)
358 src_pos.append([x_elem[j], 0.0, -half_h]); src_sign.append(-1.0); src_mode.append(j)
359 src_pos = np.array(src_pos)
361 r = cfg.listener
362 dists = np.linalg.norm(src_pos - r, axis=1)
363 slit_dist = float(np.linalg.norm(cfg.slit_point - r))
365 max_delay = int(math.ceil(max(dists.max(), slit_dist) / c * fs_hi)) + 2
366 out = np.zeros(T + max_delay)
368 # Surface (body) sources -- volume acceleration = patch area * wddot
369 for idx in range(len(dists)):
370 j = src_mode[idx]
371 wdd = basis.phi_elem[:, j] @ Qdd # (T,) surface acceleration
372 amp = src_sign[idx] * patch * rho0 / (4.0 * math.pi * dists[idx])
373 k = int(round(dists[idx] / c * fs_hi))
374 out[k:k + T] += amp * wdd
376 # Slit monopole -- volume acceleration = slit area * plug acceleration
377 amp = cfg.slit_area_m2 * rho0 / (4.0 * math.pi * slit_dist)
378 k = int(round(slit_dist / c * fs_hi))
379 out[k:k + T] += amp * Xdd
381 return out[:T]
384def check_setup(cfg: WoodblockConfig, basis: ModalBasis) -> list:
385 """Warn about configurations where the radiation model stops being valid."""
386 msgs = []
387 c = cfg.speed_of_sound
388 spacing = cfg.length_m / cfg.radiating_elements
389 f_max = basis.f[-1]
390 if spacing > 0.25 * c / f_max:
391 need = int(math.ceil(cfg.length_m / (0.25 * c / f_max)))
392 msgs.append(f"element spacing is coarse at {f_max:.0f} Hz; "
393 f"consider --radiating_elements {need}")
394 d = float(np.linalg.norm(cfg.listener))
395 if d < 0.5 * cfg.length_m:
396 msgs.append(f"listener is {d*100:.0f} cm from the block centre, inside "
397 f"the source region; point-source elements are crude here")
398 return msgs
401# ----------------------------------------------------------------------------
402# Rendering
403# ----------------------------------------------------------------------------
405def render(cfg: WoodblockConfig):
406 Qdd, Xdd, basis, diag = simulate(cfg)
407 fs_hi = cfg.sample_rate * cfg.oversample
409 p = radiate(cfg, Qdd, Xdd, basis, fs_hi) # pascals at the listener
410 diag["peak_pa"] = float(np.abs(p).max())
411 diag["peak_spl_db"] = 20.0 * math.log10(diag["peak_pa"] / 20e-6 + 1e-30)
412 diag["warnings"] = check_setup(cfg, basis)
414 # Anti-aliased decimation to the target rate
415 y = resample_poly(p, 1, cfg.oversample)
417 # Pre-roll so the attack transient does not sit against sample 0. (Do NOT
418 # try to de-DC by subtracting the head of the buffer -- the strike is at
419 # t=0, so the head is the loudest part of the file.)
420 pad = int(cfg.pre_delay_ms * 1e-3 * cfg.sample_rate)
421 if pad > 0:
422 y = np.concatenate([np.zeros(pad), y])[:int(cfg.duration_s * cfg.sample_rate)]
424 fade = min(256, y.size // 8)
425 if fade > 1:
426 y[-fade:] *= np.linspace(1.0, 0.0, fade) ** 2
428 if cfg.normalize:
429 peak = np.max(np.abs(y))
430 if peak > 0:
431 y = y / peak * (10.0 ** (cfg.peak_dbfs / 20.0))
432 else:
433 y = y / cfg.full_scale_pa
434 if np.max(np.abs(y)) > 1.0:
435 diag["warnings"].append(
436 f"clipping: peak is {diag['peak_pa']:.2f} Pa but full_scale_pa "
437 f"is {cfg.full_scale_pa:.2f}")
438 return y, basis, diag
441def write_wav(path: str, y: np.ndarray, cfg: WoodblockConfig) -> None:
442 y = np.clip(y, -1.0, 1.0)
443 with wave.open(path, "wb") as w:
444 w.setnchannels(1)
445 w.setsampwidth(cfg.bit_depth // 8)
446 w.setframerate(cfg.sample_rate)
447 if cfg.bit_depth == 16:
448 w.writeframes((y * 32767.0).astype("<i2").tobytes())
449 else:
450 i32 = (y * 8388607.0).astype("<i4")
451 raw = np.frombuffer(i32.tobytes(), dtype=np.uint8).reshape(-1, 4)
452 w.writeframes(raw[:, :3].tobytes())
455# ----------------------------------------------------------------------------
456# CLI
457# ----------------------------------------------------------------------------
459def build_parser() -> argparse.ArgumentParser:
460 p = argparse.ArgumentParser(
461 description="Simulate a woodblock strike and sample the pressure "
462 "at a point in space.")
463 p.add_argument("-o", "--output", default="woodblock.wav")
464 defaults = WoodblockConfig()
465 for f in dataclasses.fields(defaults):
466 val = getattr(defaults, f.name)
467 if isinstance(val, tuple):
468 continue
469 cast = (lambda s: bool(int(s))) if isinstance(val, bool) else type(val)
470 p.add_argument(f"--{f.name}", type=cast, default=None,
471 help=f"(default: {val})")
472 return p
475def main() -> None:
476 args = build_parser().parse_args()
477 overrides = {k: v for k, v in vars(args).items()
478 if k != "output" and v is not None}
479 cfg = WoodblockConfig(**overrides)
481 y, basis, diag = render(cfg)
482 write_wav(args.output, y, cfg)
484 d = float(np.linalg.norm(cfg.listener))
485 print(f"modes kept : {basis.n}")
486 print("partials Hz: " +
487 ", ".join(f"{f:.0f}" for f in basis.f[:6]) +
488 ("..." if basis.n > 6 else ""))
489 print("ratios : " +
490 ", ".join(f"{f / basis.f[0]:.2f}" for f in basis.f[:6]))
491 print(f"helmholtz : {diag['helmholtz_hz']:.0f} Hz")
492 print(f"contact : {diag['contact_ms']:.3f} ms, peak {diag['peak_force_N']:.0f} N")
493 print(f"listener : {d * 100:.0f} cm from block centre")
494 print(f"peak press.: {diag['peak_pa']:.3f} Pa ({diag['peak_spl_db']:.1f} dB SPL)")
495 print(f"length : {y.size / cfg.sample_rate:.3f} s @ {cfg.sample_rate} Hz")
496 for w in diag["warnings"]:
497 print(f"WARNING : {w}")
498 print(f"wrote : {args.output}")
501if __name__ == "__main__":
502 main()