#!/usr/bin/env python3 """ woodblock.py -- physically-based simulation of a single woodblock strike, sampled as acoustic pressure at a point in space. Model chain: mallet (Hunt-Crossley contact) -> body modes (free-free bar, modal reduction) <-> air cavity (Helmholtz resonator, two-way coupled) -> distributed surface sources + slit monopole -> free-field pressure at the listener point -> anti-aliased decimation -> .wav There is no microphone model: no directivity, no capsule response, no gain staging. The output is p(r, t) in pascals, sampled at one point, optionally normalized for convenience when writing the file. Everything interesting is in WoodblockConfig. Any field can be overridden from the command line, e.g. python woodblock.py --listener_z_m 0.05 --strike_position 0.30 python woodblock.py --normalize 0 --full_scale_pa 50 -o absolute.wav Requires: numpy, scipy """ from __future__ import annotations import argparse import dataclasses import math import wave from dataclasses import dataclass import numpy as np from scipy.signal import resample_poly # ---------------------------------------------------------------------------- # Configuration # ---------------------------------------------------------------------------- @dataclass class WoodblockConfig: # --- Body geometry ------------------------------------------------------- length_m: float = 0.200 # long axis of the block width_m: float = 0.055 # across the grain wall_thickness_m: float = 0.010 # the hollowing means the *top wall* is the # vibrating element, not the whole block. # Pitch scales linearly with this. block_height_m: float = 0.050 # outer height. Sets the spacing between the # antiphase top/bottom faces, i.e. the # dipole arm -- bigger = more low end. # --- Wood material ------------------------------------------------------- density: float = 850.0 # kg/m^3 (rosewood ~850, maple ~700) youngs_modulus: float = 1.6e10 # Pa, along the grain # Frequencies scale as sqrt(E/rho), so stiffer/lighter = brighter + higher. hollow_softening: float = 0.90 # 0.5-1.0; the cut slot removes stiffness participating_mass_kg: float = 0.35 # modal mass used for normalization. # The bar formula sets the *frequencies*, # but a real box vibrates with far more # mass than an isolated thin wall, and the # two are otherwise locked together. This # is a lumped stand-in for that; set 0 to # fall back to the thin-wall beam mass. # Scales amplitude as 1/mass. inharmonicity: float = 1.00 # exponent on partial ratios. >1 stretches # the upper partials (more "knock") n_modes: int = 26 max_mode_hz: float = 20000.0 # discard modes above this # --- Damping: zeta_n = 0.5 * (eta0 + eta1 * omega_n) --------------------- loss_factor: float = 0.010 # eta0, frequency-independent wood loss loss_slope: float = 1.2e-7 # eta1, extra damping for high partials # (bigger = duller, faster decay) # --- Mallet / strike ----------------------------------------------------- strike_position: float = 0.42 # 0..1 along the bar. Avoid exactly 0.5: # dead centre is a node of every # antisymmetric mode and sounds too pure mallet_mass: float = 0.025 # kg mallet_velocity: float = 2.5 # m/s at first contact -- dynamics mallet_stiffness: float = 3.0e8 # K in f = K*delta^alpha. Bigger = harder, # shorter contact, brighter click. # ~1e7 soft rubber ... ~1e10 wood on wood contact_exponent: float = 1.5 # 1.5 Hertzian (wood/plastic), 2-3 for felt contact_damping: float = 0.8 # Hunt-Crossley lambda, in s/m (NOT a # time). 0 = lossless contact; 0.5-2 is # typical. Damps rebound, softens attack # --- Air cavity (Helmholtz) --------------------------------------------- cavity_volume_m3: float = 2.6e-4 # enclosed air volume slit_area_m2: float = 8.0e-4 # opening cross-section slit_length_m: float = 0.012 # wall thickness at the slot cavity_damping: float = 0.06 # zeta_H; viscous + radiation losses cavity_coupling: float = 1.0 # scales body <-> cavity coupling. # 0 = dry bar, 1 = nominal, >1 = hollower cavity_wall_span: tuple = (0.30, 0.70) # fraction of bar facing the cavity # --- Acoustic field ------------------------------------------------------ # Block-centred coordinates, in metres: +x along the bar's long axis, # +y across the width, +z up out of the struck face. The block spans # x in [-length/2, +length/2] and its faces sit at z = +-block_height/2. listener_x_m: float = 0.15 listener_y_m: float = 0.25 listener_z_m: float = 0.30 # The slit opening, treated as a monopole at this point. slit_x_m: float = 0.0 slit_y_m: float = 0.0275 slit_z_m: float = 0.0 radiating_elements: int = 48 # surface discretisation along the bar. # Must be fine enough that element # spacing << wavelength (see check_setup) speed_of_sound: float = 343.0 air_density: float = 1.204 # --- Rendering ----------------------------------------------------------- sample_rate: int = 48000 oversample: int = 16 # solver runs at oversample * sample_rate duration_s: float = 0.45 pre_delay_ms: float = 3.0 # silence before the attack normalize: bool = True # scale so the peak hits peak_dbfs peak_dbfs: float = -3.0 # used only when normalize is True full_scale_pa: float = 20.0 # when normalize is False, this many # pascals maps to digital full scale bit_depth: int = 24 # 16 or 24 def __post_init__(self): assert 0.0 < self.strike_position < 1.0 assert self.bit_depth in (16, 24) @property def listener(self) -> np.ndarray: return np.array([self.listener_x_m, self.listener_y_m, self.listener_z_m]) @property def slit_point(self) -> np.ndarray: return np.array([self.slit_x_m, self.slit_y_m, self.slit_z_m]) # ---------------------------------------------------------------------------- # Modal basis: free-free Euler-Bernoulli bar # ---------------------------------------------------------------------------- # Roots of cos(bL)cosh(bL) = 1 for the free-free beam. _BETA_L = np.array([4.73004074, 7.85320462, 10.99560784, 14.13716549]) def _beta_l(n: int) -> float: """n is 1-based. Above the 4th root the asymptotic form is exact to 1e-9.""" if n <= len(_BETA_L): return float(_BETA_L[n - 1]) return (2 * n + 1) * math.pi / 2.0 def _mode_shape(bl: float, s: np.ndarray) -> np.ndarray: """Free-free mode shape on s in [0,1], unnormalized.""" sigma = (math.cosh(bl) - math.cos(bl)) / (math.sinh(bl) - math.sin(bl)) x = bl * s return (np.cosh(x) + np.cos(x)) - sigma * (np.sinh(x) + np.sin(x)) class ModalBasis: """Mass-normalized bending modes and their couplings. phi_c[n] : mode shape at the strike point (1/sqrt(kg)) A[n] : volume displacement into the cavity (m^2/sqrt(kg)) phi_elem[n,j] : mode shape at surface element j (1/sqrt(kg)) """ def __init__(self, cfg: WoodblockConfig): L, W, H = cfg.length_m, cfg.width_m, cfg.wall_thickness_m area = W * H mass = cfg.participating_mass_kg or (cfg.density * area * L) c_long = math.sqrt(cfg.youngs_modulus / cfg.density) # f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)), I = W H^3/12 prefactor = c_long * H / (2.0 * math.pi * math.sqrt(12.0) * L * L) s = np.linspace(0.0, 1.0, 4096) M = cfg.radiating_elements s_elem = (np.arange(M) + 0.5) / M # element centres along the bar freqs, phi_c, A, phi_elem = [], [], [], [] lo, hi = cfg.cavity_wall_span cav_mask = (s >= lo) & (s <= hi) for n in range(1, cfg.n_modes + 1): bl = _beta_l(n) f = prefactor * bl * bl * cfg.hollow_softening if n > 1 and cfg.inharmonicity != 1.0: f0 = prefactor * _beta_l(1) ** 2 * cfg.hollow_softening f = f0 * (f / f0) ** cfg.inharmonicity if f > cfg.max_mode_hz: break shape = _mode_shape(bl, s) # Normalize so that (1/L) * integral(phi^2 dx) = 1, then divide by # sqrt(total mass) to get the mass-normalized mode. rms = math.sqrt(np.trapezoid(shape ** 2, s)) shape = shape / (rms * math.sqrt(mass)) freqs.append(f) phi_c.append(float(np.interp(cfg.strike_position, s, shape))) A.append(W * L * float(np.trapezoid(np.where(cav_mask, shape, 0.0), s))) phi_elem.append(np.interp(s_elem, s, shape)) self.f = np.array(freqs) self.omega = 2.0 * math.pi * self.f self.phi_c = np.array(phi_c) self.A = np.array(A) * cfg.cavity_coupling self.phi_elem = np.array(phi_elem) # (N, M) self.s_elem = s_elem self.zeta = 0.5 * (cfg.loss_factor + cfg.loss_slope * self.omega) self.n = len(freqs) self.mass = mass # ---------------------------------------------------------------------------- # Time-domain solver # ---------------------------------------------------------------------------- _LAST_FC = [0.0] # scratch slot for contact-force diagnostics def simulate(cfg: WoodblockConfig): """Integrate the coupled system. Returns (modal accelerations (N,T), slit-plug acceleration (T,), basis, diag). """ basis = ModalBasis(cfg) N = basis.n if N == 0: raise ValueError("No modes below max_mode_hz -- check geometry/material.") fs_hi = cfg.sample_rate * cfg.oversample dt = 1.0 / fs_hi n_steps = int(cfg.duration_s * fs_hi) om2 = basis.omega ** 2 two_zeta_om = 2.0 * basis.zeta * basis.omega phi_c, A = basis.phi_c, basis.A rho0, c = cfg.air_density, cfg.speed_of_sound V, S, Lslit = cfg.cavity_volume_m3, cfg.slit_area_m2, cfg.slit_length_m # End correction on the slit: ~0.85 * equivalent radius per open end. Leff = Lslit + 2.0 * 0.85 * math.sqrt(S / math.pi) omega_H = c * math.sqrt(S / (V * Leff)) two_zeta_H_om = 2.0 * cfg.cavity_damping * omega_H p_scale = rho0 * c ** 2 / V cav_drive = c ** 2 / (V * Leff) K = cfg.mallet_stiffness alpha = cfg.contact_exponent mu = cfg.contact_damping inv_mh = 1.0 / cfg.mallet_mass # State: [q(N), qdot(N), xi, xidot, y, ydot] IX, IXD, IY, IYD = 2 * N, 2 * N + 1, 2 * N + 2, 2 * N + 3 state = np.zeros(2 * N + 4) state[IYD] = cfg.mallet_velocity def deriv(s): q, qd = s[:N], s[N:2 * N] xi, xid, y, yd = s[IX], s[IXD], s[IY], s[IYD] # Mallet contact w = float(phi_c @ q) wd = float(phi_c @ qd) delta = y - w if delta > 0.0: # Hunt-Crossley: f = K d^a (1 + lambda * ddot). lambda has units of # s/m, so the correction is dimensionless and stays O(1) over a # sub-millisecond collision. (Writing the loss as mu*d(d^a)/dt # instead makes mu a relaxation *time*; any plausible-looking value # then dwarfs the elastic term and fires one huge impulse.) fc = K * delta ** alpha * (1.0 + mu * (yd - wd)) if fc < 0.0: fc = 0.0 else: fc = 0.0 _LAST_FC[0] = fc # Cavity pressure from combined volume change of walls and slit plug vol = S * xi + float(A @ q) pc = p_scale * vol qdd = -two_zeta_om * qd - om2 * q + fc * phi_c - A * pc xidd = -two_zeta_H_om * xid - omega_H ** 2 * xi - cav_drive * float(A @ q) ydd = -fc * inv_mh d = np.empty_like(s) d[:N] = qd d[N:2 * N] = qdd d[IX] = xid d[IXD] = xidd d[IY] = yd d[IYD] = ydd return d, qdd, xidd Qdd = np.empty((N, n_steps)) Xdd = np.empty(n_steps) force = np.zeros(n_steps) for k in range(n_steps): k1, qdd, xidd = deriv(state) Qdd[:, k] = qdd Xdd[k] = xidd force[k] = _LAST_FC[0] k2, _, _ = deriv(state + 0.5 * dt * k1) k3, _, _ = deriv(state + 0.5 * dt * k2) k4, _, _ = deriv(state + dt * k3) state += (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4) contact = force > 0.0 diag = { "contact_force": force, "peak_force_N": float(force.max()), "contact_ms": float(contact.sum() * dt * 1e3), "helmholtz_hz": omega_H / (2.0 * math.pi), } return Qdd, Xdd, basis, diag # ---------------------------------------------------------------------------- # Radiation: free-field pressure at a point # ---------------------------------------------------------------------------- def radiate(cfg: WoodblockConfig, Qdd: np.ndarray, Xdd: np.ndarray, basis: ModalBasis, fs_hi: float): """Sum retarded monopole contributions from every source to the listener. p(r,t) = sum_j rho0 / (4 pi d_j) * Qddot_j(t - d_j/c) The struck face and the opposite face move together (each cross-section translates), so they inject equal and opposite volume -- each slice of the bar is a local dipole. Summed over the bar the net volume displacement is exactly zero, which is the correct statement that a bending bar has no monopole component. The audible near-field asymmetry comes entirely from the different path lengths and delays, not from any fudge factor. """ rho0, c = cfg.air_density, cfg.speed_of_sound T = Qdd.shape[1] L, W = cfg.length_m, cfg.width_m half_h = 0.5 * cfg.block_height_m M = cfg.radiating_elements patch = W * L / M # area represented by one element # Surface element positions along the bar, on both faces x_elem = (basis.s_elem - 0.5) * L src_pos, src_sign, src_mode = [], [], [] for j in range(M): src_pos.append([x_elem[j], 0.0, +half_h]); src_sign.append(+1.0); src_mode.append(j) src_pos.append([x_elem[j], 0.0, -half_h]); src_sign.append(-1.0); src_mode.append(j) src_pos = np.array(src_pos) r = cfg.listener dists = np.linalg.norm(src_pos - r, axis=1) slit_dist = float(np.linalg.norm(cfg.slit_point - r)) max_delay = int(math.ceil(max(dists.max(), slit_dist) / c * fs_hi)) + 2 out = np.zeros(T + max_delay) # Surface (body) sources -- volume acceleration = patch area * wddot for idx in range(len(dists)): j = src_mode[idx] wdd = basis.phi_elem[:, j] @ Qdd # (T,) surface acceleration amp = src_sign[idx] * patch * rho0 / (4.0 * math.pi * dists[idx]) k = int(round(dists[idx] / c * fs_hi)) out[k:k + T] += amp * wdd # Slit monopole -- volume acceleration = slit area * plug acceleration amp = cfg.slit_area_m2 * rho0 / (4.0 * math.pi * slit_dist) k = int(round(slit_dist / c * fs_hi)) out[k:k + T] += amp * Xdd return out[:T] def check_setup(cfg: WoodblockConfig, basis: ModalBasis) -> list: """Warn about configurations where the radiation model stops being valid.""" msgs = [] c = cfg.speed_of_sound spacing = cfg.length_m / cfg.radiating_elements f_max = basis.f[-1] if spacing > 0.25 * c / f_max: need = int(math.ceil(cfg.length_m / (0.25 * c / f_max))) msgs.append(f"element spacing is coarse at {f_max:.0f} Hz; " f"consider --radiating_elements {need}") d = float(np.linalg.norm(cfg.listener)) if d < 0.5 * cfg.length_m: msgs.append(f"listener is {d*100:.0f} cm from the block centre, inside " f"the source region; point-source elements are crude here") return msgs # ---------------------------------------------------------------------------- # Rendering # ---------------------------------------------------------------------------- def render(cfg: WoodblockConfig): Qdd, Xdd, basis, diag = simulate(cfg) fs_hi = cfg.sample_rate * cfg.oversample p = radiate(cfg, Qdd, Xdd, basis, fs_hi) # pascals at the listener diag["peak_pa"] = float(np.abs(p).max()) diag["peak_spl_db"] = 20.0 * math.log10(diag["peak_pa"] / 20e-6 + 1e-30) diag["warnings"] = check_setup(cfg, basis) # Anti-aliased decimation to the target rate y = resample_poly(p, 1, cfg.oversample) # Pre-roll so the attack transient does not sit against sample 0. (Do NOT # try to de-DC by subtracting the head of the buffer -- the strike is at # t=0, so the head is the loudest part of the file.) pad = int(cfg.pre_delay_ms * 1e-3 * cfg.sample_rate) if pad > 0: y = np.concatenate([np.zeros(pad), y])[:int(cfg.duration_s * cfg.sample_rate)] fade = min(256, y.size // 8) if fade > 1: y[-fade:] *= np.linspace(1.0, 0.0, fade) ** 2 if cfg.normalize: peak = np.max(np.abs(y)) if peak > 0: y = y / peak * (10.0 ** (cfg.peak_dbfs / 20.0)) else: y = y / cfg.full_scale_pa if np.max(np.abs(y)) > 1.0: diag["warnings"].append( f"clipping: peak is {diag['peak_pa']:.2f} Pa but full_scale_pa " f"is {cfg.full_scale_pa:.2f}") return y, basis, diag def write_wav(path: str, y: np.ndarray, cfg: WoodblockConfig) -> None: y = np.clip(y, -1.0, 1.0) with wave.open(path, "wb") as w: w.setnchannels(1) w.setsampwidth(cfg.bit_depth // 8) w.setframerate(cfg.sample_rate) if cfg.bit_depth == 16: w.writeframes((y * 32767.0).astype(" argparse.ArgumentParser: p = argparse.ArgumentParser( description="Simulate a woodblock strike and sample the pressure " "at a point in space.") p.add_argument("-o", "--output", default="woodblock.wav") defaults = WoodblockConfig() for f in dataclasses.fields(defaults): val = getattr(defaults, f.name) if isinstance(val, tuple): continue cast = (lambda s: bool(int(s))) if isinstance(val, bool) else type(val) p.add_argument(f"--{f.name}", type=cast, default=None, help=f"(default: {val})") return p def main() -> None: args = build_parser().parse_args() overrides = {k: v for k, v in vars(args).items() if k != "output" and v is not None} cfg = WoodblockConfig(**overrides) y, basis, diag = render(cfg) write_wav(args.output, y, cfg) d = float(np.linalg.norm(cfg.listener)) print(f"modes kept : {basis.n}") print("partials Hz: " + ", ".join(f"{f:.0f}" for f in basis.f[:6]) + ("..." if basis.n > 6 else "")) print("ratios : " + ", ".join(f"{f / basis.f[0]:.2f}" for f in basis.f[:6])) print(f"helmholtz : {diag['helmholtz_hz']:.0f} Hz") print(f"contact : {diag['contact_ms']:.3f} ms, peak {diag['peak_force_N']:.0f} N") print(f"listener : {d * 100:.0f} cm from block centre") print(f"peak press.: {diag['peak_pa']:.3f} Pa ({diag['peak_spl_db']:.1f} dB SPL)") print(f"length : {y.size / cfg.sample_rate:.3f} s @ {cfg.sample_rate} Hz") for w in diag["warnings"]: print(f"WARNING : {w}") print(f"wrote : {args.output}") if __name__ == "__main__": main()