An algorithm for spike detection

Source code link.

This is a matlab code for detecting spikes in a recording with multiple trials. These are the four steps invovled:

1. Low Pass Filtering used on the raw signal to smoothen out the recording noise. For a sampling rate of 20 kHz, a frequency of 500Hz for the low pass filter works for the purpose.

2. Second derivative based thresholding: To identify candidate peaks that can be marked as an action potential, we look at points in the smoothened signal that cross a second derivative threshold and mark the nearest local maxima.

3. Feature extraction: feature vector is 3-dimensional, populated by the minima of second derivative value reached in this region, absolute voltage value at the peak, and the amplitude of the spike measured from peak to the nearest of the two local minimas.

4. Clustering: using linkage based agglomerative clustering function of matlab.

Requirements

• Matlab R2019b or higher.

Inputs

• dataSet : a matrix with each row as one trial of the recording

• params

  • params.Fs --- Sampling frequency of the recording in Hz. (Default Value : 20000 Hz)

  • params.lowPassFreq --- Frequency of the low pass filter applied to filter recording noise in Hz. (Default Value : 500 Hz)

  • params.spikeSdThresh --- Number of standard deviations of second derivative values upto which candidate peaks are considered. (Default value of 0.5 is used)

  • params.clustSep --- Minimum value of the ratio (inter-cluster distance)/(intra-cluster distance 1 + intra-cluster distance 2) for considering trials together, else each trial is considered separately. (Default value of 0.7 is used)

  • params.maxclust --- matrix of size 1xnumTrial describing number of spike cluster to expect in each trial. (Default value of 2 used for each trial)

  • params.spikeDebug --- Flag for debugging, to be used for generating plots and choosing optimal parameter values. (Default value is 0)

Outputs

• spikeTimes --- a matrix with each row containing the time values of spike occurences (in seconds) for a trial.

Troubleshoot

The following strategies are helpful in cases of poor detection accuracy:

• params.spikeSdThresh

  • The default value of 0.5 times the standard deviation works well. In case you expect flatter action potential these parameter value can be lowered.

Example : The following recording shows an action potential that has a second derivative value that is 0.3 times the standard deviation. Also shown are corresponding changes in clustering.

• params.clustSep

  • The default ratio value of 0.7 is used. In case where trials are being considered together but one wants to consider them separately, increase this parameter.

Example: The following clusters are too close, and result in false negatives.

• params.maxClust

  • Default value of 2 is used, in case you need to increase true positive rate you can increase this parameter to 3,4 ,etc. according to expected number of spike clusters.
    Example : In the following example, true spikes get added to the cluster when maxClust is changed from 2 to 3.