7b92c9cAdd dataset generation scripts, ML fitting code and reference recordingsJeremy Magland 1"""Finds the simulated clips closest to a real recording.
3 python3 ml/match.py REAL.wav [--data DATASET_DIR] [--top 10]
5The recording is converted to mono at 48 kHz, its onset is found (the first
6sample of the impact transient), and the 3 s from the onset are compared with
7every simulated clip, which starts at the impact. Three distances are used,
8each insensitive to the recording level and to the broad coloring of the
9microphone and the room:
111. the fine structure of the mean log spectrum over the first second (the
12 log spectrum minus its smooth envelope, over 30 Hz to 5 kHz), which places
13 the partials;
142. the same over the whole clip;
153. the decay: the log energy in 24 bands (100 Hz to 5 kHz) at 12 times over
16 the clip, each band normalized to its own start, which captures how fast
17 each region of the spectrum dies away.
19The ranking uses a fourth, perceptual distance: the mean absolute difference
20in dB between log-mel-like spectrograms (48 log-spaced bands from 60 Hz to
216 kHz on 100 ms frames, relative to each clip's maximum, floored at -60 dB),
22which weighs timbre and decay as they are heard; the three distances above
23are reported alongside.
24"""
26import argparse
27import json
28import os
29import sys
30import wave
32import numpy as np
33from scipy.ndimage import median_filter
34from scipy.signal import resample_poly
36FS = 48000
37N_FFT = 8192
38HOP = 1024
39DUR = 3.0
42def read_wav_any(path):
43 with wave.open(path) as w:
44 ch, sw, fs, n = w.getnchannels(), w.getsampwidth(), w.getframerate(), w.getnframes()
45 raw = w.readframes(n)
46 if sw == 2:
47 x = np.frombuffer(raw, dtype="<i2").astype(np.float32) / 32768
48 elif sw == 3:
49 b = np.frombuffer(raw, dtype=np.uint8).reshape(-1, 3)
50 x = ((b[:, 0].astype(np.int32) | (b[:, 1].astype(np.int32) << 8) | (b[:, 2].astype(np.int32) << 16)) << 8 >> 8).astype(np.float32) / 2 ** 23
51 elif sw == 4:
52 x = np.frombuffer(raw, dtype="<i4").astype(np.float32) / 2 ** 31
53 else:
54 raise ValueError(f"unsupported sample width {sw}")
55 if ch > 1:
56 x = x.reshape(-1, ch).mean(1)
57 if fs != FS:
58 from math import gcd
59 g = gcd(fs, FS)
60 x = resample_poly(x, FS // g, fs // g).astype(np.float32)
61 return x
64def onset(x, frac=0.1, pre=0.002):
65 """First sample where the short-time envelope reaches `frac` of its peak, less `pre` seconds."""
66 env = np.abs(x)
67 k = int(0.001 * FS)
68 env = np.convolve(env, np.ones(k) / k, mode="same")
69 i = int(np.argmax(env >= frac * env.max()))
70 return max(0, i - int(pre * FS))
73def segment(x, start):
74 n = int(DUR * FS)
75 seg = np.zeros(n, np.float32)
76 part = x[start:start + n]
77 seg[: len(part)] = part
78 return seg / (np.abs(seg).max() + 1e-9)
81def stft_db(seg):
82 win = np.hanning(N_FFT + 1)[:-1].astype(np.float32)
83 n_frames = 1 + (len(seg) - N_FFT) // HOP
84 frames = np.lib.stride_tricks.as_strided(seg, shape=(n_frames, N_FFT), strides=(seg.strides[0] * HOP, seg.strides[0]))
85 return 20 * np.log10(np.abs(np.fft.rfft(frames * win, axis=1)) + 1e-9) # (frames, bins)
88def melspec(seg, n_bands=48, f_lo=60.0, f_hi=6000.0, frame=0.1, floor=-60.0):
89 """Log-mel-like spectrogram (log-spaced bands) in dB relative to the clip's
90 maximum, floored, on 100 ms frames: what the ear hears of timbre and decay."""
91 n = int(frame * FS)
92 m = len(seg) // n
93 frames = seg[: m * n].reshape(m, n) * np.hanning(n)
94 P = np.abs(np.fft.rfft(frames, axis=1)) ** 2
95 f = np.fft.rfftfreq(n, 1 / FS)
96 edges = np.geomspace(f_lo, f_hi, n_bands + 1)
97 B = np.stack([P[:, (f >= a) & (f < b)].sum(1) for a, b in zip(edges[:-1], edges[1:])], 1)
98 db = 10 * np.log10(B + 1e-18)
99 db -= db.max()
100 return np.maximum(db, floor)
103def features(seg):
104 S = stft_db(seg)
105 f = np.fft.rfftfreq(N_FFT, 1 / FS)
106 band = (f >= 30) & (f <= 5000)
107 fps = FS / HOP
108 def fine(mean_db):
109 env = median_filter(mean_db, size=int(200 / (FS / N_FFT)) | 1, mode="nearest")
110 v = (mean_db - env)[band]
111 return v / (np.linalg.norm(v) + 1e-9)
112 early = fine(S[: int(1.0 * fps)].mean(0))
113 whole = fine(S.mean(0))
114 # Decay profile: 24 log-spaced bands, 12 time steps, each band relative to its first step.
115 edges = np.geomspace(100, 5000, 25)
116 P = 10 ** (S / 10)
117 bands = np.stack([P[:, (f >= a) & (f < b)].sum(1) for a, b in zip(edges[:-1], edges[1:])], 1) # (frames, 24)
118 steps = np.array_split(np.arange(S.shape[0]), 12)
119 prof = np.stack([bands[s].mean(0) for s in steps]) # (12, 24)
120 prof = 10 * np.log10(prof + 1e-12)
121 prof = prof - prof[0]
122 return early, whole, prof.ravel(), melspec(seg)
125def main():
126 ap = argparse.ArgumentParser()
127 ap.add_argument("real")
128 ap.add_argument("--data", default="dataset")
129 ap.add_argument("--top", type=int, default=10)
130 ap.add_argument("--onset-frac", type=float, default=0.1)
131 args = ap.parse_args()
132 x = read_wav_any(args.real)
133 st = onset(x, args.onset_frac)
134 print(f"{args.real}: {len(x) / FS:.2f} s, onset at {st / FS:.3f} s (peak {np.abs(x).max():.3f})")
135 ref = features(segment(x, st))
136 rows = {}
137 with open(os.path.join(args.data, "index.jsonl")) as fh:
138 for line in fh:
139 if line.strip():
140 r = json.loads(line)
141 rows[r["id"]] = r
142 ids, D = [], []
143 for id_ in sorted(rows):
144 p = os.path.join(args.data, "clips", id_ + ".wav")
145 if not os.path.exists(p):
146 continue
147 seg = segment(read_wav_any(p), 0)
148 fe = features(seg)
149 d1 = 1 - float(ref[0] @ fe[0])
150 d2 = 1 - float(ref[1] @ fe[1])
151 d3 = float(np.sqrt(np.mean((ref[2] - fe[2]) ** 2)))
152 d4 = float(np.mean(np.abs(ref[3] - fe[3])))
153 ids.append(id_); D.append((d1, d2, d3, d4))
154 D = np.array(D)
155 # Main score: the mean absolute dB difference of the log-mel spectrograms,
156 # which weighs timbre and decay the way they are heard; the other three
157 # distances are reported for reference.
158 score = D[:, 3]
159 order = np.argsort(score, kind="stable")
160 print(f"{len(ids)} simulated clips compared\n")
161 print(f"{'id':>6} {'mel dB':>6} {'d_early':>7} {'d_whole':>7} {'d_decay':>7} | {'E GPa':>5} {'rho':>5} {'h mm':>5} {'W x H in':>11} {'v':>4} {'strike':>10} {'sig0':>5} {'sig1e4':>6} {'p2p mm':>6} features")
162 for i in order[: args.top]:
163 r = rows[ids[i]]
164 print(f"{ids[i]:>6} {score[i]:6.2f} {D[i, 0]:7.3f} {D[i, 1]:7.3f} {D[i, 2]:7.2f} | {r['E'] / 1e9:5.0f} {r['rho']:5.0f} {r['h'] * 1e3:5.2f} {r['Lx'] / 0.0254:5.1f} x {r['Ly'] / 0.0254:4.1f} {r['strikeVelocity']:4.1f} {r['strikeFracX']:4.2f},{r['strikeFracY']:4.2f} {r['sigma0']:5.2f} {r['sigma1'] * 1e4:6.1f} {r['peakToPeak'] * 1e3:6.1f} {r['features'][:70]}")
165 print(f"\nmedian over all clips: mel {np.median(D[:, 3]):.2f} dB, d_early {np.median(D[:, 0]):.3f}, d_whole {np.median(D[:, 1]):.3f}, d_decay {np.median(D[:, 2]):.2f}")
166 return [ids[i] for i in order[: args.top]]
169if __name__ == "__main__":
170 main()