6th Annual International IEEE EMBS Conference on Neural Engineering San Diego, California, 6 - 8 November, 2013
Automatic Selection of Neuronal Spike Detection Threshold via Smoothed Teager Energy Histogram Andrew Keong Ng*, Member, IEEE, Kai Keng Ang, and Cuntai Guan, Senior Member, IEEE
Abstract— Spike detection is a prerequisite to analyzing neuronal activity. While simple spike detectors are favorable for hardware implementation, manual setting of spike detection threshold can be tedious and time-consuming, especially for extracellular recordings of multiple neuronal activity. This paper, therefore, investigates and proposes an automatic threshold selection using smoothed Teager energy histogram (STEH), with consideration of signal prewhitening, histogram bin width, and histogram equalization. Results from spikes with signal-to-noise ratio = 0.9-2.6 dB reveals that (1) prewhitening of neural signals can enhance true detection rate (TDR) by 1020% at constant false alarm (FA) ranging 3-12 spikes/s; (2) Freedman-Diaconis choice of STEH bin width delivers higher TDR (1.52 ± 1.41%) and FA (0.48 ± 0.25 spikes/s) than squareroot choice; and (3) histogram equalization can raise average TDR by 2.84% and FA by 1.01 spikes/s. Thresholds determined by STEH with signal prewhitening, Freedman-Diaconis choice or square-root choice of bin width, and histogram equalization fall around knee point of receiver operating characteristic curve, yielding average TDR = 87.88% and FA = 1.82 spikes/s for Freedman-Diaconis choice, and TDR = 86.49% and FA = 1.41 spikes/s for square-root choice of STEH bin width.
I. INTRODUCTION In an effort to develop high-performance brain-computer interfaces and neuroprosthetic devices, extracellular recordings of neuronal activity have been extensively used in neuroscience research to understand how the brain processes information . Spike detection, the process of identifying neuronal action potentials (spikes) immersed in background noise, is a prerequisite to analyzing and interpreting the recorded neural signals [2,3]. Furthermore, as the process outputs only a sparse collection of spike waveforms, it reduces transmission data rate per channel and increases number of channels in a wireless multichannel extracellular recording system [4,5]. Spike detection comprises of two main steps. First, it preprocesses neural signals to remove unwanted noise and intensify spikes relative to background noise. Second, it applies a detection threshold to extract spikes from the preprocessed signals. During the days with no digital computers, spikes were detected through a simple voltage trigger with a voltage threshold set by the user. When a voltage signal crossed the threshold, a pulse would be generated to signify the presence of a spike, and the spike waveform would be extracted [2,3]. While the amplitude thresholding is simple and easy to implement, its *A.K. Ng, K.K. Ang, and C. Guan are with the Institute for Infocomm Research, Agency for Science, Technology and Research (A*STAR), 1 Fusionopolis Way, #21-01 Connexis (South Tower), Singapore 138632, Singapore (e-mail: [email protected]
978-1-4673-1969-0/13/$31.00 ©2013 IEEE
performance deteriorates with decreasing signal-to-noise ratio (SNR) [4,5]. Energy thresholding is another technique for spike detection. It considers the instantaneous energy of neural signals by squaring the signal amplitude or by using the product of signal amplitude and frequency derived from nonlinear Teager energy operator (TEO) . Other spike detectors driven by stationary wavelet transformation and template matching were also proposed, but they have high computational cost, similarity assumption between spike shapes and wavelet bases, or requirement of a priori information about spike shapes, which may not be suited for hardware implementation in a real environment [3‒5]. Amplitude- and energy-driven spike detectors are usually preferred because they are computationally simple, easy to implement, and can deliver relatively good performance [4,5]. However, in situations involving extracellular recordings of activity across multiple neurons, manual setting of amplitude or energy threshold can be tedious and time-consuming as the threshold is determined through trial and error by the user. Moreover, being subjective in nature, the threshold can have large inter-observer variability, which may lead to poor reliability and consistency. Hence, an automatic threshold selection is more realistic and desirable. To optimize spike detection threshold selection, Rizk et al.  compared four noise estimators under different noise levels and suggested mean deviation operator as the most appropriate noise estimator. Conversely, Semmaoui et al.  proposed an adaptive threshold based on the relationship between low-order statistics of noise at a smoothed TEO (STEO) detector input and the mean and standard deviation of the detector output; the optimal multiplier value was found to lie between 5.5 and 7.5. Using the STEO detector, Malarvili et al.  chose the second bin value after the highest bin of smoothed Teager energy histogram (STEH) as the detection threshold, but they gave little explanation and information about the STEH. This paper further investigates the selection of spike detection threshold via STEH by considering (1) signal prewhitening, (2) histogram bin width, and (3) histogram equalization, which were not considered earlier, to the best of our knowledge. In doing so, we hope to provide more insight into STEH as an alternative means to autonomously select threshold for detecting spikes. This paper is structured as follows. In Section II, we describe the experimental dataset and the spike detection framework, with emphasize on STEH for threshold selection. Results for various spike detectors, with and without signal prewhitening and STEH, are presented in Section III and discussed with concluding remarks in Section IV. All analysis was performed using MATLAB (version 2011b) unless otherwise stated.
II. MATERIALS AND METHODS A. Dataset Two datasets, labeled as dataset 1 and dataset 2, were respectively obtained from OSort simulations 2 and 3 with noise standard deviation of 0.15, which are good representations of neural signals in a real environment . The datasets, each 100 s at sampling frequency of 25 kHz, were generated from 150 well-separated spikes. Background noise was created by randomly choosing and scaling the spike waveforms before adding it to noise traces. Dataset 1 has three neurons, each simulated by a renewal Poisson process with a refractory period of 3 ms and a mean firing rate of 5, 7, and 4 Hz, which correspond to an average SNR of 1.4, 1.4, and 2.3 dB. Dataset 2 has five neurons with mean firing rate of 5, 7, 4, 6, and 9 Hz, and corresponding average SNR of 1.4, 1.3, 0.9, 1.6, and 2.6 dB. There are 1568 and 2986 ground truth spikes in datasets 1 and 2, respectively. B. Spike Detection During an extracellular recording, neural signals are often corrupted by two types of noises: neural noise induced by neuronal activity far from the recording electrode, as well as thermal and electrical noise from electronic circuitry in the recording system . In this study, we attempted to suppress the neural noise, which is colored and correlated, through a 4th-order linear predictive coding prewhitening filter that can decorrelate the noise and whiten the signals . The filter coefficients were computed from pure noise segment using autocorrelation method and Levinson-Durbin recursion . Moreover, we bandpass filtered the signals through a 4th-order zero-phase lag Butterworth filter (300‒3000 Hz) to eliminate low-frequency local field potentials and high-frequency thermal and electrical noise. To further accentuate spikes, we pre-emphasized the filtered signals using four different techniques, as elaborated below, along with their respective threshold setting. 1) Absolute value: A spike can have both positive and negative amplitudes, thus applying a threshold to the signal absolute value is equivalent to applying positive and negative thresholds. As suggested in , the threshold was Thr M ,
0 . 6745
where M = 4 is the threshold multiplier, σ is the estimate of noise standard deviation, and xn is the signal at time n. 2) Teager energy operator: A spike is described a shortlived burst with high amplitude and frequency. TEO ,
(3) where denotes the convolution operator, and wn denotes the Hamming window of length = 5 . The threshold,
3) Higher-order differential energy operator: The concept of instantaneous energy can be broadened to higher dimensions by measuring the cross energy between the signal and its higher derivatives . The 3rd- and 4th-order differential energy operators, which correspond to an energy velocity operator , V x n x n x n 1 x n 1 x n 2 x n 1 x n x n 2 x n 1 2 ,
and an energy acceleration operator , (6) A x n x n x n 2 x n 1 x n 3 . were assessed, for the first time as energy-driven spike detectors. Similar to the TEO, which is a special case of 2ndorder differential energy operator, the thresholds for energy velocity and energy acceleration operators were scaled version of the mean of their respective operator. To realize a dynamic threshold setting, we exploited STEH because STEO detector outperforms other amplitudeand energy-driven spike detectors, as rendered in Figure 1. For better understanding of STEH, we examined two choices of computing STEH bin width, namely Freedman-Diaconis choice , 1 3 , w 2 iqr S x n N (7) which depends on the interquartile range of STEO and produces number of bins, b max
x n min S x n w 1 ,
(8) as well as the common square-root choice whose number of bins is the square-root of N data samples, b = sqrt (N). S
In addition, we incorporated histogram equalization  to augment the separability between spike and noise distributions, similar to enhancing image contrast. The equalized STEH at bin index k is cdf k cdf k 1 k H k cdf k
k 1 k 1
where cdf is the cumulative distribution function of STEH. Eventually, an adaptive threshold ThrSTEH for spike detection can be automatically determined by the cutoff value T that maximizes the entropies between spike HS and noise HN distributions. Thr STEH arg max H S T H N T
(2) is sensitive to spikes because it can heighten local peaks in both amplitude and frequency. Since spikes are nonstationary phenomenon, the expectation operator of TEO cannot be substituted by a time-domain averaging but a frequencydomain windowing, which gives rise to STEO ,
S x n ,
was a scaled version of the mean STEO , with N data samples and threshold multiplier M = 8 .
x n x n x n 1 x n 1 ,
S xn xn wn ,
H S T
g 0 P T
H N T
g T 1
p g P T
p g P T
p g P T
analogous to thresholding gray-level image via histogrambased entropy .
3 6 9 12 False alarm (spikes/s)
90 80 70 60 50
90 80 70 60 50
3 6 9 12 False alarm (spikes/s)
True detection rate (%)
True detection rate (%)
Dataset 1 True detection rate (%)
True detection rate (%)
6 8 10 Threshold multiplier
90 80 70 60 50
6 8 10 Threshold multiplier
Figure 2. True detection rate (%) versus threshold multiplier for STEO detector without (grey solid line) and with (black solid line) signal prewhitening. Asterisks indicate performance of STEO with threshold multiplier M = 8, while circles and squares respectively indicate STEH with and without histogram equalization. Shaded markers represent Freedman-Diaconis choice of bin width, while unshaded markers represent square-root choice of bin width.
Figure 1. Receiver operating characteristic curves of four spike detectors without signal prewhitening: absolute value (dash-dot line), energy velocity operator (dash line), energy acceleration operator (dot line), and STEO (grey solid line); and a STEO detector with signal prewhitening (black solid line). Asterisks indicate performance of STEO with threshold multiplier M = 8, while circles and squares respectively indicate STEH with and without histogram equalization. Shaded markers represent Freedman-Diaconis choice of bin width, while unshaded markers represent square-root choice of bin width.
For each detected spike, 64 data samples (2.56 ms) were extracted, with the peak at sample 20. To lessen the effects of sampling jitter and improve spike alignment, we first interpolated the data samples four times using fast Fourier transform method, then re-aligned such that the peak was at sample 93, and finally downsampled to retain the original spike with sampling frequency of 25 kHz.
histogram equalization can raise TDR by 2.84 ± 2.47% and FA by 1.01 ± 0.34 spikes/s. For STEO detector with signal prewhitening, thresholds determined by STEH with histogram equalization lie around knee point of its ROC curve. Under such condition, STEH with Freedman-Diaconis choice yields TDR = 87.88 ± 0.10% and FA = 1.82 ± 0.55 spikes/s, whereas STEH with square-root choice yields TDR = 86.49 ± 0.07% and FA = 1.41 ± 0.36 spikes/s, as detailed in Table I.
C. Performance Measures Receiver operating characteristic (ROC) curves  were constructed to appraise the performance of spike detectors, with and without signal prewhitening and STEH. The ROC curve displays true detection rate (TDR, percentage of correctly detected spikes) on the vertical axis and false alarm (FA, number of falsely detected spikes per second) on the horizontal axis. A spike was deemed correctly detected if the time difference between the spike and the ground truth was within ± 0.4 ms, and there was no double counting.
To further evaluate the competency of automated threshold selection using STEH, we compared its performance and corresponding threshold multiplier value with that of manual selection (M = 8 for STEO detector). STEH with signal prewhitening, Freedman-Diaconis choice or square-root choice of bin width, and histogram equalization achieve relatively lower TDR and FA for dataset 1 but higher TDR and FA for dataset 2. These findings are due to different threshold values determined by STEH, which are equivalent to M = 8.3‒8.6 for dataset 1 and M = 6.9‒7.3 for dataset 2 in (4), as rendered in Figure 2.
Using the datasets 1 and 2 without signal prewhitening, Figure 1 demonstrates that STEO detector is superior to absolute value detector and comparable to energy velocity and energy acceleration operators, which implies STEO as a simple and efficient spike detector. Signal prewhitening can substantially improve detection outcomes. As illustrated by STEO detector in Figure 1, at constant FA ranging 3‒12 spikes/s, there is 10‒20% TDR difference between STEO detector with and without signal prewhitening.
This paper investigates and proposes the use of STEH to automatically select spike detection threshold. Results from stimulated datasets containing spikes of SNR = 0.9‒2.6 dB reveal that (1) prewhitening of neural signals can enhance TDR by 10‒20% at constant FA spanning 3‒12 spikes/s; (2) Freedman-Diaconis choice of STEH bin width yields higher TDR (1.52 ± 1.41%) and FA (0.48 ± 0.25 spikes/s) than square-root choice; and (3) histogram equalization can increase TDR by 2.84 ± 2.47% and FA by 1.01 ± 0.34 spikes/s. Thresholds determined via STEH with signal prewhitening, Freedman-Diaconis choice or square-root choice, and histogram equalization are around knee point of ROC curve.
Figure 1 and Table I present the detection performance of using STEH to select spike detection threshold. Comparing the two choices of computing STEH bin width for both datasets with and without signal prewhitening, we noticed that Freedman-Diaconis choice consistently attains higher TDR (mean ± standard deviation of 1.52 ± 1.41%) and FA (0.48 ± 0.25 spikes/s) than square-root choice. Besides that,
Neural signals are often contaminated by two types of noises during an extracellular recording. Among these noises, neural noise is relatively difficult to handle because it
DETECTION PERFORMANCE OF THRESHOLD SELECTION VIA SMOOTHED TEAGER ENERGY HISTOGRAM UNDER DIFFERENT CONDITIONS Signal Prewhitening Square-Root Freedman-Diaconis HE w/o HE HE w/o HE 86.54 78.38 87.95 83.04 1.15 0.19 1.43 0.46 86.44 84.86 87.81 86.57 1.66 1.02 2.21 1.70
TDR (%) FA (spikes/s) TDR (%) FA (spikes/s)
Dataset 1 Dataset 2
w/o Signal Prewhitening Square-Root Freedman-Diaconis HE w/o HE HE w/o HE 70.03 68.62 70.34 68.81 5.22 4.08 5.77 4.20 76.19 73.68 76.89 75.52 5.72 4.38 6.22 5.26
Amplitude of STEO
TDR refers to true detection rate in percentage (%); FA, false alarm in number of spikes per second (spikes/s); HE, histogram equalization; w/o, without
Figure 3 exemplifies a neural signal acquired from a monkey, in collaboration with Singapore Institute for Clinical Sciences and National University of Singapore.
6 4 2
Figure 3. A preprocessed real neural signal of smoothed Teager energy operator (STEO) detector with signal prewhitening, and threshold values for threshold multiplier M = 8 (dot line) and for automated selection via smoothed Teager energy histogram with Freedman-Diaconis choice of bin width and histogram equalization (solid line).
is colored and correlated . In this study, we adopted a linear predictive coding prewhitening filter to decorrelate or pose statistical independence between noise, enhancing neural signal quality and improving spike detection performance. We also pre-emphasized the filtered neural signals using four different amplitude- and energy-based techniques and noted that STEO can best accentuate spikes relative to background noise, which is attributable to its nonlinear energy-tracking ability using both signal amplitude and frequency information . Consequently, STEO detector delivers superior spike detection performance, which is in agreement with the literature [5,6]. Integrating STEH for autonomous selection of spike detection threshold can further boost the usefulness of STEO detector in a real environment. Thresholds determined by STEH, with histogram equalization that maximizes separation between spike and noise distributions, fall around knee point of the ROC curve of STEO detector with signal prewhitening. STEH bin width estimated by FreedmanDiaconis choice produces higher TDR and FA than squareroot choice because the former choice leads to smaller bin width and more bins, thereby offering more details about the spike and noise distributions for the subsequent histogram equalization and histogram-based entropy thresholding. Owing to the fact that false alarms will be discarded during spike sorting, we propose using STEH with Freedman-Diaconis choice of bin width (TDR = 87.88 ± 0.10% and FA = 1.82 ± 0.55 spikes/s), in addition to signal prewhitening and histogram equalization, for selecting spike detection threshold because its resultant threshold is more permissive, yielding more true detections with limited false alarms. While this study has demonstrated the feasibility of utilizing STEH to automatically select spike detection threshold, the results are based on stimulated datasets. Therefore, we are currently working on real datasets, and
  
       
M.A.L. Nicolelis, Methods for Neural Ensemble Recordings. Boca Raton, FL: CRC Press, 2008, pp. 57–82. M.S. Lewicki, “A review methods for spike sorting: The detection and classification of neural action potentials,” Network: Comput. Neural Syst., vol. 9, pp. R58–R78, 1998. S. Gibson, J.W. Judy, and D. Markovic, “Spike sorting: The first in decoding the brain,” Signal Process. Mag., vol. 29, pp. 124–143, 2012. I. Obeid and P.D. Wolf, “Evaluation of spike-detection algorithms for a brain-machine interface application,” IEEE Trans. Biomed. Eng., vol. 51, pp. 905–911, 2004. S. Gibson, J.W. Judy, and D. Markovic, “Technology-aware algorithm design for neural spike detection, feature extraction, and dimensionality reduction,” IEEE Trans. Neural Syst. Rehabil. Eng., vol. 18, pp. 469–478, 2010. S. Mukhopadhyay and G.C. Ray, “A new interpretation of nonlinear energy operator and its efficacy in spike detection,” IEEE Trans. Biomed. Eng., vol. 45, pp. 180–187, 1998. M. Rizk and P.D. Wolf, “Optimizing the automatic selection of spike detection thresholds using a multiple of the noise level,” Med. Biol. Eng. Comput., vol. 47, pp. 955–966, 2009. H. Semmaoui, J. Drolet, A. Lakhssassi, and M. Sawan, “Setting adaptive spike detection threshold for smoothed TEO based on robust statistics theory,” IEEE Trans. Biomed. Eng., vol. 59, pp. 474–482, 2012. M.B. Malarvili, H. Hassanpour, M. Mesbah, and B. Boashash, “A histogram-based electroencephalogram spike detection,” in Proc. 8th Int. Symp. Signal Process. Appl., 2005, pp. 207–210. U. Rutishauser, E.M. Schuman, and A.N. Mamelak, “Online detection and sorting of extracellularly recorded action potential in human medial temporal lobe recordings, in vivo,” J. Neurosci. Methods, vol. 154, pp. 204–224, pp. 2006. K.G. Oweiss, Statistical Signal Processing For Neuroscience and Neurotechnology. New York: Academic Press, 2010, pp. 15–74. L. Lennart, System Identification: Theory for the User. New Jersey: Prentice Hall, 1999, pp. 197–246. R.Q. Quiroga, Z. Nadasdy, and Y. Ben-Shaul, “Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering,” Neural Comput., vol. 16, pp. 1661–1687, 2004. P. Maragos and A. Potamianos, “Higher order differential energy operators,” IEEE Signal Process. Lett., vol. 2, pp. 152–154, 1995. D. Freedman and P. Diaconis, “On the histogram as a density estimator: L2 theory,” Probab. Theory Relat. Field, vol. 57, pp. 453– 476, 1981. K.Y. Kwon and K. Oweiss, “Wavelet footprints for detection and sorting of extracellular neural action potentials,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2011, pp. 609–612. J.N. Kapur, P.K. Sahoo, and A.K.C. Wong, “A new method for graylevel picture thresholding using entropy of the histrogram,” Comput. Vis. Graph Image Process., vol. 29, pp. 273–285, 1985. J.A. Swets, Signal Detection Theory and ROC Analysis in Psychology and Diagnostics: Collected Papers. New Jersey: Erlbaum, 1996, pp. 7–56.