"""Finds the simulated clips closest to a real recording. python3 ml/match.py REAL.wav [--data DATASET_DIR] [--top 10] The recording is converted to mono at 48 kHz, its onset is found (the first sample of the impact transient), and the 3 s from the onset are compared with every simulated clip, which starts at the impact. Three distances are used, each insensitive to the recording level and to the broad coloring of the microphone and the room: 1. the fine structure of the mean log spectrum over the first second (the log spectrum minus its smooth envelope, over 30 Hz to 5 kHz), which places the partials; 2. the same over the whole clip; 3. the decay: the log energy in 24 bands (100 Hz to 5 kHz) at 12 times over the clip, each band normalized to its own start, which captures how fast each region of the spectrum dies away. The ranking uses a fourth, perceptual distance: the mean absolute difference in dB between log-mel-like spectrograms (48 log-spaced bands from 60 Hz to 6 kHz on 100 ms frames, relative to each clip's maximum, floored at -60 dB), which weighs timbre and decay as they are heard; the three distances above are reported alongside. """ import argparse import json import os import sys import wave import numpy as np from scipy.ndimage import median_filter from scipy.signal import resample_poly FS = 48000 N_FFT = 8192 HOP = 1024 DUR = 3.0 def read_wav_any(path): with wave.open(path) as w: ch, sw, fs, n = w.getnchannels(), w.getsampwidth(), w.getframerate(), w.getnframes() raw = w.readframes(n) if sw == 2: x = np.frombuffer(raw, dtype="> 8).astype(np.float32) / 2 ** 23 elif sw == 4: x = np.frombuffer(raw, dtype=" 1: x = x.reshape(-1, ch).mean(1) if fs != FS: from math import gcd g = gcd(fs, FS) x = resample_poly(x, FS // g, fs // g).astype(np.float32) return x def onset(x, frac=0.1, pre=0.002): """First sample where the short-time envelope reaches `frac` of its peak, less `pre` seconds.""" env = np.abs(x) k = int(0.001 * FS) env = np.convolve(env, np.ones(k) / k, mode="same") i = int(np.argmax(env >= frac * env.max())) return max(0, i - int(pre * FS)) def segment(x, start): n = int(DUR * FS) seg = np.zeros(n, np.float32) part = x[start:start + n] seg[: len(part)] = part return seg / (np.abs(seg).max() + 1e-9) def stft_db(seg): win = np.hanning(N_FFT + 1)[:-1].astype(np.float32) n_frames = 1 + (len(seg) - N_FFT) // HOP frames = np.lib.stride_tricks.as_strided(seg, shape=(n_frames, N_FFT), strides=(seg.strides[0] * HOP, seg.strides[0])) return 20 * np.log10(np.abs(np.fft.rfft(frames * win, axis=1)) + 1e-9) # (frames, bins) def melspec(seg, n_bands=48, f_lo=60.0, f_hi=6000.0, frame=0.1, floor=-60.0): """Log-mel-like spectrogram (log-spaced bands) in dB relative to the clip's maximum, floored, on 100 ms frames: what the ear hears of timbre and decay.""" n = int(frame * FS) m = len(seg) // n frames = seg[: m * n].reshape(m, n) * np.hanning(n) P = np.abs(np.fft.rfft(frames, axis=1)) ** 2 f = np.fft.rfftfreq(n, 1 / FS) edges = np.geomspace(f_lo, f_hi, n_bands + 1) B = np.stack([P[:, (f >= a) & (f < b)].sum(1) for a, b in zip(edges[:-1], edges[1:])], 1) db = 10 * np.log10(B + 1e-18) db -= db.max() return np.maximum(db, floor) def features(seg): S = stft_db(seg) f = np.fft.rfftfreq(N_FFT, 1 / FS) band = (f >= 30) & (f <= 5000) fps = FS / HOP def fine(mean_db): env = median_filter(mean_db, size=int(200 / (FS / N_FFT)) | 1, mode="nearest") v = (mean_db - env)[band] return v / (np.linalg.norm(v) + 1e-9) early = fine(S[: int(1.0 * fps)].mean(0)) whole = fine(S.mean(0)) # Decay profile: 24 log-spaced bands, 12 time steps, each band relative to its first step. edges = np.geomspace(100, 5000, 25) P = 10 ** (S / 10) bands = np.stack([P[:, (f >= a) & (f < b)].sum(1) for a, b in zip(edges[:-1], edges[1:])], 1) # (frames, 24) steps = np.array_split(np.arange(S.shape[0]), 12) prof = np.stack([bands[s].mean(0) for s in steps]) # (12, 24) prof = 10 * np.log10(prof + 1e-12) prof = prof - prof[0] return early, whole, prof.ravel(), melspec(seg) def main(): ap = argparse.ArgumentParser() ap.add_argument("real") ap.add_argument("--data", default="dataset") ap.add_argument("--top", type=int, default=10) ap.add_argument("--onset-frac", type=float, default=0.1) args = ap.parse_args() x = read_wav_any(args.real) st = onset(x, args.onset_frac) print(f"{args.real}: {len(x) / FS:.2f} s, onset at {st / FS:.3f} s (peak {np.abs(x).max():.3f})") ref = features(segment(x, st)) rows = {} with open(os.path.join(args.data, "index.jsonl")) as fh: for line in fh: if line.strip(): r = json.loads(line) rows[r["id"]] = r ids, D = [], [] for id_ in sorted(rows): p = os.path.join(args.data, "clips", id_ + ".wav") if not os.path.exists(p): continue seg = segment(read_wav_any(p), 0) fe = features(seg) d1 = 1 - float(ref[0] @ fe[0]) d2 = 1 - float(ref[1] @ fe[1]) d3 = float(np.sqrt(np.mean((ref[2] - fe[2]) ** 2))) d4 = float(np.mean(np.abs(ref[3] - fe[3]))) ids.append(id_); D.append((d1, d2, d3, d4)) D = np.array(D) # Main score: the mean absolute dB difference of the log-mel spectrograms, # which weighs timbre and decay the way they are heard; the other three # distances are reported for reference. score = D[:, 3] order = np.argsort(score, kind="stable") print(f"{len(ids)} simulated clips compared\n") 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") for i in order[: args.top]: r = rows[ids[i]] 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]}") 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}") return [ids[i] for i in order[: args.top]] if __name__ == "__main__": main()