3688ceacymsim: Python port of cymbalsim with a modal linear model and a PyTorch von Karman modelJeremy Magland 1"""Modal synthesis of a strike. Each mode m of the linearized plate is a damped
2oscillator q'' + 2 sigma q' + omega^2 q = f, with sigma the sum of its three
3loss rates. The stick is the Hertz contact of the full model, stepped against
4the modal system with an exact update of every oscillator under a piecewise
5constant force. Once the stick has left, every mode is a damped sinusoid with a
6known complex amplitude, and the delayed Rayleigh sum over the plate collapses
7to one complex gain per mode. During the attack, while the modes are being
8driven, the delayed sum is evaluated directly from a table of the nodal
9accelerations.
10"""
12from __future__ import annotations
14import math
15from dataclasses import dataclass
17import numpy as np
19from ..constants import STICK_RELEASE
20from ..params import Listener, Strike
21from ..plate import Plate
22from ..sound import Sound
23from .modes import Modes
26def strike_gains(modes: Modes, strike: Strike) -> np.ndarray:
27 """The contact: modal gains g_m = sum over the patch of (m_k / M_patch) phi_m(k),
28 so that the modal force is -F g_m and the contact coordinate is sum_m g_m q_m."""
29 pl = modes.plate
30 patch = pl.strike_patch(strike)
31 frac = pl.node_mass[patch.mask] / patch.mass
32 return modes.w[:, patch.mask] @ frac
36class Contact:
37 """The state of every mode at the end of the contact phase and the tables recorded during it."""
39 t_end: float
40 """Time at which the stick is withdrawn, STICK_RELEASE after the contact last opened [s]."""
41 t_contact: float
42 """Time at which the contact last opened [s]; every mode decays freely from here on."""
43 amplitude: np.ndarray
44 """Complex modal amplitudes A at ``t_end``: q(t) = Re[A exp(s (t - t_end))], s = -sigma + i omega_d."""
45 qdd: np.ndarray
46 """Modal accelerations sampled at ``table_rate`` from t = 0 to t_end, ``(table_len, M)``."""
47 table_rate: float
48 q_frames: np.ndarray
49 """Modal displacements at the movie frame times before ``t_end``, ``(n_frames, M)``."""
50 force: np.ndarray
51 """Stick force history at the fine step [N]."""
52 dt: float
53 peak_force: float
54 contact_time: float
55 """Total time the contact was closed [s]."""
58def simulate_contact(modes: Modes, sigma: np.ndarray, g: np.ndarray, strike: Strike, fps: float, fs: float, dt: float = 2e-6, table_rate: float | None = None) -> Contact:
59 """Steps the stick against the modal system until the contact has been open
60 for STICK_RELEASE, recording the tables the synthesis needs. Every mode is
61 advanced by its exact state transition under a constant force over the step;
62 the stick by symplectic Euler; the force is the Hertz law of the penetration
63 at the start of the step."""
64 table_rate = table_rate or 2 * fs
65 omega = modes.omega
66 wd = np.sqrt(omega**2 - sigma**2)
67 e, c, s = np.exp(-sigma * dt), np.cos(wd * dt), np.sin(wd * dt)
68 p11 = e * (c + sigma * s / wd)
69 p12 = e * s / wd
70 p21 = -e * omega**2 * s / wd
71 p22 = e * (c - sigma * s / wd)
72 M = modes.count
73 q, qd = np.zeros(M), np.zeros(M)
74 t_max = 0.2
75 max_steps = math.ceil(t_max / dt)
76 table_cap = math.ceil(t_max * table_rate) + 2
77 frame_cap = math.ceil(t_max * fps) + 2
78 qdd = np.zeros((table_cap, M))
79 q_frames = np.zeros((frame_cap, M))
80 force = np.zeros(max_steps + 1)
81 table_len = n_frames = 0
82 peak_force = contact_time = 0.0
83 z, zd = 0.0, -strike.velocity # stick tip, positive away from the plate; it arrives moving toward it
84 last_open, ever, F = 0.0, False, 0.0
85 step = 0
86 for step in range(max_steps):
87 t = step * dt
88 delta = float(g @ q) - z # penetration: the plate has come past the tip
89 F = strike.stiffness * delta**strike.exponent if delta > 0 else 0.0
90 if F > 0:
91 ever = True
92 last_open = t + dt
93 contact_time += dt
94 peak_force = max(peak_force, F)
95 force[step] = F
96 while table_len < table_cap and table_len / table_rate <= t + 1e-12:
97 qdd[table_len] = -2 * sigma * qd - omega**2 * q - F * g
98 table_len += 1
99 while n_frames < frame_cap and n_frames / fps <= t + 1e-12:
100 q_frames[n_frames] = q
101 n_frames += 1
102 if ever and t - last_open >= STICK_RELEASE:
103 break
104 qp = -F * g / omega**2
105 x = q - qp
106 q = p11 * x + p12 * qd + qp
107 qd = p21 * x + p22 * qd
108 zd += F / strike.mass * dt
109 z += zd * dt
110 else:
111 step = max_steps
112 t_end = step * dt
113 amplitude = q - 1j * (qd + sigma * q) / wd
114 return Contact(t_end, min(last_open, t_end), amplitude, qdd[:table_len], table_rate, q_frames[:n_frames], force[: step + 1], dt, peak_force, contact_time)
117def listener_weights(pl: Plate, listener: Listener) -> tuple[np.ndarray, np.ndarray, float]:
118 """Listener terms with the common delay removed, so that only the spread of the delays across the plate remains."""
119 weight, delay = pl.listener_terms(listener)
120 sel = weight > 0
121 delay = np.where(sel, delay - delay[sel].min(), 0.0)
122 return weight, delay, float(delay.max())
125def tail_start(c: Contact, max_delay: float) -> float:
126 """The time from which the tail formula holds at the listener: every node's
127 retarded time is past the last contact. The table path covers the output before it."""
128 return min(c.t_end, c.t_contact + max_delay + 2 / c.table_rate)
131def radiation_gains(modes: Modes, sigma: np.ndarray, weight: np.ndarray, delay: np.ndarray, batch: int = 256) -> np.ndarray:
132 """The complex radiation gain of each mode for the listener:
133 G_m = sum_k weight_k phi_m(k) exp(-s_m delay_k), with s_m = -sigma + i omega_d."""
134 sel = weight > 0
135 w, tau = weight[sel], delay[sel]
136 wd = np.sqrt(modes.omega**2 - sigma**2)
137 G = np.zeros(modes.count, dtype=complex)
138 for a in range(0, modes.count, batch):
139 b = min(modes.count, a + batch)
140 phi = modes.w[a:b][:, sel]
141 E = np.exp((sigma[a:b, None] - 1j * wd[a:b, None]) * tau[None, :])
142 G[a:b] = (phi * w[None, :] * E).sum(axis=1)
143 return G
146def attack_pressure(modes: Modes, c: Contact, weight: np.ndarray, delay: np.ndarray, fs: float, t_start: float, out: np.ndarray) -> None:
147 """The listener pressure over the attack from the nodal acceleration table:
148 p(t) = sum_k weight_k wdd_k(t - delay_k), linearly interpolated, added to
149 ``out`` for the samples before ``t_start``."""
150 n_out = min(len(out), math.ceil(t_start * fs - 1e-9))
151 if n_out <= 0:
152 return
153 length = min(len(c.qdd), math.ceil(t_start * c.table_rate) + 2)
154 sel = weight > 0
155 table = modes.w[:, sel].T @ c.qdd[:length].T # (R, length) nodal accelerations
156 x = (np.arange(n_out)[None, :] / fs - delay[sel][:, None]) * c.table_rate # (R, n_out)
157 i = np.floor(x).astype(int)
158 valid = (x >= 0) & (i + 1 < length)
159 i0 = np.clip(i, 0, length - 1)
160 i1 = np.clip(i + 1, 0, length - 1)
161 fr = x - i
162 rows = np.arange(sel.sum())[:, None]
163 val = (1 - fr) * table[rows, i0] + fr * table[rows, i1]
164 out[:n_out] += (weight[sel][:, None] * np.where(valid, val, 0.0)).sum(axis=0)
167def modal_tail(modes: Modes, sigma: np.ndarray, c: Contact, gains: np.ndarray, fs: float, t_start: float, out: np.ndarray, block: int = 4096) -> None:
168 """The listener pressure from ``t_start`` on: p(t) = sum_m Re[G_m s_m^2 A_m
169 exp(s_m (t - t_end))], added to ``out`` in blocks. Modes that have decayed
170 by 140 dB are dropped."""
171 n = len(out)
172 n0 = math.ceil(t_start * fs - 1e-9)
173 if n0 >= n:
174 return
175 wd = np.sqrt(modes.omega**2 - sigma**2)
176 s = -sigma + 1j * wd
177 coef = gains * s**2 * c.amplitude * np.exp(s * (n0 / fs - c.t_end)) # the value at sample n0
178 Z = np.exp(np.outer(np.arange(block) / fs, s)) # (block, M): z^b
179 scale = np.abs(coef)
180 for start in range(n0, n, block):
181 stop = min(n, start + block)
182 cur = coef * np.exp(s * ((start - n0) / fs))
183 active = np.abs(cur) > 1e-7 * scale
184 if not active.any():
185 break
186 out[start:stop] += (Z[: stop - start, active] @ cur[active]).real
189def movie_frames(modes: Modes, sigma: np.ndarray, c: Contact, fps: float, n_frames: int) -> np.ndarray:
190 """Movie frames of the transverse displacement, ``(n_frames, ny, nx)`` float32:
191 sum_m phi_m(k) q_m(t_f), with q from the contact tables before t_end and from the free decay after it."""
192 wd = np.sqrt(modes.omega**2 - sigma**2)
193 t = np.arange(n_frames) / fps - c.t_end
194 Q = np.exp(-sigma[None, :] * t[:, None]) * (c.amplitude.real[None, :] * np.cos(wd[None, :] * t[:, None]) - c.amplitude.imag[None, :] * np.sin(wd[None, :] * t[:, None]))
195 k = min(len(c.q_frames), n_frames)
196 Q[:k] = c.q_frames[:k]
197 g = modes.plate.grid
198 sel = g.mask.reshape(-1)
199 frames = np.zeros((n_frames, g.n), dtype=np.float32)
200 frames[:, sel] = Q @ modes.w.reshape(modes.count, -1)[:, sel]
201 return frames.reshape(n_frames, g.ny, g.nx)
204def synthesize(modes: Modes, strike: Strike = Strike(), listener: Listener = Listener(), *, duration: float = 4.0, fs: int = 48000, fps: float = 120.0, movie: bool = False) -> Sound:
205 """The sound of a strike on the plate by modal synthesis.
207 Returns a ``Sound`` with the listener pressure [Pa] at ``fs`` over
208 ``duration`` seconds after the impact and, with ``movie``, the transverse
209 displacement at ``fps`` frames per second. The strike and the listener can
210 be changed at no cost beyond this call; the modes are the expensive part.
211 """
212 pl = modes.plate
213 sigma = modes.sigma
214 g = strike_gains(modes, strike)
215 c = simulate_contact(modes, sigma, g, strike, fps, fs)
216 weight, delay, max_delay = listener_weights(pl, listener)
217 t_start = tail_start(c, max_delay)
218 n = round(duration * fs)
219 audio = np.zeros(n)
220 attack_pressure(modes, c, weight, delay, fs, t_start, audio)
221 modal_tail(modes, sigma, c, radiation_gains(modes, sigma, weight, delay), fs, t_start, audio)
222 frames = movie_frames(modes, sigma, c, fps, round(duration * fps)) if movie else None
223 info = {
224 "model": "linear",
225 "modes": modes.count,
226 "peak_force": c.peak_force,
227 "contact_time": c.contact_time,
228 "t_end": c.t_end,
229 "force": c.force,
230 "force_dt": c.dt,
231 }
232 return Sound(audio.astype(np.float32), fs, pl, strike, listener, frames, fps, info)