1"""Spectrogram comparison of the recording, the simulation and the matched outputs.
3Usage: python3 plot.py OUT_PREFIX (reads the wavs written by match.py)
4"""
5import sys, wave, numpy as np
6import matplotlib
7matplotlib.use("Agg")
8import matplotlib.pyplot as plt
9from scipy.signal import stft
11FS = 48000
14def load(path):
15 w = wave.open(path)
16 return np.frombuffer(w.readframes(w.getnframes()), dtype="<i2").astype(float) / 32768
19out = sys.argv[1]
20files = [("recording", f"{out}_reference_aligned.wav"), ("simulation", f"{out}_sim_aligned.wav"),
21 ("matched: decay rates", f"{out}_decay.wav"), ("matched: full envelope", f"{out}_envelope.wav")]
22fig, axes = plt.subplots(2, 2, figsize=(13, 8), sharex=True, sharey=True)
23for ax, (title, path) in zip(axes.flat, files):
24 x = load(path)
25 f, t, Z = stft(x, FS, nperseg=4096, noverlap=4096 - 512)
26 S = 20 * np.log10(np.abs(Z) + 1e-9)
27 S -= S.max()
28 ax.pcolormesh(t, f, S, vmin=-90, vmax=0, cmap="magma", shading="auto")
29 ax.set_yscale("log")
30 ax.set_ylim(25, 12000)
31 ax.set_title(title)
32 ax.set_xlabel("time [s]")
33 ax.set_ylabel("frequency [Hz]")
34fig.suptitle("Reference hit versus simulated strike (dB relative to each clip's peak, 90 dB range)")
35fig.tight_layout()
36fig.savefig(f"{out}_spectrograms.png", dpi=110)
37print(f"wrote {out}_spectrograms.png")
39# Band decay curves for a few bands.
40a = np.load(f"{out}_analysis.npz")
41t, Er, Es, centers, floor = a["t"], a["Er"], a["Es"], a["centers"], a["floor"]
42fig, axes = plt.subplots(2, 4, figsize=(14, 6), sharex=True)
43for ax, fc in zip(axes.flat, [35, 89, 178, 356, 713, 1425, 2851, 5702]):
44 b = int(np.argmin(np.abs(centers - fc)))
45 ax.plot(t, Er[b] - Er[b, np.searchsorted(t, 0.1)], label="recording")
46 ax.plot(t, Es[b] - Es[b, np.searchsorted(t, 0.1)], label="simulation")
47 ax.axhline(floor[b] - Er[b, np.searchsorted(t, 0.1)], color="gray", ls=":", lw=1, label="recording noise floor")
48 ax.set_title(f"{centers[b]:.0f} Hz band")
49 ax.set_ylim(-80, 10)
50 ax.grid(alpha=0.3)
51axes[0, 0].legend(fontsize=8)
52for ax in axes[1]:
53 ax.set_xlabel("time after onset [s]")
54for ax in axes[:, 0]:
55 ax.set_ylabel("level rel. 0.1 s [dB]")
56fig.tight_layout()
57fig.savefig(f"{out}_band_decay.png", dpi=110)
58print(f"wrote {out}_band_decay.png")