1"""Spectral diagnostics of the radiated pressure: per-octave band envelopes, onset times and the
2spectral centroid, which together show whether energy climbs from low to high frequency (the
3nonlinear upward cascade that gives a crash its build-up).
4"""
5from __future__ import annotations
7from dataclasses import dataclass
9import numpy as np
10from scipy.signal import stft
12DEFAULT_BANDS = (60, 125, 250, 500, 1000, 2000, 4000, 8000, 16000)
16class BandEnvelope:
17 lo: float
18 hi: float
19 onset: float # first time the band reaches -20 dB of its own peak [s]
20 peak_time: float # time of the band's peak [s]
21 level_db: float # band peak relative to the loudest band [dB]
22 envelope: np.ndarray # band energy vs time (linear)
26class CascadeAnalysis:
27 times: np.ndarray
28 freqs: np.ndarray
29 power: np.ndarray # STFT power (n_freq, n_time)
30 bands: list[BandEnvelope]
31 centroid: np.ndarray # spectral centroid vs time [Hz]
33 def centroid_at(self, t: float) -> float:
34 return float(self.centroid[min(np.searchsorted(self.times, t), len(self.times) - 1)])
37def analyze_cascade(p: np.ndarray, fs: int, bands=DEFAULT_BANDS, nperseg: int = 256,
38 noverlap: int = 224) -> CascadeAnalysis:
39 """Band envelopes at roughly 6 ms resolution (for fs = 44.1 kHz and the default hop)."""
40 f_, t_, Z = stft(p, fs, nperseg=nperseg, noverlap=noverlap)
41 Pw = np.abs(Z) ** 2
42 tot_max = Pw.sum(0).max() or 1.0
43 out = []
44 for lo, hi in zip(bands[:-1], bands[1:]):
45 m = (f_ >= lo) & (f_ < hi)
46 env = Pw[m].sum(0)
47 if env.max() <= 0:
48 continue
49 env_db = 10 * np.log10(env / env.max() + 1e-30)
50 out.append(BandEnvelope(lo, hi, float(t_[np.argmax(env_db > -20)]), float(t_[np.argmax(env)]),
51 float(10 * np.log10(env.max() / tot_max)), env))
52 cen = (Pw * f_[:, None]).sum(0) / (Pw.sum(0) + 1e-30)
53 return CascadeAnalysis(t_, f_, Pw, out, cen)
56def format_cascade_report(a: CascadeAnalysis, label: str = "") -> str:
57 lines = [f"[{label}] band onset(-20 dB) peak level vs loudest band"]
58 for b in a.bands:
59 lines.append(f" {b.lo:5.0f}-{b.hi:5.0f} Hz {b.onset * 1e3:7.1f} ms {b.peak_time * 1e3:7.1f} ms"
60 f" {b.level_db:6.1f} dB")
61 for tq in (0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5):
62 if tq < a.times[-1]:
63 lines.append(f" spectral centroid at {tq * 1e3:5.0f} ms: {a.centroid_at(tq):7.0f} Hz")
64 return "\n".join(lines)
67def plot_cascade(a: CascadeAnalysis, path: str, title: str = "", fs: float = 44100.0, t_max: float = 0.5) -> str:
68 """Spectrogram and band envelopes (log-log). Requires matplotlib."""
69 import matplotlib
70 matplotlib.use("Agg")
71 import matplotlib.pyplot as plt
73 tot_max = a.power.sum(0).max() or 1.0
74 sel = a.times <= t_max
75 fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
76 ax[0].pcolormesh(a.times[sel] * 1e3, a.freqs, 10 * np.log10(a.power[:, sel] / a.power.max() + 1e-12),
77 vmin=-90, vmax=0, shading="auto")
78 ax[0].set_yscale("log")
79 ax[0].set_ylim(50, fs / 2)
80 ax[0].set_xlabel("t [ms]")
81 ax[0].set_ylabel("f [Hz]")
82 ax[0].set_title(title)
83 for b in a.bands:
84 ax[1].plot(a.times[sel] * 1e3, 10 * np.log10(b.envelope[sel] / tot_max + 1e-30), label=f"{b.lo:.0f}-{b.hi:.0f} Hz")
85 ax[1].set_ylim(-90, 5)
86 ax[1].set_xscale("log")
87 ax[1].set_xlim(2, t_max * 1e3)
88 ax[1].set_xlabel("t [ms]")
89 ax[1].set_ylabel("band energy [dB re loudest]")
90 ax[1].legend(fontsize=7, ncol=2)
91 ax[1].grid(alpha=0.3)
92 ax[1].set_title("band envelopes (a cascade shows staggered rises, high bands last)")
93 plt.tight_layout()
94 plt.savefig(path, dpi=120)
95 plt.close(fig)
96 return path