- Email: [email protected]
Spike sorting: a novel shift and amplitude invariant technique
Neurocomputing 44–46 (2002) 1133 – 1139
www.elsevier.com/locate/neucom
Spike sorting: a novel shift and amplitude invariant technique Karim G. Oweissa; ∗, David J. Andersona; b a Electrical
Engineering and Computer Science Department, University of Michigan, Ann Arbor, MI, USA b Biomedical Engineering Department, University of Michigan, Ann Arbor, MI, USA
Abstract This paper deals with the spike classi.cation problem encountered in multi-unit recordings of neural activity in the brain. We recently developed a new methodology for estimating and classifying multi-units recorded by means of multichannel silicon probes from the observed spike trains (Proceedings of the ICASSP’01, May 2001, pp. 2813–2816; Proceedings of the IEEE 35th Asilomar Conference on Signals, Systems and Computers, Paci.c Grove, CA, November 2001). In this work, we demonstrate the robustness of the technique to single unit spike amplitude variation often encountered in burst activity or long term chronic recordings. In low signal-to-noise ratio scenarios where variability in spike threshold crossings during classical detection is always a problem, we show that the technique is extremely robust to shifts in spike event times. Results showing the e=ciency of the algorithm from simulated and experimental data are presented. c 2002 Elsevier Science B.V. All rights reserved. Keywords: Multichannel recording; Spike sorting; Array processing; Wavelet transform
1. Introduction Extracellular recordings of neural activity in the brain has been the subject of many research eAorts for many decades. The recent advances in computational power, data acquisition systems and multichannel microelectrode systems design, and in particular silicon substrate microprobes [2], increased the feasibility of recording multiple neural cells simultaneously from small populations of neurons. As the number of channels in a recording device increases, the number of units recorded increases in direct proportion. ∗
Corresponding author. E-mail address: [email protected] (K.G. Oweiss).
c 2002 Elsevier Science B.V. All rights reserved. 0925-2312/02/$ - see front matter PII: S 0 9 2 5 - 2 3 1 2 ( 0 2 ) 0 0 4 3 6 - 8
1134
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
Since spike waveforms tend to be correlated, the spike sorting stage following recording and detection is increasingly challenged by a harder spike discrimination task. There are numerous techniques for spike sorting that have been developed in the past. We will not attempt to review them since this is not the objective of this paper, but the interested reader is referred to [6,8]. The most widely known methods rely on principle component analysis and template matching [1]. In the ideal case of a stable recording and reasonable signal-to-noise ratio (SNR), the detection of spikes is achieved using amplitude threshold crossing or energy measures associated with a certain time window extracted surrounding the threshold crossing point. Usually, a .xed number of points before and after the threshold crossing point is set a priori, and all spikes are extracted accordingly. If the detection is always perfect, spikes belonging to the same cell are perfectly aligned to each other and in most cases have very close peak-to-peak amplitude and exact shape. This is not always true in practice for many reasons. First, the additive noise always alters the threshold crossing points among single cell spike waveforms, which implies that aligning is necessary before classi.cation can be done. Second, assuming the alignment was done, the nonstationarity in the noise characteristics and the variation in spike peak-to-peak amplitude over time due to electrode drift make the .rst problem persistent. Third, spike waveforms belonging to the same cell exhibit variations in amplitude and in waveform shape in diAerent scenarios such as burst activity or during long term chronic recordings where animal motion is inevitable [4,11]. All existing spike sorting techniques do not take into account these facts and thus exhibit large error margins with increasing neural activity. In some cases, some techniques rely on semi-automatic clustering permitting the user to judge the clustering performance in a trial to remedy these shortcomings. The transform domain array processing framework that we developed in [9,10], showed very promising results especially in cases where other existing techniques have failed to yield good performance. In this work, we demonstrate with various examples its ability to deal with these shortcomings, namely, single unit spike waveform variations, in a very e=cient way. In what follows, we brieHy describe the main steps of the algorithm and then focus on a single cell’s spike waveform and introduce variations in amplitude and waveform shape as well as shifts in spike event time occurrence. We show that the technique yields a signi.cant robustness against those variations. 2. Multiresolution analysis of signal subspace invariance technique (MASSIT): a brief overview The MASSIT technique is based on performing a multiresolution analysis of the spike waveform by means of the discrete wavelet packet transform (DWPT) [7]. The DWPT is better understood using a dyadic tree concept as shown in Fig. 1, where the root constitute the signal level and nodes represent successive projections onto subspaces spanned by the corresponding wavelet basis functions [11]. The wavelet expansions obtained in the DWPT tree is overcomplete and redundant. Best signal representation can be obtained by pruning the wavelet packet tree according to a best basis selection criterion. This is usually achieved by associating a cost function to each node of the
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
1135
Fig. 1. 3-level Binary Tree illustrating the DWPT decomposition. Nodes are labeled from left to right and from top to bottom starting at 0 (signal level, i.e. observation data matrix X ).
tree and comparing the parent node’s cost to its children’s cost. A node is kept for further splitting if it yields a lower cost than its children [3,5]. We will denote the M × N observation data matrix by X , where M is the number of channels, N is the number of snapshots recorded by the array in time interval T . The L-level DWPT of X , denoted Wx will constitute the set {Wx(1) ; Wx(2) ; : : : ; Wx(L) }, where Wx(i) is a set of matrices that constitute the DWPT of X at level i. In each decomposition level i, Wx(i) comprises the set {Wx(i)0 ; Wx(i)1 ; : : : ; Wx(i)K }, where Wx(i)k is the DWPT at the kth node, and k varies from 0 to K = 2i − 1 diAerent nodes obtained in level i. Using this convention, we note that i = 0 corresponds to the signal level and i equals L corresponds to the coarsest level of decomposition which is determined by a dyadic relation to the data vector length N , i.e., L = log2 N . We perform a singular value decomposition (SVD) of the spatial covariance matrix of Wx(i)k (i = 1; 2; : : : ; L and k = 0; 1; : : : ; 2i − 1) to obtain the ordered singular values 1 ; 2 ; : : : ; M for each node in the tree. We associate a characteristic best tree to each spike by searching for a subtree of the full tree with nodes having the largest subset of maximum 1 ’s. This corresponds to selecting frequency subbands where the spike’s energy is most concentrated by projecting the spike waveform onto subspaces spanned by wavelet basis functions that form wavelet low pass and band pass .lters having the closest match to the spike waveform. We omit here the mathematical details of the algorithm which can be found in [9,10]. 3. Results The MASSIT is implemented using an undecimated version of the DWPT. The reason is two fold: 1. The transform in this case is shift invariant, since all shifts in the input waveform are propagated within the transform coe=cients such that they are intertwined together. 2. The spatial covariance matrix estimate needed for performing the SVD will be based on averaging N outer products in each node at the jth level rather than N=2j in the
1136
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
orthogonal DWPT case. This yields a closer estimate to the true spatial covariance matrix. We present here results from two simulations. The .rst one simulated a single unit burst activity consisting of a spike train with equal inter-event time durations and decreasing spike magnitude and shape change across time. The elements of the spike train were obtained from a single unit template waveform that was successively low pass .ltered using an FIR .lter of diAerent lengths. The purpose is two fold: 1. Introduce a slight variation in the spike shape and amplitude to mimic a real data scenario. 2. Introduce a time shift due to the .lter delay to mimic diAerent event detection times within the detection window. Spike detection was performed using a classical threshold crossing criterion estimated based on the noise variance from data with no signal present. Once an event is detected, a window surrounding the event is extracted such that no other spikes are contained within the window. Fig. 2a shows the simulated 4 channel data set with spatially correlated experimental noise. The spike train had 10 events that were extracted and displayed in Fig. 2b for illustration of waveform shape changes within a single cell burst activity. By carefully examining Fig. 2d, consider nodes (31) and (34) for example across all spikes. The percentage decrease in 1 for node (34) is much bigger than that of node (31), which can be interpreted that node (31) is less sensitive to changes in the spike’s energy than node (34), yet they are both included in the characteristic tree. The same argument can be generalized by comparing the change in 1 ’s for other nodes. Generally speaking, a node is more sensitive to changes in spike’s energy if the corresponding wavelet basis has a close match to the spike waveform shape. We carried out a second experiment in which 2 cells having exactly the same spike shape but with diAerent spatial amplitude distribution across channels were simulated. Fig. 3 shows a 4 channel data where each cell .red twice in a 100 ms time interval. Spikes 1 and 3 belong to cell 1 while spikes 2 and 4 belong to cell 2. Both cells were scaled according to the weight vectors shown in Fig. 3b. Obviously, maximum ambiguity occurs on channel 2 where the cells were almost equally weighted. No sorting technique would be able to resolve the ambiguity based on single channel data, especially channel 2. Moreover, applying the .rst stage of MASSIT, the eigenvalue pro.le across nodes directly mapping the spectral energy distribution would reveal exactly the same characteristic tree for both cells yielding a sorting error. However, MASSIT’s second stage oAers extra information that can be used to resolve this ambiguity. Since the largest singular value is associated with an eigenvector that reveals information on the spatial distribution of the signal’s energy, then simply examining the coe=cients of this eigenvector would yield the spatial distribution as a basis for discrimination. Fig. 3c and d and Fig. 3e illustrate this idea. The principal eigenvector corresponding to 1 across nodes is plotted for spikes 1 and 2 as detected on channels 4 and 1, respectively.
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
1137
Fig. 2. (a) Noisy simulated burst activity across 4 channels, and (b) extracted spikes after detection from channel 1 of the 4 channel data in (a) superimposed together for waveform comparison. (c) Characteristic best basis tree obtained from MASSIT, (d) pro.le of eigen values across nodes of the tree in (c), and (e) signal subspace (spike footprints) across channels. Note the decrease in eigenvalue in nodes at coarse levels in the tree (nodes 16 and above) while the signal subspace in (e) remains almost invariant.
1138
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
Fig. 3. (a) 4-channel data set with 2 cells each .ring twice with exact similar spike shapes but diAerent spatial amplitude distribution according to (b). Note that maximum ambiguity in sorting occurs on channel 2 where the spikes were closely scaled. (c) Characteristic tree for all 4 spikes detected in (a). (d) Eigenvalue pro.le across nodes. (e) Estimated Signal subspace.
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
1139
4. Conclusion Promising results have been presented to show the capability of the MASSIT algorithm to deal with variations of single cell discharge patterns, in particular, amplitude and shift variations in burst activity and long term chronic recordings. The second advantage superceeds existing techniques in that it is based on performing sensor array processing techniques using signal and noise subspace analysis methods to simultaneously characterize spatial, temporal and spectral characteristics of multi unit neural activity in a very e=cient, cost inexpensive way. Acknowledgements This work was supported by NIH under contract number P 41 RR09754. References [1] M. Abeles, M.H. Goldstein Jr., Multispike train analysis, Proc. IEEE 65 (1977) 762–773. [2] D.J. Anderson, K. Naja., S.J. Tanghe, D.A. Evans, K.L. Levy, J.F. Hetke, X. Xue, J.J. Zappia, K.D. Wise, Batch-fabricated thin-.lm electrodes for stimulation of the central auditory system, IEEE Trans. Bio-Med. Electron. 36 (7) (1989) 693–704. [3] R. Coifman, M. Wickerhauser, Entropy based algorithms for best basis selection, IEEE Trans. Inform. Theory 38 (1992) 713–718. [4] K. Harris, et al., Accuracy of Tetrode spike separation as determined by simultaneous intracellular and extracellular measurements, J. Neurophysiol. 84 (2000) 401–414. [5] H. Krim, et al., Best basis algorithm for signal enhancement, in: Proceedings of the IEEE ICASSP’95, May 1995, Detroit, MI, pp. 1561–1564. [6] M.S. Lewicki, A review of methods for spike sorting: the detection and classi.cation of neuronal action potentials, Network: Comput. Neural Syst. 9 (1998) R53–R78. [7] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, New York, 2nd Edition, 1999, 413p. [8] M. Nicolelis (Ed.), Methods for Neural Ensemble Recordings, CRC Press, Boca Raton, FL, 1999 (Chapter 4). [9] K.G. Oweiss, D.J. Anderson, A new technique for blind source separation using subband subspace analysis in correlated multichannel signal environments, in: Proceedings of the ICASSP’01, May 2001, Salt Lake City, UT, pp. 2813–2816. [10] K.G. Oweiss, D.J. Anderson, MASSIT-multiresolution analysis of signal subspace invariance techniques: a novel technique for blind source separation, in: Proceedings of the IEEE 35th Asilomar Conference on Signals, Systems and Computers, Paci.c Grove, CA, November 2001, pp. 819–823. [11] M. Quirk, M. Wilson, Interaction between spike waveform classi.cation and temporal sequence detection, J. Neurosci. Meth. 94 (1999) 41–52.
www.elsevier.com/locate/neucom
Spike sorting: a novel shift and amplitude invariant technique Karim G. Oweissa; ∗, David J. Andersona; b a Electrical
Engineering and Computer Science Department, University of Michigan, Ann Arbor, MI, USA b Biomedical Engineering Department, University of Michigan, Ann Arbor, MI, USA
Abstract This paper deals with the spike classi.cation problem encountered in multi-unit recordings of neural activity in the brain. We recently developed a new methodology for estimating and classifying multi-units recorded by means of multichannel silicon probes from the observed spike trains (Proceedings of the ICASSP’01, May 2001, pp. 2813–2816; Proceedings of the IEEE 35th Asilomar Conference on Signals, Systems and Computers, Paci.c Grove, CA, November 2001). In this work, we demonstrate the robustness of the technique to single unit spike amplitude variation often encountered in burst activity or long term chronic recordings. In low signal-to-noise ratio scenarios where variability in spike threshold crossings during classical detection is always a problem, we show that the technique is extremely robust to shifts in spike event times. Results showing the e=ciency of the algorithm from simulated and experimental data are presented. c 2002 Elsevier Science B.V. All rights reserved. Keywords: Multichannel recording; Spike sorting; Array processing; Wavelet transform
1. Introduction Extracellular recordings of neural activity in the brain has been the subject of many research eAorts for many decades. The recent advances in computational power, data acquisition systems and multichannel microelectrode systems design, and in particular silicon substrate microprobes [2], increased the feasibility of recording multiple neural cells simultaneously from small populations of neurons. As the number of channels in a recording device increases, the number of units recorded increases in direct proportion. ∗
Corresponding author. E-mail address: [email protected] (K.G. Oweiss).
c 2002 Elsevier Science B.V. All rights reserved. 0925-2312/02/$ - see front matter PII: S 0 9 2 5 - 2 3 1 2 ( 0 2 ) 0 0 4 3 6 - 8
1134
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
Since spike waveforms tend to be correlated, the spike sorting stage following recording and detection is increasingly challenged by a harder spike discrimination task. There are numerous techniques for spike sorting that have been developed in the past. We will not attempt to review them since this is not the objective of this paper, but the interested reader is referred to [6,8]. The most widely known methods rely on principle component analysis and template matching [1]. In the ideal case of a stable recording and reasonable signal-to-noise ratio (SNR), the detection of spikes is achieved using amplitude threshold crossing or energy measures associated with a certain time window extracted surrounding the threshold crossing point. Usually, a .xed number of points before and after the threshold crossing point is set a priori, and all spikes are extracted accordingly. If the detection is always perfect, spikes belonging to the same cell are perfectly aligned to each other and in most cases have very close peak-to-peak amplitude and exact shape. This is not always true in practice for many reasons. First, the additive noise always alters the threshold crossing points among single cell spike waveforms, which implies that aligning is necessary before classi.cation can be done. Second, assuming the alignment was done, the nonstationarity in the noise characteristics and the variation in spike peak-to-peak amplitude over time due to electrode drift make the .rst problem persistent. Third, spike waveforms belonging to the same cell exhibit variations in amplitude and in waveform shape in diAerent scenarios such as burst activity or during long term chronic recordings where animal motion is inevitable [4,11]. All existing spike sorting techniques do not take into account these facts and thus exhibit large error margins with increasing neural activity. In some cases, some techniques rely on semi-automatic clustering permitting the user to judge the clustering performance in a trial to remedy these shortcomings. The transform domain array processing framework that we developed in [9,10], showed very promising results especially in cases where other existing techniques have failed to yield good performance. In this work, we demonstrate with various examples its ability to deal with these shortcomings, namely, single unit spike waveform variations, in a very e=cient way. In what follows, we brieHy describe the main steps of the algorithm and then focus on a single cell’s spike waveform and introduce variations in amplitude and waveform shape as well as shifts in spike event time occurrence. We show that the technique yields a signi.cant robustness against those variations. 2. Multiresolution analysis of signal subspace invariance technique (MASSIT): a brief overview The MASSIT technique is based on performing a multiresolution analysis of the spike waveform by means of the discrete wavelet packet transform (DWPT) [7]. The DWPT is better understood using a dyadic tree concept as shown in Fig. 1, where the root constitute the signal level and nodes represent successive projections onto subspaces spanned by the corresponding wavelet basis functions [11]. The wavelet expansions obtained in the DWPT tree is overcomplete and redundant. Best signal representation can be obtained by pruning the wavelet packet tree according to a best basis selection criterion. This is usually achieved by associating a cost function to each node of the
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
1135
Fig. 1. 3-level Binary Tree illustrating the DWPT decomposition. Nodes are labeled from left to right and from top to bottom starting at 0 (signal level, i.e. observation data matrix X ).
tree and comparing the parent node’s cost to its children’s cost. A node is kept for further splitting if it yields a lower cost than its children [3,5]. We will denote the M × N observation data matrix by X , where M is the number of channels, N is the number of snapshots recorded by the array in time interval T . The L-level DWPT of X , denoted Wx will constitute the set {Wx(1) ; Wx(2) ; : : : ; Wx(L) }, where Wx(i) is a set of matrices that constitute the DWPT of X at level i. In each decomposition level i, Wx(i) comprises the set {Wx(i)0 ; Wx(i)1 ; : : : ; Wx(i)K }, where Wx(i)k is the DWPT at the kth node, and k varies from 0 to K = 2i − 1 diAerent nodes obtained in level i. Using this convention, we note that i = 0 corresponds to the signal level and i equals L corresponds to the coarsest level of decomposition which is determined by a dyadic relation to the data vector length N , i.e., L = log2 N . We perform a singular value decomposition (SVD) of the spatial covariance matrix of Wx(i)k (i = 1; 2; : : : ; L and k = 0; 1; : : : ; 2i − 1) to obtain the ordered singular values 1 ; 2 ; : : : ; M for each node in the tree. We associate a characteristic best tree to each spike by searching for a subtree of the full tree with nodes having the largest subset of maximum 1 ’s. This corresponds to selecting frequency subbands where the spike’s energy is most concentrated by projecting the spike waveform onto subspaces spanned by wavelet basis functions that form wavelet low pass and band pass .lters having the closest match to the spike waveform. We omit here the mathematical details of the algorithm which can be found in [9,10]. 3. Results The MASSIT is implemented using an undecimated version of the DWPT. The reason is two fold: 1. The transform in this case is shift invariant, since all shifts in the input waveform are propagated within the transform coe=cients such that they are intertwined together. 2. The spatial covariance matrix estimate needed for performing the SVD will be based on averaging N outer products in each node at the jth level rather than N=2j in the
1136
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
orthogonal DWPT case. This yields a closer estimate to the true spatial covariance matrix. We present here results from two simulations. The .rst one simulated a single unit burst activity consisting of a spike train with equal inter-event time durations and decreasing spike magnitude and shape change across time. The elements of the spike train were obtained from a single unit template waveform that was successively low pass .ltered using an FIR .lter of diAerent lengths. The purpose is two fold: 1. Introduce a slight variation in the spike shape and amplitude to mimic a real data scenario. 2. Introduce a time shift due to the .lter delay to mimic diAerent event detection times within the detection window. Spike detection was performed using a classical threshold crossing criterion estimated based on the noise variance from data with no signal present. Once an event is detected, a window surrounding the event is extracted such that no other spikes are contained within the window. Fig. 2a shows the simulated 4 channel data set with spatially correlated experimental noise. The spike train had 10 events that were extracted and displayed in Fig. 2b for illustration of waveform shape changes within a single cell burst activity. By carefully examining Fig. 2d, consider nodes (31) and (34) for example across all spikes. The percentage decrease in 1 for node (34) is much bigger than that of node (31), which can be interpreted that node (31) is less sensitive to changes in the spike’s energy than node (34), yet they are both included in the characteristic tree. The same argument can be generalized by comparing the change in 1 ’s for other nodes. Generally speaking, a node is more sensitive to changes in spike’s energy if the corresponding wavelet basis has a close match to the spike waveform shape. We carried out a second experiment in which 2 cells having exactly the same spike shape but with diAerent spatial amplitude distribution across channels were simulated. Fig. 3 shows a 4 channel data where each cell .red twice in a 100 ms time interval. Spikes 1 and 3 belong to cell 1 while spikes 2 and 4 belong to cell 2. Both cells were scaled according to the weight vectors shown in Fig. 3b. Obviously, maximum ambiguity occurs on channel 2 where the cells were almost equally weighted. No sorting technique would be able to resolve the ambiguity based on single channel data, especially channel 2. Moreover, applying the .rst stage of MASSIT, the eigenvalue pro.le across nodes directly mapping the spectral energy distribution would reveal exactly the same characteristic tree for both cells yielding a sorting error. However, MASSIT’s second stage oAers extra information that can be used to resolve this ambiguity. Since the largest singular value is associated with an eigenvector that reveals information on the spatial distribution of the signal’s energy, then simply examining the coe=cients of this eigenvector would yield the spatial distribution as a basis for discrimination. Fig. 3c and d and Fig. 3e illustrate this idea. The principal eigenvector corresponding to 1 across nodes is plotted for spikes 1 and 2 as detected on channels 4 and 1, respectively.
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
1137
Fig. 2. (a) Noisy simulated burst activity across 4 channels, and (b) extracted spikes after detection from channel 1 of the 4 channel data in (a) superimposed together for waveform comparison. (c) Characteristic best basis tree obtained from MASSIT, (d) pro.le of eigen values across nodes of the tree in (c), and (e) signal subspace (spike footprints) across channels. Note the decrease in eigenvalue in nodes at coarse levels in the tree (nodes 16 and above) while the signal subspace in (e) remains almost invariant.
1138
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
Fig. 3. (a) 4-channel data set with 2 cells each .ring twice with exact similar spike shapes but diAerent spatial amplitude distribution according to (b). Note that maximum ambiguity in sorting occurs on channel 2 where the spikes were closely scaled. (c) Characteristic tree for all 4 spikes detected in (a). (d) Eigenvalue pro.le across nodes. (e) Estimated Signal subspace.
K.G. Oweiss, D.J. Anderson / Neurocomputing 44–46 (2002) 1133 – 1139
1139
4. Conclusion Promising results have been presented to show the capability of the MASSIT algorithm to deal with variations of single cell discharge patterns, in particular, amplitude and shift variations in burst activity and long term chronic recordings. The second advantage superceeds existing techniques in that it is based on performing sensor array processing techniques using signal and noise subspace analysis methods to simultaneously characterize spatial, temporal and spectral characteristics of multi unit neural activity in a very e=cient, cost inexpensive way. Acknowledgements This work was supported by NIH under contract number P 41 RR09754. References [1] M. Abeles, M.H. Goldstein Jr., Multispike train analysis, Proc. IEEE 65 (1977) 762–773. [2] D.J. Anderson, K. Naja., S.J. Tanghe, D.A. Evans, K.L. Levy, J.F. Hetke, X. Xue, J.J. Zappia, K.D. Wise, Batch-fabricated thin-.lm electrodes for stimulation of the central auditory system, IEEE Trans. Bio-Med. Electron. 36 (7) (1989) 693–704. [3] R. Coifman, M. Wickerhauser, Entropy based algorithms for best basis selection, IEEE Trans. Inform. Theory 38 (1992) 713–718. [4] K. Harris, et al., Accuracy of Tetrode spike separation as determined by simultaneous intracellular and extracellular measurements, J. Neurophysiol. 84 (2000) 401–414. [5] H. Krim, et al., Best basis algorithm for signal enhancement, in: Proceedings of the IEEE ICASSP’95, May 1995, Detroit, MI, pp. 1561–1564. [6] M.S. Lewicki, A review of methods for spike sorting: the detection and classi.cation of neuronal action potentials, Network: Comput. Neural Syst. 9 (1998) R53–R78. [7] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, New York, 2nd Edition, 1999, 413p. [8] M. Nicolelis (Ed.), Methods for Neural Ensemble Recordings, CRC Press, Boca Raton, FL, 1999 (Chapter 4). [9] K.G. Oweiss, D.J. Anderson, A new technique for blind source separation using subband subspace analysis in correlated multichannel signal environments, in: Proceedings of the ICASSP’01, May 2001, Salt Lake City, UT, pp. 2813–2816. [10] K.G. Oweiss, D.J. Anderson, MASSIT-multiresolution analysis of signal subspace invariance techniques: a novel technique for blind source separation, in: Proceedings of the IEEE 35th Asilomar Conference on Signals, Systems and Computers, Paci.c Grove, CA, November 2001, pp. 819–823. [11] M. Quirk, M. Wilson, Interaction between spike waveform classi.cation and temporal sequence detection, J. Neurosci. Meth. 94 (1999) 41–52.