function [times, labels] = silence_gap_knn(x, fs, training)
% Segment on quiet gaps, then classify with kNN over spectral shape.
%
% Ported from the audio-spike-sorting benchmark, where it was developed
% against synthesized recordings whose background is digitally silent.
% One adaptation for live microphone takes: the segmentation threshold
% is raised from a fraction of the recording peak to a robust noise
% floor estimate (a multiple of the median absolute sample), since a
% mic's background is never exactly zero. On digitally silent audio the
% median is zero and the original threshold takes over, so benchmark
% behavior is unchanged.
%
% Detection: samples above the threshold are grouped into regions, gaps
% of a few milliseconds inside a region are closed, sub-10 ms glitches
% are dropped, and a region much longer than one snippet is split into
% snippet-length events, which handles back-to-back hits. The event time
% is the region start plus the pre-onset offset. Each event is then
% classified by distance-weighted k-nearest-neighbours over
% level-invariant spectral features of the post-onset segment.
%
% Overlapping events remain the limitation: simultaneous hits merge into
% one region and only the length-based split separates them.
nTypes = numel(training.snips);
L = size(training.snips{1}, 1);
pre = training.pre;
n = numel(x);

% ---- detection: contiguous above-noise regions ---------------------------
eps0 = max(1e-4 * max(abs(x)), 6 * median(abs(x)));
gapClose = round(0.002 * fs);
minDur = round(0.010 * fs);

starts = zeros(0, 1);
ends = zeros(0, 1);
i = 1;
while i <= n
  if abs(x(i)) > eps0
    s = i;
    e = i;
    j = i + 1;
    quiet = 0;
    while j <= n && quiet <= gapClose
      if abs(x(j)) > eps0
        e = j;
        quiet = 0;
      else
        quiet = quiet + 1;
      end
      j = j + 1;
    end
    if e - s + 1 >= minDur
      starts(end+1, 1) = s;
      ends(end+1, 1) = e;
    end
    i = e + gapClose + 1;
  else
    i = i + 1;
  end
end

% split regions covering several back-to-back snippets
ev = zeros(0, 1);
for r = 1:numel(starts)
  nev = max(1, round((ends(r) - starts(r) + 1) / L));
  for q = 1:nev
    ev(end+1, 1) = starts(r) + (q - 1) * L;
  end
end

% ---- training features --------------------------------------------------
nTrain = 0;
for k = 1:nTypes
  nTrain = nTrain + size(training.snips{k}, 2);
end
F = zeros(nTrain, 7);
Fl = zeros(nTrain, 1);
r = 0;
for k = 1:nTypes
  S = training.snips{k};
  for j = 1:size(S, 2)
    r = r + 1;
    F(r, :) = spectral_features(S(pre+1:end, j), fs);
    Fl(r) = k;
  end
end
mu = mean(F, 1);
sd = std(F, 0, 1) + 1e-9;
Fz = (F - repmat(mu, nTrain, 1)) ./ repmat(sd, nTrain, 1);

% ---- classify each event ------------------------------------------------
kNN = min(5, nTrain);
times = zeros(0, 1);
labels = zeros(0, 1);
for e = 1:numel(ev)
  s = ev(e);
  seg = x(min(n, s + pre):min(n, s + L - 1));
  f = (spectral_features(seg, fs) - mu) ./ sd;
  d = sqrt(sum((Fz - repmat(f, nTrain, 1)).^2, 2));
  [ds, di] = sort(d);
  votes = zeros(nTypes, 1);
  for v = 1:kNN
    votes(Fl(di(v))) = votes(Fl(di(v))) + 1 / (ds(v) + 1e-6);
  end
  [~, lab] = max(votes);
  times(end+1, 1) = (s - 1 + pre) / fs;
  labels(end+1, 1) = lab;
end
end

function f = spectral_features(seg, fs)
% Six log band-energy ratios (mean removed, so absolute level cancels)
% and the log spectral centroid.
nfft = 2048;
m = min(numel(seg), nfft);
s2 = seg(1:m);
win = 0.5 - 0.5 * cos(2 * pi * (0:m-1)' / max(m - 1, 1));
S = fft(s2 .* win, nfft);
p = abs(S(1:nfft/2+1)).^2;
hz = fs / nfft;
edges = [0 200 500 1000 2000 4000 8000];
b = zeros(1, 6);
for q = 1:6
  k0 = max(2, floor(edges(q) / hz) + 1);
  k1 = min(nfft/2 + 1, floor(edges(q+1) / hz));
  b(q) = log10(sum(p(k0:k1)) + 1e-12);
end
b = b - mean(b);
kv = (0:numel(p)-1)' * hz;
cen = sum(p .* kv) / (sum(p) + 1e-12);
f = [b, log10(cen + 1)];
end