function [times, labels] = onset_knn(x, fs, training)
% Detection and classification kept separate: events are found on a
% smoothed energy envelope (rise above a robust threshold), then each
% event is classified by distance-weighted k-nearest-neighbours over
% spectral-shape features of the training snippets. Features are level
% ratios, so the random per-event gain carries no information.
nTypes = numel(training.snips);
L = size(training.snips{1}, 1);
pre = training.pre;

% ---- detection ----------------------------------------------------------
w = round(0.004 * fs);
env = conv(x .* x, ones(w, 1) / w, 'same');
env = sqrt(env);
n = numel(env);
thr = max(8 * median(env), 0.05 * max(env));
refr = round(L / 2);

starts = zeros(0, 1);
i = 2;
while i < n
  if env(i) > thr && env(i-1) <= thr
    % crossing; find the local peak, then walk back to the foot of the rise
    j2 = min(n, i + refr - 1);
    [pk, rel] = max(env(i:j2));
    p = i + rel - 1;
    s = p;
    while s > 1 && env(s-1) < env(s) && env(s) > 0.05 * pk
      s = s - 1;
    end
    starts(end+1, 1) = s;
    % hold off until the envelope falls well below threshold (hysteresis),
    % so ripple in a long ringing tail is not reported again
    q = p;
    while q < n && env(q) > 0.5 * thr
      q = q + 1;
    end
    i = max(p + refr, q);
  else
    i = i + 1;
  end
end

% ---- training features --------------------------------------------------
nTrain = 0;
for k = 1:nTypes
  nTrain = nTrain + size(training.snips{k}, 2);
end
F = zeros(nTrain, 6);
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, :) = event_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 detected event --------------------------------------
kNN = min(5, nTrain);
times = zeros(0, 1);
labels = zeros(0, 1);
for e = 1:numel(starts)
  s = starts(e);
  seg = x(s:min(n, s + (L - pre) - 1));
  f = (event_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) / fs;
  labels(end+1, 1) = lab;
end
end

function f = event_features(seg, fs)
% Spectral shape of the event: four log band-energy ratios (mean removed,
% so absolute level cancels), log spectral centroid, and zero-crossing
% rate over the first 25 ms.
nfft = 1024;
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)).^2;
hz = fs / nfft;
edges = [0 200 800 2500 8000];
b = zeros(1, 4);
for q = 1:4
  k0 = max(2, floor(edges(q) / hz) + 1);
  k1 = min(nfft/2, floor(edges(q+1) / hz));
  b(q) = log10(sum(p(k0:k1)) + 1e-12);
end
b = b - mean(b);
kv = (0:nfft/2-1)' * hz;
cen = sum(p .* kv) / (sum(p) + 1e-12);
zn = min(numel(seg), round(0.025 * fs));
zc = sum(abs(diff(seg(1:zn) >= 0)));
f = [b, log10(cen + 1), zc / max(zn, 1) * fs / 1000];
end