Skip to content
dsp-kitchen

From spikes to units: snippets, features, clusters

Each detected spike becomes a snippet: a couple of milliseconds around it, on the channels near it on the probe. Snippets of one unit look alike (same shape, same channels), so they sit together once each is reduced to a few numbers: principal components of the snippets. A clustering of those points proposes the units. Here the truth is known, so the proposal can be checked: each detection is labelled by the true spike it matches.

Two clusterings: k-means (you say how many units) and HDBSCAN (finds dense groups by itself and calls the rest noise, -1).

import numpy as np


import dspk
from dspk.synapse import clustering, detection, extraction, io
from plating import signals


s = signals.synthetic(channels=8, seconds=4.0, units=3, rate_hz=15.0, noise_uv=10.0)
band = dspk.ops.filter.iir(s.data, s.sample_rate, "bandpass", [300.0, 6000.0])
probe = io.Probe.from_positions([[0, 20 * c] for c in range(s.channels)])
spikes = detection.deduplicate(detection.threshold(band, s.sample_rate), probe, s.sample_rate)


# Snippets on every channel (k = all), rows put back in channel order so features line up
snip = extraction.snippets(band, spikes, probe, s.sample_rate, k=s.channels, pre_ms=0.6, post_ms=1.0)
order = np.argsort(snip["channels"], axis=1)
waves = np.take_along_axis(snip["waveforms"], order[:, :, None], axis=1)  # (spikes, channels, samples)
flat = waves.reshape(waves.shape[0], -1)


# Principal components: dspk's PCA takes (dimensions, observations)
pca = dspk.ops.linalg.PCA.fit(flat.T, 3)
pcs = pca.transform(flat.T).T  # (spikes, 3)
print(f"{flat.shape[0]} snippets of {flat.shape[1]} values → 3 PCs ({pca.explained_variance_ratio.sum():.0%} of the variance)")


# The truth of each detection: the unit of the nearest true spike
offset = int(np.argmin(s.templates[0])) - s.templates.shape[1] // 2
truth_t = s.spike_samples.astype(np.int64) + offset
t = spikes["samples"][snip["kept"]].astype(np.int64)
nearest = np.abs(t[:, None] - truth_t[None, :]).argmin(axis=1)
truth = s.spike_units[nearest]


km = clustering.kmeans(pcs, s.params["units"])
hd = clustering.hdbscan(pcs, min_cluster_size=10)


def agreement(labels: np.ndarray) -> float:
    """Spikes whose cluster's majority unit is their own unit."""
    ok = 0
    for c in np.unique(labels[labels >= 0]):
        members = truth[labels == c]
        ok += np.bincount(members).max()
    return ok / labels.size


print(f"k-means agreement {agreement(km['labels']):.0%} · HDBSCAN agreement {agreement(hd):.0%} ({np.unique(hd[hd >= 0]).size} clusters, {np.mean(hd < 0):.0%} noise)")

From spikes to units: snippets, features, clusters

The whole script is playground/08_detection/02_features_and_clusters.py. Run it with uv run python playground/08_detection/02_features_and_clusters.py.