Methods

What each stage computes, why that choice and not another, and the conditions under which it gives a wrong answer. Parameter names here are the keyword arguments on the API page.

Filtering#

The default band is 300–6000 Hz, fourth-order Butterworth, applied as second-order sections forwards and then backwards. The low corner removes the local field potential, which would otherwise wander the baseline by hundreds of microvolts and make a fixed threshold meaningless. The high corner is set by the sample rate as much as by the signal: at 30 kHz there is little above 6 kHz that is not noise.

The forward-backward pass cancels phase distortion, which matters because the next stage measures the position of a peak to better than one sample. A causal filter of this order shifts the peak by roughly 0.2 ms at 300 Hz, and worse, shifts it by a different amount depending on the spike's own spectrum — so two spikes from the same neuron would be aligned differently according to how much baseline noise rode along with each. The cost of the zero-phase pass is that it needs the whole chunk in hand, so this cannot run on a live stream.

The noise estimate#

The threshold is a multiple of the noise level, so the noise level has to be estimated from a recording that contains the very spikes you are trying to find. A plain standard deviation is pulled upwards by them: on the session in the tutorial it reads 9.1 µV against 7.4 µV for the robust estimate, which moves the threshold by 8 µV and silently discards the smallest units.

axonforge uses the median of the absolute value, scaled to a Gaussian:

sigma = median(|x|) / 0.6745

The constant is the inverse of the standard normal quantile at 0.75, so the expression returns the standard deviation for data that really is Gaussian, while ignoring the tails where the spikes live. It is computed per channel over the first span_s seconds by default, which assumes the noise is stationary — see the limitation below.

Threshold and dead time#

An event is recorded at the sample where the signal crosses k × sigma in the chosen polarity. The default k = 4.5 is a compromise that the validation page quantifies: at 4.0 recall improves by about two points and false positives roughly double; at 5.0 the smallest ground-truth unit disappears entirely.

After each crossing, further crossings on the same channel are suppressed for dead_time_ms. This is not a model of the refractory period — it is there because a single spike crosses the threshold several times on its way down and back. One millisecond is longer than any spike and shorter than any real refractory period, which is exactly the window you want.

Alignment#

Extraction cuts a window from pre_ms before to post_ms after the crossing, which is 60 samples at 30 kHz for the defaults. The crossing sample is a poor reference point, because where exactly the signal crosses a fixed level depends on amplitude: a larger spike of identical shape crosses earlier. Aligning on the crossing therefore maps amplitude differences onto time differences, and the clusters shear.

So the window is re-centred. In mode="peak" the three samples around the extremum are fitted with a parabola and the vertex is taken as the true peak time, giving a resolution finer than the 33 µs sample interval. In mode="min_slope" the steepest point of the falling edge is used instead, which is more stable for spikes whose peak is clipped or flat.

Features#

Aligned waveforms are projected onto their first principal components, three by default. This is a choice of convenience and convention rather than of optimality: the components are chosen to explain variance, and variance is not the same thing as the difference between units. A component that captures the drift of the baseline across the session will dominate a component that separates two cells of similar amplitude.

That is why explained_variance_ratio_ is worth printing. If the first component alone carries more than roughly 80 % of the variance, the decomposition has found amplitude and little else, and the clustering that follows will sort by distance from the electrode rather than by cell identity.

Clustering and model order#

A Gaussian mixture with full covariance is fitted for each k from 1 to max_k, restarts times from different random initialisations, and the fit with the best Bayesian information criterion is kept. Full covariance is the default because clusters in principal-component space are elongated along the amplitude direction far more often than they are round, and a spherical model splits one elongated cluster into two.

The criterion penalises parameters, so it resists the split that maximum likelihood alone would always prefer. It is not a truth oracle. On our synthetic sets the criterion picks the correct k in 17 of 20 repetitions at a signal-to-noise ratio of 6 and in 11 of 20 at a ratio of 3, and when it is wrong it over-splits rather than merges. Over-splitting is the more recoverable error — UnitSet.merge exists and is recorded in the provenance — which is part of why this criterion rather than a less conservative one.

Standing limitations#

• Stationary noise is assumed. One estimate per channel is made and reused for the whole recording. A session whose noise floor rises during the second half gets a threshold that is effectively lowered in sigma terms, and the event rate rises with it. Estimate per segment and pass sigma explicitly if this matters.
• No drift correction. Waveforms are treated as samples of a fixed set of shapes. Electrode movement violates this, and the symptom is a falling presence ratio rather than an error.
• Overlaps are not resolved. Two spikes within the dead time are one event. The overlap measurement shows the second spike is lost entirely below about 0.3 ms and half-recovered at 0.5 ms.
• Single-channel detection. A spike seen on four channels of a group produces up to four events, de-duplicated by time and channel group afterwards. Spatial detection across a dense probe is a different algorithm, and we have not written it.