This is a very basic but very fast window discriminator for detecting and extracting spikes in a time series. It was developed for analyzing extracellular neural recordings, but also works with intracellular data and probably many other kinds of time series.
Here’s how it works:
detection diagram
The algorithm iterates through the time series. When the signal crosses the threshold (1) going away from zero, the algorithm then looks for a peak (2) that occurs within some number of samples from the threshold crossing. The number of samples can be adjusted to filter out broad spikes that are likely to be artifacts. If a peak occurs, its sample index is added to an array. These times can be used as-is, or they can be used to extract samples to either side of the peak for further analysis (e.g. spike sorting).
The algorithm uses a streaming pattern (i.e., it processes chunks of data and keeps its state between chunks), so it’s suitable for realtime operations. Many signals of interest will require highpass filtering to remove slow variations.
The algorithm is written in cython. You can get a python package from PyPI:
pip install quickspikes
Or to build from a copy of the repository:
pip install .
To detect peaks, you instantiate the detector with parameters that match the events you want to detect, and then send the detector chunks of data. For example, an extracellular recording at 20 kHz stored in 16-bit integers may have a noise floor around 2000, and the spikes will be on the order of 20 samples wide:
import quickspikes as qs
det = qs.detector(1000, 30)
times = det.send(samples)You can continue sending chunks of data by calling send(). The
detector will keep its state between calls, so you can detect spikes
that occur on chunk boundaries. For example, if you’re receiving data
from some kind of generator, you could use a pattern like this:
for chunk in my_data_generator():
times = det.send(chunk)
# process timesConversely, if the data are not contiguous, you should reinitialize the detector for each chunk.
You can adjust the detector’s threshold at any point, for example to compensate for shifts in the mean and standard deviation of the signal:
reldet = qs.detector(2.5, 30)
reldet.scale_thresh(samples.mean(), samples.std())
times = reldet.send(samples)To detect negative-going events, you’ll need to invert the signal.
There are also some functions you can use to extract and align spike
waveforms. Given a list of times returned from the detector.send()
method, to extract waveforms starting 30 samples before the peak and
ending 270 samples after:
f_times = qs.filter_times(times, 30, samples.size - 270)
spikes = qs.peaks(samples, f_times, 30, 270)
times, aligned = qs.realign_spikes(f_times, spikes, upsample=3, jitter=4)Note that the list of event times may need to be filtered to avoid
trying to access data points outside the bounds of the input time
series. If you care about these events, you’ll need to pad your input
signal. The realign_spikes function uses a sinc-based resampling to
more accurately locate the peak of the event.
There is also a reference copy of an ANSI C implementation and an
f2py wrapper in f2py/. This algorithm is slightly less efficient
and flexible, but may give better results if included directly in a C
codebase.
There are some specialized tools for working with intracellular data. Intracellular recordings present some special challenges for detecting spikes. The data are not centered, and spike waveforms can change over the course of stimulation. The approach used here is to find the largest spike in the recording, (typically the first) and locate the onset of the spike based on when the derivative begins to rapidly increase. The peak and onset are used to set a threshold for extracting subsequent spikes.
The following is an example for a recording at 50 kHz of a neuron with spikes that have a rise time of about 1 ms (50 samples). Spikes waveforms will start 7 ms before the peak and end 40 ms after, and will be trimmed to avoid any overlap with subsequent spikes.
from quickspikes.intracellular import SpikeFinder
detector = SpikeFinder(n_rise=50, n_before=350, n_after=2000)
detector.calculate_threshold(samples)
times, spikes = zip(*detector.extract_spikes(samples, min_amplitude=10))Free for use under the terms of the GNU General Public License. See [[COPYING]] for details.
If you use this code in an academic work, citations are appreciated. There is no methods paper describing the algorithm, but the most relevant reference is:
C. D. Meliza and D. Margoliash (2012). Emergence of selectivity and tolerance in the avian auditory cortex. Journal of Neuroscience, doi:10.1523/JNEUROSCI.0845-12.2012