7facb2fAdd 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.730Jeremy Magland 1function [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.
33nTypes = numel(training.snips);
34L = size(training.snips{1}, 1);
35pre = training.pre;
36n = numel(x);
38% tuned constants
39gmin = 0.35; % smallest plausible fitted gain for a pursuit event
40redmin = 0.5; % smallest energy reduction for a duplicate event
41prox = round(0.0156 * fs); % same-type exclusion radius, ~250 samples
42maxExtra = 6; % coverage-pursuit iterations per block
43maxDup = 4; % duplicate-pursuit iterations per block
45% ---- training bank ------------------------------------------------------
46nTrain = 0;
47for k = 1:nTypes
48 nTrain = nTrain + size(training.snips{k}, 2);
49end
50B = zeros(L, nTrain);
51bt = zeros(nTrain, 1);
52c = 0;
53for 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
60end
61n2 = zeros(nTrain, 1);
62for j = 1:nTrain
63 n2(j) = sum(B(:, j) .* B(:, j));
64end
66% ---- kNN training features for single-event blocks ----------------------
67F = zeros(nTrain, 8);
68for j = 1:nTrain
69 F(j, :) = spectral_features(B(pre+1:end, j), fs);
70end
71mu = mean(F, 1);
72sd = std(F, 0, 1) + 1e-9;
73Fz = (F - repmat(mu, nTrain, 1)) ./ repmat(sd, nTrain, 1);
75% ---- silence segmentation ------------------------------------------------
76eps0 = 1e-4 * max(abs(x));
77gapClose = round(0.002 * fs);
78minDur = round(0.010 * fs);
79bs = zeros(0, 1);
80be = zeros(0, 1);
81i = 1;
82while 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
105end
107% ---- per-block detection and classification ------------------------------
108times = zeros(0, 1);
109labels = zeros(0, 1);
110for 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
135end
136end
138function [evOn, evType] = decompose(blk, B, bt, n2, L, gmin, redmin, prox, maxExtra, maxDup)
139% Anchored matching pursuit over one multi-event block.
140D = numel(blk);
141nl = D - L + 1;
142nTrain = numel(bt);
144% conjugate template spectra for the FFT correlation scans
145nfft = 2^nextpow2(D);
146BF = zeros(nfft, nTrain) + 0i;
147for j = 1:nTrain
148 BF(:, j) = conj(fft(B(:, j), nfft));
149end
151% events as parallel arrays, seeded with the two exact boundary anchors
152evOn = [1; nl];
153evType = zeros(2, 1);
154evIdx = zeros(2, 1);
155evGain = zeros(2, 1);
156residual = blk;
158[residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);
160% coverage pursuit: every block sample must lie in some event span
161for 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);
186end
188% duplicate pursuit: events fully covered by other spans
189mask = zeros(nl, 1) + 1;
190for 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);
210end
211end
213function [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);
216for 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
245end
246end
248function [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.
251nfft = size(BF, 1);
252nl = numel(mask);
253R = fft(residual, nfft);
254o = 0; kk = 0; idx = 0; g = 0; red = 0;
255for 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
268end
269end
271function 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.
274nfft = 2048;
275m = min(numel(seg), nfft);
276s2 = seg(1:m);
277win = 0.5 - 0.5 * cos(2 * pi * (0:m-1)' / max(m - 1, 1));
278S = fft(s2 .* win, nfft);
279p = abs(S(1:nfft/2+1)).^2;
280hz = fs / nfft;
281edges = [0 100 200 400 800 1600 3200 6400];
282b = zeros(1, 7);
283for 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);
287end
288b = b - mean(b);
289kv = (0:numel(p)-1)' * hz;
290cen = sum(p .* kv) / (sum(p) + 1e-12);
291f = [b, log10(cen + 1)];
292end