1function [times, labels] = silence_gap_knn(x, fs, training)
2% Segment on quiet gaps, then classify with kNN over spectral shape.
3%
4% Ported from the audio-spike-sorting benchmark, where it was developed
5% against synthesized recordings whose background is digitally silent.
6% One adaptation for live microphone takes: the segmentation threshold
7% is raised from a fraction of the recording peak to a robust noise
8% floor estimate (a multiple of the median absolute sample), since a
9% mic's background is never exactly zero. On digitally silent audio the
10% median is zero and the original threshold takes over, so benchmark
11% behavior is unchanged.
12%
13% Detection: samples above the threshold are grouped into regions, gaps
14% of a few milliseconds inside a region are closed, sub-10 ms glitches
15% are dropped, and a region much longer than one snippet is split into
16% snippet-length events, which handles back-to-back hits. The event time
17% is the region start plus the pre-onset offset. Each event is then
18% classified by distance-weighted k-nearest-neighbours over
19% level-invariant spectral features of the post-onset segment.
20%
21% Overlapping events remain the limitation: simultaneous hits merge into
22% one region and only the length-based split separates them.
23nTypes = numel(training.snips);
24L = size(training.snips{1}, 1);
25pre = training.pre;
26n = numel(x);
28% ---- detection: contiguous above-noise regions ---------------------------
29eps0 = max(1e-4 * max(abs(x)), 6 * median(abs(x)));
30gapClose = round(0.002 * fs);
31minDur = round(0.010 * fs);
33starts = zeros(0, 1);
34ends = zeros(0, 1);
35i = 1;
36while i <= n
37 if abs(x(i)) > eps0
38 s = i;
39 e = i;
40 j = i + 1;
41 quiet = 0;
42 while j <= n && quiet <= gapClose
43 if abs(x(j)) > eps0
44 e = j;
45 quiet = 0;
46 else
47 quiet = quiet + 1;
48 end
49 j = j + 1;
50 end
51 if e - s + 1 >= minDur
52 starts(end+1, 1) = s;
53 ends(end+1, 1) = e;
54 end
55 i = e + gapClose + 1;
56 else
57 i = i + 1;
58 end
59end
61% split regions covering several back-to-back snippets
62ev = zeros(0, 1);
63for r = 1:numel(starts)
64 nev = max(1, round((ends(r) - starts(r) + 1) / L));
65 for q = 1:nev
66 ev(end+1, 1) = starts(r) + (q - 1) * L;
67 end
68end
70% ---- training features --------------------------------------------------
71nTrain = 0;
72for k = 1:nTypes
73 nTrain = nTrain + size(training.snips{k}, 2);
74end
75F = zeros(nTrain, 7);
76Fl = zeros(nTrain, 1);
77r = 0;
78for k = 1:nTypes
79 S = training.snips{k};
80 for j = 1:size(S, 2)
81 r = r + 1;
82 F(r, :) = spectral_features(S(pre+1:end, j), fs);
83 Fl(r) = k;
84 end
85end
86mu = mean(F, 1);
87sd = std(F, 0, 1) + 1e-9;
88Fz = (F - repmat(mu, nTrain, 1)) ./ repmat(sd, nTrain, 1);
90% ---- classify each event ------------------------------------------------
91kNN = min(5, nTrain);
92times = zeros(0, 1);
93labels = zeros(0, 1);
94for e = 1:numel(ev)
95 s = ev(e);
96 seg = x(min(n, s + pre):min(n, s + L - 1));
97 f = (spectral_features(seg, fs) - mu) ./ sd;
98 d = sqrt(sum((Fz - repmat(f, nTrain, 1)).^2, 2));
99 [ds, di] = sort(d);
100 votes = zeros(nTypes, 1);
101 for v = 1:kNN
102 votes(Fl(di(v))) = votes(Fl(di(v))) + 1 / (ds(v) + 1e-6);
103 end
104 [~, lab] = max(votes);
105 times(end+1, 1) = (s - 1 + pre) / fs;
106 labels(end+1, 1) = lab;
107end
108end
110function f = spectral_features(seg, fs)
111% Six log band-energy ratios (mean removed, so absolute level cancels)
112% and the log spectral centroid.
113nfft = 2048;
114m = min(numel(seg), nfft);
115s2 = seg(1:m);
116win = 0.5 - 0.5 * cos(2 * pi * (0:m-1)' / max(m - 1, 1));
117S = fft(s2 .* win, nfft);
118p = abs(S(1:nfft/2+1)).^2;
119hz = fs / nfft;
120edges = [0 200 500 1000 2000 4000 8000];
121b = zeros(1, 6);
122for q = 1:6
123 k0 = max(2, floor(edges(q) / hz) + 1);
124 k1 = min(nfft/2 + 1, floor(edges(q+1) / hz));
125 b(q) = log10(sum(p(k0:k1)) + 1e-12);
126end
127b = b - mean(b);
128kv = (0:numel(p)-1)' * hz;
129cen = sum(p .* kv) / (sum(p) + 1e-12);
130f = [b, log10(cen + 1)];
131end