/ proof-of-concept / audio-spike-sorting
Sign in
proof-of-concept / audio-spike-sorting
audio-spike-sorting / src / algorithms / silence_gap_knn.m
128 lines · 3.6 KBBlameHistoryRaw
1function [times, labels] = silence_gap_knn(x, fs, training)
2% Segment on silence, then classify with kNN over spectral shape.
3%
4% The synthesized recordings place each event snippet on a digitally
5% silent background, so an event is exactly one contiguous run of
6% nonzero samples, one snippet long. Detection therefore reduces to
7% finding those runs: samples above a tiny threshold (relative to the
8% recording peak) are grouped into regions, gaps of a few milliseconds
9% inside a region are closed, sub-10 ms glitches are dropped, and a
10% region much longer than one snippet is split into snippet-length
11% events, which handles back-to-back placements. The event time is the
12% region start plus the pre-onset offset, which lands on the ground
13% truth exactly. Each event is then classified by distance-weighted
14% k-nearest-neighbours over level-invariant spectral features of the
15% post-onset segment, as in onset_knn.
17% The silence assumption is the limitation: on a recording with real
18% background noise, or with overlapping events, the segmentation would
19% fail and an envelope or matched-filter detector would be needed.
20nTypes = numel(training.snips);
21L = size(training.snips{1}, 1);
22pre = training.pre;
23n = numel(x);
25% ---- detection: contiguous nonzero regions ------------------------------
26eps0 = 1e-4 * max(abs(x));
27gapClose = round(0.002 * fs);
28minDur = round(0.010 * fs);
30starts = zeros(0, 1);
31ends = zeros(0, 1);
32i = 1;
33while i <= n
34 if abs(x(i)) > eps0
35 s = i;
36 e = i;
37 j = i + 1;
38 quiet = 0;
39 while j <= n && quiet <= gapClose
40 if abs(x(j)) > eps0
41 e = j;
42 quiet = 0;
43 else
44 quiet = quiet + 1;
45 end
46 j = j + 1;
47 end
48 if e - s + 1 >= minDur
49 starts(end+1, 1) = s;
50 ends(end+1, 1) = e;
51 end
52 i = e + gapClose + 1;
53 else
54 i = i + 1;
55 end
56end
58% split regions covering several back-to-back snippets
59ev = zeros(0, 1);
60for r = 1:numel(starts)
61 nev = max(1, round((ends(r) - starts(r) + 1) / L));
62 for q = 1:nev
63 ev(end+1, 1) = starts(r) + (q - 1) * L;
64 end
65end
67% ---- training features --------------------------------------------------
68nTrain = 0;
69for k = 1:nTypes
70 nTrain = nTrain + size(training.snips{k}, 2);
71end
72F = zeros(nTrain, 7);
73Fl = zeros(nTrain, 1);
74r = 0;
75for k = 1:nTypes
76 S = training.snips{k};
77 for j = 1:size(S, 2)
78 r = r + 1;
79 F(r, :) = spectral_features(S(pre+1:end, j), fs);
80 Fl(r) = k;
81 end
82end
83mu = mean(F, 1);
84sd = std(F, 0, 1) + 1e-9;
85Fz = (F - repmat(mu, nTrain, 1)) ./ repmat(sd, nTrain, 1);
87% ---- classify each event ------------------------------------------------
88kNN = min(5, nTrain);
89times = zeros(0, 1);
90labels = zeros(0, 1);
91for e = 1:numel(ev)
92 s = ev(e);
93 seg = x(min(n, s + pre):min(n, s + L - 1));
94 f = (spectral_features(seg, fs) - mu) ./ sd;
95 d = sqrt(sum((Fz - repmat(f, nTrain, 1)).^2, 2));
96 [ds, di] = sort(d);
97 votes = zeros(nTypes, 1);
98 for v = 1:kNN
99 votes(Fl(di(v))) = votes(Fl(di(v))) + 1 / (ds(v) + 1e-6);
100 end
101 [~, lab] = max(votes);
102 times(end+1, 1) = (s - 1 + pre) / fs;
103 labels(end+1, 1) = lab;
104end
105end
107function f = spectral_features(seg, fs)
108% Six log band-energy ratios (mean removed, so absolute level cancels)
109% and the log spectral centroid.
110nfft = 2048;
111m = min(numel(seg), nfft);
112s2 = seg(1:m);
113win = 0.5 - 0.5 * cos(2 * pi * (0:m-1)' / max(m - 1, 1));
114S = fft(s2 .* win, nfft);
115p = abs(S(1:nfft/2+1)).^2;
116hz = fs / nfft;
117edges = [0 200 500 1000 2000 4000 8000];
118b = zeros(1, 6);
119for q = 1:6
120 k0 = max(2, floor(edges(q) / hz) + 1);
121 k1 = min(nfft/2 + 1, floor(edges(q+1) / hz));
122 b(q) = log10(sum(p(k0:k1)) + 1e-12);
123end
124b = b - mean(b);
125kv = (0:numel(p)-1)' * hz;
126cen = sum(p .* kv) / (sum(p) + 1e-12);
127f = [b, log10(cen + 1)];
128end
↑↓ move↵openescclose