7b92c9cAdd dataset generation scripts, ML fitting code and reference recordingsJeremy Magland 1"""A similarity score for struck-sheet sounds that compares the partials with a
2tolerance in frequency, and their levels and decays, without requiring the
3lines to coincide.
5Representation. Each clip is cut into K time slices (0-0.3, 0.3-1, 1-2, 2-3 s
6by default). For each slice the mean power spectrum (8192-point windows,
75.9 Hz bins) is mapped onto a log-frequency axis (96 bins per octave, 30 Hz to
85 kHz) and smoothed there with a Gaussian whose width is the frequency
9tolerance (2 percent by default), then converted to dB relative to the clip's
10loudest value and floored. A line thus becomes a bump of width `tol`: two
11lines a fraction of `tol` apart overlap almost completely, lines several `tol`
12apart do not overlap at all, and in between the overlap falls off smoothly.
13The level of a line is the height of its bump, and its decay is how the bump
14shrinks from slice to slice, so the level, the position, and the decay of
15every partial are all in the K x n_bins picture.
17Score. The mean absolute difference in dB between the two pictures, over the
18cells where either clip is above the floor (so silent regions do not dilute
19it), with the later slices weighted equally to the first. Zero for identical
20clips; smooth in the parameters of a simulation; insensitive to the level and
21(through the log axis) treating a 2 percent shift at 100 Hz like one at 1 kHz.
23 python3 ml/linesim.py TARGET.wav CLIP.wav [CLIP.wav ...] [--tol 0.02]
24 python3 ml/linesim.py TARGET.wav --data DATASET_DIR --top 10 # rank a dataset
25"""
27import argparse
28import json
29import os
30import sys
32import numpy as np
34sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
35from match import read_wav_any, onset, segment, FS # noqa: E402
37N_FFT = 8192
38HOP = 1024
39F_LO, F_HI = 30.0, 5000.0
40BINS_PER_OCTAVE = 96
41SLICES = ((0, 0.3), (0.3, 1.0), (1.0, 2.0), (2.0, 3.0))
42FLOOR = -70.0
44_fb_cache = {}
47def logf_axis():
48 n = int(np.log2(F_HI / F_LO) * BINS_PER_OCTAVE)
49 return F_LO * 2 ** (np.arange(n) / BINS_PER_OCTAVE)
52def smoothing_matrix(tol):
53 """(n_log, n_lin): maps a linear power spectrum onto the log axis and smooths it with a Gaussian of relative width tol."""
54 if tol in _fb_cache:
55 return _fb_cache[tol]
56 f_lin = np.arange(N_FFT // 2 + 1) * FS / N_FFT
57 centers = logf_axis()
58 W = np.exp(-0.5 * ((np.log(f_lin[None, 1:] / centers[:, None])) / tol) ** 2)
59 W[:, :] /= np.maximum(W.sum(1, keepdims=True), 1e-12)
60 W = np.concatenate([np.zeros((len(centers), 1)), W], 1).astype(np.float32)
61 _fb_cache[tol] = W
62 return W
65def picture(seg, tol=0.02):
66 """(K, n_log) dB picture of the clip; see the module docstring."""
67 win = np.hanning(N_FFT + 1)[:-1].astype(np.float32)
68 n_frames = 1 + (len(seg) - N_FFT) // HOP
69 frames = np.lib.stride_tricks.as_strided(seg, shape=(n_frames, N_FFT), strides=(seg.strides[0] * HOP, seg.strides[0]))
70 P = np.abs(np.fft.rfft(frames * win, axis=1)) ** 2
71 fps = FS / HOP
72 W = smoothing_matrix(tol)
73 out = []
74 for a, b in SLICES:
75 sl = P[int(a * fps): max(int(a * fps) + 1, int(b * fps))].mean(0)
76 out.append(W @ sl)
77 db = 10 * np.log10(np.array(out) + 1e-20)
78 db -= db.max()
79 return np.maximum(db, FLOOR)
82def distance(A, B):
83 """Mean |dB difference| over the cells where either picture is above the floor, averaged over slices."""
84 per_slice = []
85 for a, b in zip(A, B):
86 m = (a > FLOOR + 1) | (b > FLOOR + 1)
87 per_slice.append(float(np.abs(a - b)[m].mean()) if m.any() else 0.0)
88 return float(np.mean(per_slice)), per_slice
91def load(path, align):
92 x = read_wav_any(path)
93 return segment(x, onset(x) if align else 0)
96def main():
97 ap = argparse.ArgumentParser()
98 ap.add_argument("target")
99 ap.add_argument("clips", nargs="*")
100 ap.add_argument("--data", default=None)
101 ap.add_argument("--top", type=int, default=10)
102 ap.add_argument("--tol", type=float, default=0.02)
103 ap.add_argument("--no-align", action="store_true", help="do not detect the onset of the target (simulated clips start at the impact)")
104 args = ap.parse_args()
105 T = picture(load(args.target, not args.no_align), args.tol)
106 if args.clips:
107 for c in args.clips:
108 d, parts = distance(T, picture(load(c, False), args.tol))
109 print(f"{os.path.basename(c):>16}: {d:6.2f} dB slices " + " ".join(f"{p:.2f}" for p in parts))
110 if args.data:
111 rows = {json.loads(l)["id"]: json.loads(l) for l in open(os.path.join(args.data, "index.jsonl")) if l.strip()}
112 ids = sorted(i for i in rows if os.path.exists(os.path.join(args.data, "clips", i + ".wav")))
113 cache = os.path.join(args.data, "ml", f"linesim-tol{args.tol:g}.npz")
114 if os.path.exists(cache) and list(np.load(cache)["ids"]) == ids:
115 pics = np.load(cache)["pics"]
116 else:
117 pics = np.stack([picture(load(os.path.join(args.data, "clips", i + ".wav"), False), args.tol) for i in ids]).astype(np.float32)
118 np.savez(cache, ids=np.array(ids), pics=pics)
119 d = np.array([distance(T, p)[0] for p in pics])
120 order = np.argsort(d)
121 print(f"{len(ids)} clips; median distance {np.median(d):.2f} dB")
122 for k in order[: args.top]:
123 r = rows[ids[k]]
124 print(f"{ids[k]} {d[k]:6.2f} dB | E {r['E'] / 1e9:.0f} rho {r['rho']:.0f} h {r['h'] * 1e3:.2f} {r['Lx'] / 0.0254:.1f}x{r['Ly'] / 0.0254:.1f} in v {r['strikeVelocity']:.1f} sig0 {r['sigma0']:.2f} sig1 {r['sigma1'] * 1e4:.1f} p2p {r['peakToPeak'] * 1e3:.1f}")
125 return [ids[k] for k in order[: args.top]]
128if __name__ == "__main__":
129 main()