3688ceacymsim: Python port of cymbalsim with a modal linear model and a PyTorch von Karman modelJeremy Magland 1"""The modal synthesis: a single mode follows the closed-form damped sinusoid, the
2radiation gain of a piston equals the analytic Rayleigh sum, the attack and the
3tail join, a lighter hit is proportionally quieter."""
5import numpy as np
7from cymsim import Listener, Plate, Strike, compute_modes, make_grid, synthesize
8from cymsim.constants import LN1000, RHO_AIR
9from cymsim.linear.synth import listener_weights, movie_frames, radiation_gains, simulate_contact, strike_gains, tail_start
11from conftest import STEEL
14def small_model():
15 pl = Plate.uniform(make_grid(0.2, 0.15, 0.0125), STEEL, 1e-3)
16 return pl, compute_modes(pl, 4000, damping=False)
19def test_single_mode_follows_the_closed_form():
20 pl, ms = small_model()
21 one = ms.select([0])
22 one.sigma_material[:] = LN1000 # a decay time of exactly 1 s
23 sigma = one.sigma
24 assert abs(sigma[0] - LN1000) < 1e-9
25 strike = Strike()
26 g = strike_gains(one, strike)
27 fps = 100
28 c = simulate_contact(one, sigma, g, strike, fps, 48000)
29 assert c.peak_force > 0 and c.contact_time > 0
30 frames = movie_frames(one, sigma, c, fps, 10)
31 w = 2 * np.pi * one.freqs[0]
32 wd = np.sqrt(w * w - sigma[0] ** 2)
33 j, i = pl.grid.ny // 2, pl.grid.nx // 3
34 for f in range(1, 10):
35 t = f / fps - c.t_end
36 if t < 0:
37 continue
38 q = np.exp(-sigma[0] * t) * (c.amplitude[0].real * np.cos(wd * t) - c.amplitude[0].imag * np.sin(wd * t))
39 assert abs(frames[f, j, i] - one.w[0, j, i] * q) < 1e-9 + 1e-6 * abs(one.w[0, j, i] * q)
42def test_piston_radiation_gain_equals_the_analytic_sum():
43 pl, ms = small_model()
44 one = ms.select([0])
45 phi = 1 / np.sqrt(pl.mass)
46 one.w[0] = np.where(pl.grid.mask, phi, 0.0)
47 one.freqs[0] = 1000.0
48 listener = Listener(distance=1.0, angle=0)
49 weight, delay, _ = listener_weights(pl, listener)
50 G = radiation_gains(one, np.zeros(1), weight, delay)[0]
51 w = 2 * np.pi * 1000
52 sel = weight > 0
53 expected = (weight[sel] * phi * np.exp(-1j * w * delay[sel])).sum()
54 assert abs(G - expected) < 1e-9 * abs(expected)
55 assert abs(abs(G) / (RHO_AIR * pl.area / (2 * np.pi) * phi) - 1) < 0.05
58def test_attack_and_tail_join_and_the_tail_decays():
59 pl, ms = small_model()
60 ms.sigma_material[:] = 3e-4 * ms.omega / 2
61 strike, listener = Strike(velocity=2), Listener()
62 fs = 48000
63 s = synthesize(ms, strike, listener, duration=0.2, fs=fs, movie=True)
64 assert len(s.pressure) == round(0.2 * fs) and s.frames.shape == (24, pl.grid.ny, pl.grid.nx)
65 g = strike_gains(ms, strike)
66 c = simulate_contact(ms, ms.sigma, g, strike, 120, fs)
67 _, _, max_delay = listener_weights(pl, listener)
68 n0 = int(np.ceil(tail_start(c, max_delay) * fs - 1e-9))
69 steps = np.abs(np.diff(s.pressure[n0 + 1 : n0 + 200]))
70 jump = abs(s.pressure[n0] - s.pressure[n0 - 1])
71 assert jump < 3 * steps.max()
72 rms = lambda a: np.sqrt((a.astype(float) ** 2).mean()) # noqa: E731
73 assert rms(s.pressure[-960:]) < rms(s.pressure[:960])
74 assert s.info["peak_force"] == c.peak_force and s.peak_spl > 60
77def test_a_lighter_hit_is_proportionally_quieter():
78 _, ms = small_model()
79 a = synthesize(ms, Strike(velocity=0.2), duration=0.05).pressure
80 b = synthesize(ms, Strike(velocity=0.1), duration=0.05).pressure
81 ratio = np.abs(a).max() / np.abs(b).max()
82 # The Hertz contact is nonlinear (a faster hit is shorter and brighter), so the ratio exceeds 2 somewhat.
83 assert 1.8 < ratio < 3.5
86def test_sound_outputs(tmp_path):
87 _, ms = small_model()
88 s = synthesize(ms, duration=0.1)
89 s.to_wav(tmp_path / "a.wav")
90 import wave
92 with wave.open(str(tmp_path / "a.wav")) as w:
93 assert w.getnchannels() == 1 and w.getframerate() == 48000 and w.getnframes() == 4800
94 assert np.abs(s.samples("peak")).max() <= 1 and abs(s.peak_dbfs("peak") + 1) < 0.01
95 sp = s.spectrogram(512, 128)
96 assert sp.db.shape == (257, len(sp.times)) and sp.db.max() == 0
97 assert "linear" in repr(s)