Add overlap_pursuit as a built-in: anchored matching pursuit over merged silence blocks (boundary onsets exact, coverage-forced interior events, duplicate pursuit with a same-type veto); 0.958 mean accuracy on the 3-type 20%-overlap export where the prior best scored 0.730
5 changed files+305−1
scripts/smoke-browser.tsmodified+1−1View file
@@ -57,7 +57,7 @@ try {
5757
5858 // accuracy of the two real baselines should beat 0.5 on demo data
5959 const accs = rows.map((r) => Number(r.split('|')[1].trim().split(' ')[0]));
60- if (rows.length !== 4) throw new Error(`expected 4 result rows, got ${rows.length}`);
60+ if (rows.length !== 5) throw new Error(`expected 5 result rows, got ${rows.length}`);
6161 if (!(accs[0] > 0.5 && accs[1] > 0.5)) throw new Error(`accuracies too low: ${accs.join(', ')}`);
6262 if (errors.length) throw new Error(`page errors:\n${errors.join('\n')}`);
6363 console.log('SMOKE PASS');
scripts/test-node.tsmodified+2−0View file
@@ -15,6 +15,7 @@ import { makeSyntheticLibrary } from './synthetic-lib.ts';
1515 import matchedFilterSrc from '../src/algorithms/template_matched_filter.m?raw';
1616 import onsetKnnSrc from '../src/algorithms/onset_knn.m?raw';
1717 import silenceGapKnnSrc from '../src/algorithms/silence_gap_knn.m?raw';
18+import overlapPursuitSrc from '../src/algorithms/overlap_pursuit.m?raw';
1819 import naiveSrc from '../src/algorithms/naive_threshold.m?raw';
1920
2021 const lib = makeSyntheticLibrary();
@@ -34,6 +35,7 @@ const algos = [
3435 { name: 'template_matched_filter', source: matchedFilterSrc },
3536 { name: 'onset_knn', source: onsetKnnSrc },
3637 { name: 'silence_gap_knn', source: silenceGapKnnSrc },
38+ { name: 'overlap_pursuit', source: overlapPursuitSrc },
3739 { name: 'naive_threshold', source: naiveSrc },
3840 ];
3941
scripts/test-offline.tsmodified+2−0View file
@@ -21,6 +21,7 @@ import { makeSyntheticLibrary } from './synthetic-lib.ts';
2121 import matchedFilterSrc from '../src/algorithms/template_matched_filter.m?raw';
2222 import onsetKnnSrc from '../src/algorithms/onset_knn.m?raw';
2323 import silenceGapKnnSrc from '../src/algorithms/silence_gap_knn.m?raw';
24+import overlapPursuitSrc from '../src/algorithms/overlap_pursuit.m?raw';
2425 import naiveSrc from '../src/algorithms/naive_threshold.m?raw';
2526
2627 const lib = makeSyntheticLibrary();
@@ -34,6 +35,7 @@ const algos = [
3435 { name: 'template_matched_filter', source: matchedFilterSrc },
3536 { name: 'onset_knn', source: onsetKnnSrc },
3637 { name: 'silence_gap_knn', source: silenceGapKnnSrc },
38+ { name: 'overlap_pursuit', source: overlapPursuitSrc },
3739 { name: 'naive_threshold', source: naiveSrc },
3840 ];
3941
src/algorithms/index.tsmodified+8−0View file
@@ -7,6 +7,7 @@ import type { AlgorithmDef } from '../bench/types.ts';
77 import matchedFilterSrc from './template_matched_filter.m?raw';
88 import onsetKnnSrc from './onset_knn.m?raw';
99 import silenceGapKnnSrc from './silence_gap_knn.m?raw';
10+import overlapPursuitSrc from './overlap_pursuit.m?raw';
1011 import naiveSrc from './naive_threshold.m?raw';
1112
1213 export const BUILTIN_ALGORITHMS: AlgorithmDef[] = [
@@ -31,6 +32,13 @@ export const BUILTIN_ALGORITHMS: AlgorithmDef[] = [
3132 enabled: true,
3233 builtin: true,
3334 },
35+ {
36+ id: 'builtin:overlap_pursuit',
37+ name: 'overlap_pursuit',
38+ source: overlapPursuitSrc,
39+ enabled: true,
40+ builtin: true,
41+ },
3442 {
3543 id: 'builtin:naive_threshold',
3644 name: 'naive_threshold',
src/algorithms/overlap_pursuit.madded+292−0View file
@@ -0,0 +1,292 @@
1+function [times, labels] = overlap_pursuit(x, fs, training)
2+% Silence segmentation plus per-block matching pursuit for overlaps.
3+%
4+% The recordings place each event, an exact snippet scaled by one gain,
5+% on a digitally silent background, so activity splits into blocks: a
6+% block exactly one snippet long holds a single event, and a longer
7+% block holds several overlapping events whose first onset is the block
8+% start and whose last onset is the block end minus one snippet length.
9+% Singles are classified by the spectral kNN of silence_gap_knn. For a
10+% multi-event block the two boundary onsets are taken as anchors and
11+% fitted by least squares against every training snippet (the best
12+% instance and gain per onset, refined by a few coordinate-descent
13+% sweeps that add each event back and refit it against the residual).
14+% Two pursuit stages then add interior events. First, coverage: every
15+% sample in a block must lie inside some event's snippet span, so while
16+% uncovered samples remain, the residual is scanned by FFT
17+% cross-correlation against the whole bank and the best-reducing onset
18+% whose span touches uncovered samples is added; this is a proof an
19+% interior event exists, not a guess. Second, duplicates: remaining
20+% residual energy is scanned the same way, and a candidate is accepted
21+% while its gain and energy reduction are plausible, unless an event of
22+% the same type already sits within a few hundred samples, since
23+% instance-to-instance mismatch of an already-fitted event is the usual
24+% source of such candidates and near-coincident same-type pairs are
25+% rare in comparison.
26+%
27+% Limitations: the block bookkeeping relies on the silent background,
28+% and nearly-coincident events of the same type are reported once. The
29+% snippet bank explains an event only up to instance variability
30+% (worst for the least stereotyped type), which is what bounds how far
31+% the residual can be trusted; with added background noise the
32+% thresholds below would need retuning.
33+nTypes = numel(training.snips);
34+L = size(training.snips{1}, 1);
35+pre = training.pre;
36+n = numel(x);
37+
38+% tuned constants
39+gmin = 0.35; % smallest plausible fitted gain for a pursuit event
40+redmin = 0.5; % smallest energy reduction for a duplicate event
41+prox = round(0.0156 * fs); % same-type exclusion radius, ~250 samples
42+maxExtra = 6; % coverage-pursuit iterations per block
43+maxDup = 4; % duplicate-pursuit iterations per block
44+
45+% ---- training bank ------------------------------------------------------
46+nTrain = 0;
47+for k = 1:nTypes
48+ nTrain = nTrain + size(training.snips{k}, 2);
49+end
50+B = zeros(L, nTrain);
51+bt = zeros(nTrain, 1);
52+c = 0;
53+for k = 1:nTypes
54+ S = training.snips{k};
55+ for j = 1:size(S, 2)
56+ c = c + 1;
57+ B(:, c) = S(:, j);
58+ bt(c) = k;
59+ end
60+end
61+n2 = zeros(nTrain, 1);
62+for j = 1:nTrain
63+ n2(j) = sum(B(:, j) .* B(:, j));
64+end
65+
66+% ---- kNN training features for single-event blocks ----------------------
67+F = zeros(nTrain, 8);
68+for j = 1:nTrain
69+ F(j, :) = spectral_features(B(pre+1:end, j), fs);
70+end
71+mu = mean(F, 1);
72+sd = std(F, 0, 1) + 1e-9;
73+Fz = (F - repmat(mu, nTrain, 1)) ./ repmat(sd, nTrain, 1);
74+
75+% ---- silence segmentation ------------------------------------------------
76+eps0 = 1e-4 * max(abs(x));
77+gapClose = round(0.002 * fs);
78+minDur = round(0.010 * fs);
79+bs = zeros(0, 1);
80+be = zeros(0, 1);
81+i = 1;
82+while i <= n
83+ if abs(x(i)) > eps0
84+ s = i;
85+ e = i;
86+ j = i + 1;
87+ quiet = 0;
88+ while j <= n && quiet <= gapClose
89+ if abs(x(j)) > eps0
90+ e = j;
91+ quiet = 0;
92+ else
93+ quiet = quiet + 1;
94+ end
95+ j = j + 1;
96+ end
97+ if e - s + 1 >= minDur
98+ bs(end+1, 1) = s;
99+ be(end+1, 1) = e;
100+ end
101+ i = e + gapClose + 1;
102+ else
103+ i = i + 1;
104+ end
105+end
106+
107+% ---- per-block detection and classification ------------------------------
108+times = zeros(0, 1);
109+labels = zeros(0, 1);
110+for bi = 1:numel(bs)
111+ s = bs(bi);
112+ D = be(bi) - s + 1;
113+ if D <= L + 2
114+ % single event: onset at the block start, spectral kNN label
115+ seg = x(s+pre : min(n, s+L-1));
116+ f = (spectral_features(seg, fs) - mu) ./ sd;
117+ d = sqrt(sum((Fz - repmat(f, nTrain, 1)).^2, 2));
118+ [ds, di] = sort(d);
119+ kNN = min(5, nTrain);
120+ votes = zeros(nTypes, 1);
121+ for v = 1:kNN
122+ votes(bt(di(v))) = votes(bt(di(v))) + 1 / (ds(v) + 1e-6);
123+ end
124+ [~, lab] = max(votes);
125+ times(end+1, 1) = (s - 1 + pre) / fs;
126+ labels(end+1, 1) = lab;
127+ else
128+ blk = x(s : s+D-1);
129+ [evOn, evType] = decompose(blk, B, bt, n2, L, gmin, redmin, prox, maxExtra, maxDup);
130+ for q = 1:numel(evOn)
131+ times(end+1, 1) = (s - 1 + evOn(q) - 1 + pre) / fs;
132+ labels(end+1, 1) = evType(q);
133+ end
134+ end
135+end
136+end
137+
138+function [evOn, evType] = decompose(blk, B, bt, n2, L, gmin, redmin, prox, maxExtra, maxDup)
139+% Anchored matching pursuit over one multi-event block.
140+D = numel(blk);
141+nl = D - L + 1;
142+nTrain = numel(bt);
143+
144+% conjugate template spectra for the FFT correlation scans
145+nfft = 2^nextpow2(D);
146+BF = zeros(nfft, nTrain) + 0i;
147+for j = 1:nTrain
148+ BF(:, j) = conj(fft(B(:, j), nfft));
149+end
150+
151+% events as parallel arrays, seeded with the two exact boundary anchors
152+evOn = [1; nl];
153+evType = zeros(2, 1);
154+evIdx = zeros(2, 1);
155+evGain = zeros(2, 1);
156+residual = blk;
157+
158+[residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);
159+
160+% coverage pursuit: every block sample must lie in some event span
161+for it = 1:maxExtra
162+ cov = zeros(D, 1);
163+ for q = 1:numel(evOn)
164+ cov(evOn(q) : evOn(q)+L-1) = 1;
165+ end
166+ nUnc = D - sum(cov);
167+ if nUnc == 0
168+ break;
169+ end
170+ % lags whose span [o, o+L-1] contains an uncovered sample
171+ m = conv(1 - cov, ones(L, 1));
172+ mask = zeros(nl, 1);
173+ for o = 1:nl
174+ mask(o) = m(o + L - 1) > 0;
175+ end
176+ [o, kk, idx, g, red] = scan_best(residual, BF, bt, n2, L, mask);
177+ if idx == 0
178+ break;
179+ end
180+ residual(o:o+L-1) = residual(o:o+L-1) - g * B(:, idx);
181+ evOn(end+1, 1) = o;
182+ evType(end+1, 1) = kk;
183+ evIdx(end+1, 1) = idx;
184+ evGain(end+1, 1) = g;
185+ [residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);
186+end
187+
188+% duplicate pursuit: events fully covered by other spans
189+mask = zeros(nl, 1) + 1;
190+for it = 1:maxDup
191+ [o, kk, idx, g, red] = scan_best(residual, BF, bt, n2, L, mask);
192+ if idx == 0 || red < redmin || g < gmin
193+ break;
194+ end
195+ clash = 0;
196+ for q = 1:numel(evOn)
197+ if evOn(q) == o || (abs(evOn(q) - o) < prox && evType(q) == kk)
198+ clash = 1;
199+ end
200+ end
201+ if clash
202+ break;
203+ end
204+ residual(o:o+L-1) = residual(o:o+L-1) - g * B(:, idx);
205+ evOn(end+1, 1) = o;
206+ evType(end+1, 1) = kk;
207+ evIdx(end+1, 1) = idx;
208+ evGain(end+1, 1) = g;
209+ [residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);
210+end
211+end
212+
213+function [residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, nsw)
214+% Coordinate descent: add each event back and refit its type/instance/gain.
215+[~, ord] = sort(evOn);
216+for it = 1:nsw
217+ for z = 1:numel(ord)
218+ q = ord(z);
219+ o = evOn(q);
220+ if evIdx(q) > 0
221+ residual(o:o+L-1) = residual(o:o+L-1) + evGain(q) * B(:, evIdx(q));
222+ end
223+ seg = residual(o:o+L-1);
224+ bestRed = 0;
225+ bestJ = 0;
226+ bestG = 0;
227+ for j = 1:numel(bt)
228+ cj = sum(seg .* B(:, j));
229+ if cj > 0
230+ red = cj * cj / n2(j);
231+ if red > bestRed
232+ bestRed = red;
233+ bestJ = j;
234+ bestG = cj / n2(j);
235+ end
236+ end
237+ end
238+ if bestJ > 0
239+ residual(o:o+L-1) = residual(o:o+L-1) - bestG * B(:, bestJ);
240+ evType(q) = bt(bestJ);
241+ evIdx(q) = bestJ;
242+ evGain(q) = bestG;
243+ end
244+ end
245+end
246+end
247+
248+function [o, kk, idx, g, red] = scan_best(residual, BF, bt, n2, L, mask)
249+% Best (onset, template) over masked lags by FFT cross-correlation;
250+% reduction is the residual energy removed by the least-squares gain.
251+nfft = size(BF, 1);
252+nl = numel(mask);
253+R = fft(residual, nfft);
254+o = 0; kk = 0; idx = 0; g = 0; red = 0;
255+for j = 1:numel(bt)
256+ cf = real(ifft(R .* BF(:, j)));
257+ cl = cf(1:nl);
258+ rv = (cl .* cl) / n2(j);
259+ rv = rv .* (cl > 0) .* mask;
260+ [m, mi] = max(rv);
261+ if m > red
262+ red = m;
263+ o = mi;
264+ kk = bt(j);
265+ idx = j;
266+ g = cl(mi) / n2(j);
267+ end
268+end
269+end
270+
271+function f = spectral_features(seg, fs)
272+% Seven log band-energy ratios (mean removed, level cancels) and the
273+% log spectral centroid; the extra low bands split the two bass types.
274+nfft = 2048;
275+m = min(numel(seg), nfft);
276+s2 = seg(1:m);
277+win = 0.5 - 0.5 * cos(2 * pi * (0:m-1)' / max(m - 1, 1));
278+S = fft(s2 .* win, nfft);
279+p = abs(S(1:nfft/2+1)).^2;
280+hz = fs / nfft;
281+edges = [0 100 200 400 800 1600 3200 6400];
282+b = zeros(1, 7);
283+for q = 1:7
284+ k0 = max(2, floor(edges(q) / hz) + 1);
285+ k1 = min(nfft/2 + 1, floor(edges(q+1) / hz));
286+ b(q) = log10(sum(p(k0:k1)) + 1e-12);
287+end
288+b = b - mean(b);
289+kv = (0:numel(p)-1)' * hz;
290+cen = sum(p .* kv) / (sum(p) + 1e-12);
291+f = [b, log10(cen + 1)];
292+end