7d7cc56Add overlap pursuit classifier and update composerJeremy Magland 1function [times, labels] = overlap_pursuit(x, fs, training)
2% Quiet-gap segmentation plus per-block matching pursuit for overlaps.
3%
4% Ported from the audio-spike-sorting benchmark. Activity above a noise
5% threshold splits into regions; a region about one snippet long holds
6% a single hit, classified by the spectral kNN of silence_gap_knn. A
7% longer region holds overlapping hits and is decomposed by matching
8% pursuit: the region start is an anchor, fitted by least squares
9% against every training snippet (the best instance and gain, refined
10% by coordinate-descent sweeps that add each hit back and refit it
11% against the residual), and further hits are found by scanning the
12% residual with FFT cross-correlation against the whole bank. A
13% candidate is accepted while its gain and energy reduction are
14% plausible, unless a hit of the same voice already sits within a few
15% hundred samples, since instance mismatch of an already-fitted hit is
16% the usual source of such candidates; that veto also merges genuine
17% same-voice hits closer than ~16 ms into one.
18%
19% Adaptations for live microphone takes. On the benchmark, a region's
20% extent is the exact union of snippet spans, which justifies a second
21% anchor at the region end minus one snippet and a coverage rule that
22% forces an event onto every sample. On a live take a hit's ring-out
23% ends wherever it sinks into the noise floor, so both inferences would
24% invent hits out of decay tails; this port keeps only the region-start
25% anchor and finds every further hit by the threshold-gated correlation
26% pursuit. Concretely: the segmentation threshold is raised to a robust
27% noise floor estimate (a multiple of the median absolute sample; on
28% digitally silent audio the median is zero and the original threshold
29% takes over); a region up to 5 ms longer than one snippet is a single
30% hit rather than a near-coincident pair; each merged region is padded
31% with one snippet length of trailing context so a hit near the region
32% end keeps a full span, with onsets confined to the region proper; and
33% the pursuit's energy-reduction floor is a fraction of the weakest
34% training snippet's energy rather than an absolute number, so it does
35% not depend on the recording level.
36nTypes = numel(training.snips);
37L = size(training.snips{1}, 1);
38pre = training.pre;
39n = numel(x);
41% tuned constants
42gmin = 0.35; % smallest plausible fitted gain for a pursuit hit
43prox = round(0.0156 * fs); % same-voice exclusion radius, ~250 samples
44maxDup = 6; % pursuit iterations per merged region
46% ---- training bank ------------------------------------------------------
47nTrain = 0;
48for k = 1:nTypes
49 nTrain = nTrain + size(training.snips{k}, 2);
50end
51B = zeros(L, nTrain);
52bt = zeros(nTrain, 1);
53c = 0;
54for k = 1:nTypes
55 S = training.snips{k};
56 for j = 1:size(S, 2)
57 c = c + 1;
58 B(:, c) = S(:, j);
59 bt(c) = k;
60 end
61end
62n2 = zeros(nTrain, 1);
63for j = 1:nTrain
64 n2(j) = sum(B(:, j) .* B(:, j));
65end
66% smallest energy reduction for a pursuit hit, relative to the weakest
67% training snippet so the recording level cancels
68redmin = 0.03 * min(n2);
70% ---- kNN training features for single-event blocks ----------------------
71F = zeros(nTrain, 8);
72for j = 1:nTrain
73 F(j, :) = spectral_features(B(pre+1:end, j), fs);
74end
75mu = mean(F, 1);
76sd = std(F, 0, 1) + 1e-9;
77Fz = (F - repmat(mu, nTrain, 1)) ./ repmat(sd, nTrain, 1);
79% ---- silence segmentation ------------------------------------------------
80eps0 = max(1e-4 * max(abs(x)), 6 * median(abs(x)));
81gapClose = round(0.002 * fs);
82minDur = round(0.010 * fs);
83bs = zeros(0, 1);
84be = zeros(0, 1);
85i = 1;
86while i <= n
87 if abs(x(i)) > eps0
88 s = i;
89 e = i;
90 j = i + 1;
91 quiet = 0;
92 while j <= n && quiet <= gapClose
93 if abs(x(j)) > eps0
94 e = j;
95 quiet = 0;
96 else
97 quiet = quiet + 1;
98 end
99 j = j + 1;
100 end
101 if e - s + 1 >= minDur
102 bs(end+1, 1) = s;
103 be(end+1, 1) = e;
104 end
105 i = e + gapClose + 1;
106 else
107 i = i + 1;
108 end
109end
111% ---- per-block detection and classification ------------------------------
112times = zeros(0, 1);
113labels = zeros(0, 1);
114for bi = 1:numel(bs)
115 s = bs(bi);
116 D = be(bi) - s + 1;
117 if D <= L + round(0.005 * fs)
118 % single hit: onset at the block start, spectral kNN label
119 seg = x(s+pre : min(n, s+L-1));
120 f = (spectral_features(seg, fs) - mu) ./ sd;
121 d = sqrt(sum((Fz - repmat(f, nTrain, 1)).^2, 2));
122 [ds, di] = sort(d);
123 kNN = min(5, nTrain);
124 votes = zeros(nTypes, 1);
125 for v = 1:kNN
126 votes(bt(di(v))) = votes(bt(di(v))) + 1 / (ds(v) + 1e-6);
127 end
128 [~, lab] = max(votes);
129 times(end+1, 1) = (s - 1 + pre) / fs;
130 labels(end+1, 1) = lab;
131 else
132 % pad with trailing ring context so a hit near the region end keeps
133 % a full snippet span; onsets stay confined to the region proper
134 blk = x(s : min(n, s + D - 1 + L));
135 [evOn, evType] = decompose(blk, B, bt, n2, L, gmin, redmin, prox, maxDup, D);
136 for q = 1:numel(evOn)
137 times(end+1, 1) = (s - 1 + evOn(q) - 1 + pre) / fs;
138 labels(end+1, 1) = evType(q);
139 end
140 end
141end
142end
144function [evOn, evType] = decompose(blk, B, bt, n2, L, gmin, redmin, prox, maxDup, nOn)
145% Matching pursuit over one merged region. blk is the region plus up to
146% one snippet length of trailing context; nOn is the length of the
147% region proper, the only stretch where an onset may be placed.
148D = numel(blk);
149nl = D - L + 1;
150nTrain = numel(bt);
152% conjugate template spectra for the FFT correlation scans
153nfft = 2^nextpow2(D);
154BF = zeros(nfft, nTrain) + 0i;
155for j = 1:nTrain
156 BF(:, j) = conj(fft(B(:, j), nfft));
157end
159% events as parallel arrays, seeded with the region-start anchor
160evOn = [1];
161evType = zeros(1, 1);
162evIdx = zeros(1, 1);
163evGain = zeros(1, 1);
164residual = blk;
166[residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);
168% correlation pursuit for further hits in the region
169mask = zeros(nl, 1);
170mask(1 : min(nOn, nl)) = 1;
171for it = 1:maxDup
172 [o, kk, idx, g, red] = scan_best(residual, BF, bt, n2, L, mask);
173 if idx == 0 || red < redmin || g < gmin
174 break;
175 end
176 clash = 0;
177 for q = 1:numel(evOn)
178 if evOn(q) == o || (abs(evOn(q) - o) < prox && evType(q) == kk)
179 clash = 1;
180 end
181 end
182 if clash
183 break;
184 end
185 residual(o:o+L-1) = residual(o:o+L-1) - g * B(:, idx);
186 evOn(end+1, 1) = o;
187 evType(end+1, 1) = kk;
188 evIdx(end+1, 1) = idx;
189 evGain(end+1, 1) = g;
190 [residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);
191end
192end
194function [residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, nsw)
195% Coordinate descent: add each event back and refit its type/instance/gain.
196[~, ord] = sort(evOn);
197for it = 1:nsw
198 for z = 1:numel(ord)
199 q = ord(z);
200 o = evOn(q);
201 if evIdx(q) > 0
202 residual(o:o+L-1) = residual(o:o+L-1) + evGain(q) * B(:, evIdx(q));
203 end
204 seg = residual(o:o+L-1);
205 bestRed = 0;
206 bestJ = 0;
207 bestG = 0;
208 for j = 1:numel(bt)
209 cj = sum(seg .* B(:, j));
210 if cj > 0
211 red = cj * cj / n2(j);
212 if red > bestRed
213 bestRed = red;
214 bestJ = j;
215 bestG = cj / n2(j);
216 end
217 end
218 end
219 if bestJ > 0
220 residual(o:o+L-1) = residual(o:o+L-1) - bestG * B(:, bestJ);
221 evType(q) = bt(bestJ);
222 evIdx(q) = bestJ;
223 evGain(q) = bestG;
224 end
225 end
226end
227end
229function [o, kk, idx, g, red] = scan_best(residual, BF, bt, n2, L, mask)
230% Best (onset, template) over masked lags by FFT cross-correlation;
231% reduction is the residual energy removed by the least-squares gain.
232nfft = size(BF, 1);
233nl = numel(mask);
234R = fft(residual, nfft);
235o = 0; kk = 0; idx = 0; g = 0; red = 0;
236for j = 1:numel(bt)
237 cf = real(ifft(R .* BF(:, j)));
238 cl = cf(1:nl);
239 rv = (cl .* cl) / n2(j);
240 rv = rv .* (cl > 0) .* mask;
241 [m, mi] = max(rv);
242 if m > red
243 red = m;
244 o = mi;
245 kk = bt(j);
246 idx = j;
247 g = cl(mi) / n2(j);
248 end
249end
250end
252function f = spectral_features(seg, fs)
253% Seven log band-energy ratios (mean removed, level cancels) and the
254% log spectral centroid; the extra low bands split the two bass types.
255nfft = 2048;
256m = min(numel(seg), nfft);
257s2 = seg(1:m);
258win = 0.5 - 0.5 * cos(2 * pi * (0:m-1)' / max(m - 1, 1));
259S = fft(s2 .* win, nfft);
260p = abs(S(1:nfft/2+1)).^2;
261hz = fs / nfft;
262edges = [0 100 200 400 800 1600 3200 6400];
263b = zeros(1, 7);
264for q = 1:7
265 k0 = max(2, floor(edges(q) / hz) + 1);
266 k1 = min(nfft/2 + 1, floor(edges(q+1) / hz));
267 b(q) = log10(sum(p(k0:k1)) + 1e-12);
268end
269b = b - mean(b);
270kv = (0:numel(p)-1)' * hz;
271cen = sum(p .* kv) / (sum(p) + 1e-12);
272f = [b, log10(cen + 1)];
273end