[1] 15w A method for spike sorting and detection based on wavelet packets and Shannon's mutual information

[1] 174w Studying the dynamics of neural activity via electrical recording, relies on the ability to detect and sort neural spikes recorded from a number of neurons by the same electrode. We suggest the wavelet packets decomposition (WPD) as a tool to analyze neural spikes and extract their main features. The unique quality of the wavelet packets-adaptive coverage of both time and frequency domains using a set of localized packets, facilitate the task. The best basis algorithm utilizing the Shannon's information cost function and local discriminant basis (LDB) using mutual information are employed to select a few packets that are sufficient for both detection and sorting of spikes. The efficiency of the method is demonstrated on data recorded from in vitro 2D neural networks, placed on electrodes that read data from as many as five neurons. Comparison between our method and the widely used principal components method and a sorting technique based on the ordinary wavelet transform (WT) shows that our method is more efficient both in separating spikes from noise and in resolving overlapping spikes.

[1] 67w Current electronics enables the recording of in vitro neural networks cell cultures, hippocampal slices, and in vivo tissues by simultaneous measurement of the activity of many neurons using multichannel electrode arrays (Jimbo et al., 1993;Stenger and Mckenna, 1994;Camepari et al., 1997;Fromherz et al., 1997;Egert et al., 1998;Maher et al., 1999). It is expected to improve our understanding of how real neural networks code, learn and store information.

[2] 70w Extracting useful information from these measurements relies on the ability to detect and sort the recorded neural spikes (Lewicki, 1998). Any method for detecting and sorting of spikes is comprised of two stages: (1) extracting the important features of the spikes and thus lowering the dimension of the parametric set representing the spikes. (2) Clustering of the parametric sets into groups, where each cluster relates to a different spike type.

[3] 72w The clustering problem is widely known and there are many methods to solve it. The performance of any clustering method relies on the original differences among the groups in the data. Hence, the efficiency of the feature extraction method in signifying those differences is crucial for the success of clustering. In this paper, we present a novel feature extraction method that efficiently lowers the dimension while signifying the differences among spike types.

[4] 55w Many powerful clustering methods have been introduced in spike sorting techniques (Lewicki, 1994(Lewicki, , 1998;;Rinberg et al., 1998;Zouridakis and Tam, 2000). However, we will show that the superior performance of our feature extraction technique enables us to successfully cluster the data using merely the k -means technique, which is one of the simplest methods known.

[5] 103w A new method to extract the features of recorded neural spikes is the wavelet transform (WT) (Unser and Aldroubi, 1996;Zouridakis and Tam, 1997;Yang and Shamma, 1998;Letelier and Wever, 2000). The neural spikes, being short voltage pulses, are localized func-tions whose natural analysis is via the WT (Unser and Aldroubi, 1996;Yang and Shamma, 1998), which is the Fourier transform analog for localized functions. Recent publications have shown good performance of the WT in detection and sorting of artificial data (Zouridakis and Tam, 1997;Letelier and Wever, 2000). Yet, the WT does not perform well in separating spikes that overlap or are even close in time.

[6] 162w Here we present a more efficient method based on WPD (Hulata et al., 2000), which is a generalization of the WT. The wavelet packets (Coifman et al., 1992;Daubechies, 1992;Mallat, 1998) are a large family of localized functions, where sets of packets can be adaptively selected to span signals of various characteristics. Algorithms such as best basis and local discriminant basis (LDB) (Coifman et al., 1992;Coifman and Wickerhauser, 1992;Saito, 1994;Mallat, 1998) are used to select sets of packets that enhance the features of the signal using a small number of significant coefficients. In the case of neural spikes, these packets form a compact representation, with minimal coverage of nearby spikes. This gives rise to detection and sorting of spikes even in cases of superimposed (overlapping) spikes. We present the efficiency of our method utilizing the selection of only nine packets. These packets were tested for many different spike waveforms on different electrodes in different experiments, and appear universal for recordings of this kind.

[7] 16w All numerical procedures were coded using MATLAB (www.mathworks.com) and the WAVELAB toolbox (Donoho et al., 2000).

[8] 114w The best basis algorithm is a method to select a basis that spans the signal with the smallest number of significant coefficients. The idea is to find a way to favor either a block with parameters [i, j] or the corresponding next level blocks labeled with [i'/1, 2j] and [i'/1, 2j'/1] (as defined in Fig. 7a). To do so, we choose a cost function I, such that if I i ,j ]/I i'1,2j '/ I i'1 , 2j'1 , we favor block [i, j ], and vice versa. A cost function commonly used for this task, is the Shannon Information Function, which corresponds to minus the Entropy function (Coifman and Wickerhauser, 1992;Saito, 1994;Mallat, 1998):

[9] 13w is the energy contribution of f(t) in packet c i ,j ,k .

[10] 113w Once we have calculated I for all the blocks in the decomposition, we can perform an overall maximization of I , and the resulting set of blocks is the basis of maximum information, and thus defined as the best basis. The maximization is performed as follows: starting at the bottom level i 0/L (the lowest row in Fig. 7a), we compare every pair of blocks with their corresponding block at level i0/L(/1 and choose according to the maximum value of I. After choosing the blocks from these two levels, we go up one level and compare the blocks of levels i0/L(/1 and i0/L(/2. We repeat this procedure up to the top level i0/1.

[11] 41w There is an intuitive relation between the entropy cost function and the physical entropy. Maximum entropy corresponds to an homogeneous distribution over micro-states. Thus, minimum entropy (i.e. maximum information) should correspond to concentration of energy over a minimum of basis elements.

[12] 52w The LDB method (Saito, 1994) is a generalization of the best basis algorithm. It is used when the main point is to enhance the differences between two signals. In resemblance to the best basis, the LDB is a basis that uses the smallest number of significant basis elements to separate two signals.

[13] 49w The scheme of the LDB is as follows: given two signals, f (1) (t) and f (2) (t ), we build for each signal the block diagram (as described before). Next, we choose the cost function M to be the Shannon's mutual information (minus the mutual entropy) (Saito, 1994):

[14] 25w Next, we perform an overall maximization of the mutual information M in a way similar to the maximization of I in the best basis algorithm.

[15] 41w In some cases, to reduce the influence of fluctuations, M is calculated not on the energy distribution of a signal (q i ,j ,k ), but on an averaged measurement. We take N 1 measurements of signal type 1 {f n

[16] 9w , and calculate the energy distribution over the packets:

[17] 33w In this case, the discriminant measure between the two signal types, is calculated on G i ,j ,k 1 and G i ,j ,k 2 , and we continue in the same manner.

[1] 80w The nature of the task is presented in Fig. 1. Fig. 1a presents a recording of a typical neural spike. Fig. 1b presents an example of recorded noise. It was recorded due to the strong peak amplitude. However, it is not a neural spike as the width of the event is less than 1 ms, and it lacks a refractory period. Fig. 1c presents a recording of a number of spikes from five different neurons that fired close in time.

[2] 64w The data acquisition is based on a hardware threshold detector that is adjusted to detect voltage amplitudes, which exceed the background noise by 3 folds (the background noise is estimated during non-firing periods). Whenever an event of amplitude above the threshold is detected, a window is recorded with the event at its center. In our case the length of the window is 20 ms.

[3] 33w Summarizing the main obstacles, a good spike sorting and detection method must carry out the following: 1) detection of spikes in the presence of noise, and distinction between neural spikes and noise events.

[4] 52w 2) Sorting among different neurons read by the same electrode utilizing the differences in the spike waveforms. 3) Detection of close fired spikes, and identification of the origin in presence of superposition. Note that when superposition occurs (as in Fig. 1c), the spike's waveform is not complete and that complicates the identification.

[5] 12w The Fourier transform of a signal f(t ) is generally defined by:

[6] 68w When f(t) is measured over a finite time interval T , and sampled every time interval t , the Fourier transform can also be viewed as a set of N filters (where N 0/T /t). The filters cover the frequency region [(/p/t , p/t], each with a bandwidth of Dv 0/2p/T . The energy distribution of f (t ) over the frequency interval is j/ f/(v )j, and

[7] 6w does not contain any temporal information.

[8] 69w The WT was devised to maintain both the temporal and frequency properties of a signal. This analysis is based on a complete set of localized functions (named wavelets) that span both the time and frequency domains. The WT of f(t) is thus a function of two parameters a and b , where a represents a frequency scale, and b indicates the time location of the wavelet (Daubechies, 1992;Mallat, 1998):

[9] 77w c (t) is called the mother wavelet and c a ,b (t) is a scaled and translated wavelet. Note that c a ,b (t) replaces the e (ivt of Eq. (1). However, unlike the Fourier transform, there are numerous mother wavelets, and c (t) is chosen according to the problem. For example, Fig. 2 presents the mother wavelet we use, the 3 rd order Coiflet mother wavelet, and a number of scaled wavelets. The inverse WT is:

[10] 120w For a signal f(t ) measured over time T and sampled at a time scale of t , the WT can also be viewed as a discrete set of filters. In contrast to the Fourier transform, the set of filters of the WT are not equal in bandwidth. The filters vary between the narrowest filter Dv min 0/p/T to the widest filter Dv max 0/p/2t, which is half the frequency domain. The width of the filter is determined by the value of the scale parameter a : Dv 0/p/at. The discrete values of a are of the type a0/2 i (where i 0/1, . . ., L and L 0/log 2 (N )). Fig. 3 presents this set of filters.

[11] 39w The filters are organized such that the widest filter covers the higher part of the frequency domain, and as i increases, both the bandwidth is narrower and the filter location is lower. The filters also divide the time domain.

[12] 104w The uncertainty principal, Dv Dt ]/p, imposes that as the filter is narrower in bandwidth, the resolution in time is lowered. As a consequence, when the bandwidth is Dv min , Dt0/T and we lose all temporal information. As Dv increases, time resolution is gained. The concept of both time and frequency division is presented in Fig. 6b. In this figure, a neural spike is presented. It's Fourier transform is plotted in vertical, and the WT is shown as a 2D surface, which represents the time Á/frequency space. The surface is covered by rectangles, where the two dimensions represent the frequency and temporal resolutions.

[13] 60w Due to the localization property of the wavelets, the distribution of the signal's energy jW (a, b )j 2 is better concentrated than in the Fourier transform. A comparison between the Fourier and WTs of a neural spike is presented in Fig. 4. In contrast to the WT, the Fourier transform of the signal does not present any profound features.

[14] 109w The usage of the WT for analyzing neural spikes is not optimal. For large values of a (i.e. large scales), the time resolution of the wavelets are wider than the width of the spikes. Therefore, when recording overlapping spikes, the wavelets fail to separate them. Fig. 5 presents a window with three recorded spikes, and the reconstruction of the event using coefficients of the WT. As the wavelets are too wide, they simultaneously cover two spikes. Another example of the inadequate resolution is seen in Fig. 6a. The rectangles that bear the energy of the spike extend beyond the length of the spike both in time and frequency domains.

[15] 77w The timeÁ/frequency space can be divided using a more adaptive set of functions that better localize the spikes. That is implemented via the WPD. Fig. 6b presents the division of the time Á/frequency space according to the WPD. Note that the rectangles localize the spike better than in Fig. 6a. Fig. 5 presents the reconstruction using WPD of a window with multiple spikes. Note that the WPD is able to separate the spikes, not like the WT.

[16] 205w The WPD generalizes the filtering concept of the WT. In the WT, only one filter is defined for each scale a (with bandwidth Dv0/p/at), thus forming a set of filters each half in bandwidth than it's former. In contrast, the key characteristic of the WPD is that each bandwidth is used for forming a complete coverage of the frequency domain, thus creating multiple divisions of the domain. Fig. 7a illustrates this division. The width of the filters is related to the parameter i exactly as in the WT. Note that for each value of i , a series of filters with an equal bandwidth (numbered using the parameter j) cover the entire frequency domain. Thus, the wavelet packets constitute a family of functions characterized by three parameters: v i,j,k (t ) where i 0/1, . . ., log 2 (N ) is the level of decomposition, j0/0, . . ., 2 i (/1 designates the region that is covered on the frequency domain, and k indicates the temporal location of the packet (equivalent to b in the WT). Any selection of i and j indicates a frequency region (also called block) of width Dv 0/p/2 i t covering the range [j Dv , (j'/1)Dv ].

[17] 244w As seen in Fig. 7a, this form of filtering creates an over-complete division of the time Á/frequency space, as every region in both frequency and time domain is covered by more than one block. Any set of blocks that covers the frequency domain without overlapping, also covers the time domain (Fig. 7b presents an example of a possible set). A non-overlapping set of blocks is commonly referred to as a basis, as it spans the time Á/ frequency space. Hence, when analyzing a signal f(t), we would like to find a basis that optimally localizes the features of f(t). Note that in the Fourier domain (b) the energy is disturbed over a wide bandwidth, with no favorable frequencies. In the WT representation (c), the coefficients are lined according to their frequency range, from low to high frequencies. Note that the energy of the spike is concentrated using very few wavelets, whereas the noise is disturbed almost evenly over all the other wavelets. Reconstruction of the spike in the middle using the nine most energetic wavelets. (c) Reconstruction of the spike using nine selected wavelet packets. This quality of reconstruction and the seclusion of the spike demonstrate the efficiency of using optimized wavelet packets. Note that WT did not separate the spike from its consecutive since high scale wavelets are wide and cover more than a single spike. However, the width of the optimized selected packets enable the successful reconstruction of the analyzed spike alone.

[18] 103w Whenever attempting to detect neural spikes merely with a threshold detector, there is a trade off between the level of the threshold and acceptance of noise events. If the threshold is set at a high level, very few noise events will be detected. However, low amplitude spikes will be a priori excluded. This can be mended by lowering the threshold to include the low amplitude spikes accompanied by noise events, and using a well performing sorting method to disqualify the noise. In this sense, a good sorting method can be adopted also to detect spikes whose amplitudes are comparable with the noise level.

[19] 79w Moreover, in cases of rapid firing of spikes (such as in Fig. 1c), a simple real-time threshold detector is not capable of detecting all the peak amplitudes due to the frequent firing. In order to analyze all the spikes that appear in the recording, and not just those that appear in the center of a digitized window, we have developed a method, using wavelet packets, to mark all the peak amplitudes in a window. The procedure is as follows.

[20] 147w We decompose the 20 ms digitized window using all the shifted versions of the wavelet packet c 2,0 (which by applying the best basis technique was found to be the most characteristic packet for detecting a spike, regardless of its origin). On this decomposed signal, we search for all the samples that are local minima, that isextremum among its 12 nearest neighbors (0.5 ms to either side). Each of these peak amplitudes is suspected to be a spike, and a window of 128 samples (Â/10 ms) is isolated with the peak at its center for further WPD analysis and sorting. Fig. 9 presents an example of this procedure. This procedure is motivated by the shiftinvariant WPD (for further mathematical description see Cohen et al., 1997;Zouridakis and Tam, 1997). Note that by doing so, we confine our method to detect spikes that are at least 0.5 ms apart.

[21] 49w To demonstrate the efficiency of our method, we worked on a recording from an electrode that measured neural spikes originating from five different neurons (as presented in Fig. 1c). We manually selected 100 events for each of the neurons and of noise. Fig. 10 presents the different spikes' groups.

[22] 146w We performed the WPD on a window of 128 samples around each event. Since LDB is designed to find the separating basis between two types of signals, we generalized the LDB technique in order to differentiate between four signal types. We apply the LDB procedure on each pair of groups and retrieved a basis separating best among the pair (ten bases for all the different couples). The union of these bases is an over-redundant discriminant basis between all the groups. To extract one non-redundant basis, we first checked which packets appeared in all the bases. We found that the union of those bases was only slightly over-complete, as most of the bases contained similar packets. For the rest of the packets, we selected the level of decomposition that was most common among the ten bases. The resulting basis is a discriminating basis for all six groups.

[23] 60w Among the 128 basis elements, the nine most discriminating packets were selected both for detection and sorting. Those packets are c 4,0,6 , c 4,0,7 , c 4,1,4 , c 3,1,8 , c 3,1,9 , c 2,1,17 , c 2,1,18 , c 3,4,7 and c 5,22,0 . Fig. 8 presents the temporal shape and frequency coverage of the selected packets.

[24] 85w In Fig. 8, we present the decomposition of a typical spike. The spike is decomposed into the nine packets, and then reconstructed. Note that the reconstruction of the spike using the nine packets reveals the unique features of the spike's shape and covers most of the spike's energy (Â/93% of the window's energy). It is possible to view the union of the selected packets also as a highly optimized filter for neural spikes, that preserves its shape and filters out most of the background noise.

[25] 62w Figs. 11 and 12 presents the process of detecting multiple spikes in the same window, where the spikes are from different neurons and also overlap. Each of the detected peaks is independently decomposed using the same nine packets, and is then reconstructed. Note that the packets succeed in preserving three different spike shapes in this window even though two of them overlap.

[26] 218w In order to visualize the separation ability of the selected packets, we embed the coefficients of those packets in a 9 dimensional space. Separation means that spikes from each type are clustered together in distance from clusters of other spike types. In Fig. 13 we present a matrix plot of four of the packets' coefficients. Projections of the clusters on the different planes are plotted on each sub-figure. It is clear that the group of noise events is completely isolated from the five spike 11. WPD in a window with a number of spikes. The top level presents an event with a number of spike events. In the second level, we present a copy of the 128 samples with an event at the center. In the third level, we present the reconstruction of the spike using the nine packets. Note that we succeeded in detecting two different but overlapping spikes (marked by 1 and 3), and sorted them both from the spike marked by 2 and 4. The later two spikes are also very close to one another, hence the shape of spike 4 is a-priori not easy to classify. However, even though the overlapping of the spikes damages the typical shape of the spike, the packets capture the key features of all the spikes in the event.

[27] 52w Fig. 12. Another example of a window with a number of spikes. Note that even though spikes 1 and 2 are overlapping, we succeeded in detecting both spikes. Moreover, the shape of the reconstructed spikes is almost not affected by the overlapping. types, while the different spike types are also well separable.

[28] 227w As the different spike types and noise are clustered in the 9 dimensional space, we can automatically separate the clusters by using a classifier such as the k -means classifier (Bow, 1991;Dellaert, 1999). The k -means algorithm finds the N k most distinct clusters in a multi-dimensional space (where N k is given as an input to the algorithm). We use the Euclidean measure to calculate the distances between vectors in the 9 dimensional space. The k -means is not designed to find the optimal number of clusters in the data. In our case, N k is influenced by the number of spike types in the recording and thus may vary according to the number of neurons measured by the electrode. As a consequence, N k is determined manually for each tested electrode at the beginning of the process. This is done by applying the kmeans classifier on a short segment of the recording (first 2000 events) with different selections on N k . The minimal N k that produces satisfying results (all the spike types are separated from the noise and among each other), is selected for that electrode. In many cases, satisfying results were achieved for N k larger than the number of spike types'/one for noise. That is because the noise has different forms that are some times more pronounce than small amplitude spikes.

[29] 46w These 2000 events comprise the training set for the automatic classification of events throughout the recording. Each event is classified according to its Euclidean distance from the k centers of clusters in the training set. The cluster that is closest determines the class of the event.

[30] 120w In Fig. 14a we present the centers of the five clusters of spikes found in the training set using k -means. Note the resemblance to the mean spikes presented in Fig. 10. Table 1 presents the detection and sorting results on a test set of 400 spikes and 400 noise events. One hundred of the spikes are in overlap with a nearby spike (distance less than 1.5 ms). An overall of 92% of the spikes in the test set were detected and 87% were classified correctly. There were less than 2% false detections of noise. Regarding the overlapping spikes, 74 of the 100 were detected, 66 of them were classified correctly. The verification of the classification results was performed manually.

[31] 71w 1 46 7 2 51 4 N 2 3 89 1 2 1 79 6 8 1 6 80 5 4 N 3 1 30 4 5 27 8 5 5 15 3 2 15 N 4 92 3 19 5 84 9 16 40 21 3 50 N 5 2 3 81 8 3 14 64 13 5 69 3 3 14 FA 4 3 19 8 1 75 59

[32] 124w Each row presents the classification of neural spikes (named N i where i0 1, . . ., 5), and the last row (FA) presents noise events that were falsely detected and classified as neural spikes. For each sorting method, five columns present the classification distribution among the five classes, and the last column (X ) states the number of spikes that were not detected (classified as noise). The classification was performed using k -means for all methods. For WT we used the nine most energetic wavelets and for PC we used the top two principal components. The results in this table are verified manually. sorting efficiency will be presented elsewhere. We intend to present our work on different classifiers emphasizing the clustering problem elsewhere.

[33] 183w The comparison is performed using the same data that was described in Section 5. In Fig. 16a we show the comparison between the separability of the principal components and the wavelet packets. Separation ability is defined via the capability to capture a higher portion of a spike's energy than of noise energy. Thus, the 'bending' point in the energy distribution curve distinguishes between the separating and non-separating packets (or principal components). We find that only the first two principal components have separation ability. On the other hand, there are nine separating wavelet packets. The first two principal components contain on average 86% of the strongest spikes' energy, while the nine packets capture an average of 84%. However, the principal components contain more noise energy than the packets (Fig. 16b). Using the WT, nine wavelets seem to have separation ability, and they capture 81% of the energy. As will be shown ahead, the fact that only two principal components have separation ability in a data set of five different spikes, will ultimately prevent the PCA method from detecting and sorting all the spike types.

[34] 85w As for the WPD case, the separation of clusters is done using the k -means classifier, and the sorting of a test set is done as described previously. In Fig. 14b we present the centers of the five clusters found in the test set using the principal components. It is clear that the principal components detect three spike templates, but fail to distinguish the two weaker classes from the noise. Using the WT, the centers of the clusters follow those of the WPD (Fig. 14a).

[35] 249w Table 1 presents the comparison of detection and sorting results on a test set of 400 spikes and 400 noise events. One hundred of the spikes are in overlap with a nearby spike (distance less than 1.5 ms). For WT using nine wavelets, an overall of 91% of the spikes in the test set were detected and 75% were classified correctly. There were 7% false detections of noise. Regarding the overlapping spikes, 73 of the 100 were detected, 45 of them were classified correctly. For PCA using the two separating components, an overall of 59% of the spikes in the test set were detected and 42% were classified correctly. There were 34% false detections of noise. Regarding the overlapping spikes, 50 of the 100 were detected, 26 of them were classified correctly. The verification of the classification results was performed according to the a-priori knowledge of the spikes origin. That is, the tested data set was carefully inspected by an independent trained eye in order to determine the ascription of the spikes to the different groups. The superiority of our method compared with the PCA is outstanding in all aspects (detection and sorting of both overlapping and non-overlapping of spikes and in rejection of noise). Compared with the WT, we notice substantial improvement in noise rejection, and in the correct sorting among the five neurons (from 75 to 87%). The contribution is, as expected, especially in the case of overlapping, where the sorting performance improved from 45 to 66%.

[1] 153w We have presented a new method for detection and sorting of neural spikes using WPD. We explained how to select wavelet packets with separation ability among a large number of spikes using the LDB algorithm. We have stated that the performance of this feature extraction method enables us to use a simple clustering method such as the k -means algorithm. We demonstrated the efficiency of our method on in vitro neural networks recordings, and compared the performance to those of the principal components method and regular WT. The main advantage of our method is the ability to detect and sort overlapping spikes. Furthermore, our method is capable of sorting a large number of neurons on a single electrode. The method is general and can be adapted for recordings of in vitro tissues as well as in vivo preparations. The use of more sophisticated classifiers (e.g. support vector machines), will enhance the sorting performance.

[2] 71w During long term recordings, the neurons move with respect to the electrodes. As a result, the shapes of the recorded spikes variate. Our method is capable of adaptively following these variations, and we shall present this concept elsewhere.This feature extraction method also enables us to investigate the inner structure of the neural action potential, as it enables to reliably describe the spike as a linear combination of nine smooth wavelet packets.

[1] 95w We apply our method on measurements of in vitro cell cultures using multi electrodes arrays (Segev et al., 2001). The neurons are randomly attached to the surface. The capacitance coupling between the neurons' membrane and the electrode enables the recording of the spontaneous firing of neural spikes (Jimbo et al., 1993(Jimbo et al., , 2000;;Kamioka et al., 1996;Camepari et al., 1997;Maher et al., 1999). The expression of the neural spikes in the recording is affected by the distance from the center of the electrode and the quality of the coupling (Lewicki, 1998;Rinberg et al., 1998).

[2] 69w The recording of the neuronal activity is held consecutively for days, and we have cases of as many as five neurons read by a single electrode. The signals collected from the electrodes are digitized at a rate of 12 kHz. The WPD of an event is performed on a window of 128 samples (Â/10 ms) with the spike at it's center (as a consequence, L 0/log 2 (N )0/7).

[3] 163w The first step is to select an appropriate mother wavelet. The one we have selected is the 3 rd order Coiflet. Fig. 8 presents different packets of this mother wavelet. These Coiflets are nearly symmetric and have little overlapping with their neighbors both on time and frequency domains. The compactness is important, as it does not spread the signal's energy into neighboring packets (see Mallat, 1998 for discussion). The decomposition of the spike into nine selected packets. Each row presents the time and spectral behavior of a packet (left and middle, respectively), and the packet multiplied by the coefficient of decomposition (hence c 4,0,6 and c 4,0,7 carry opposite signs etc.). (c) The reconstruction of the spike from the nine packets' coefficients. It is clear that the packets grab the main shape of the spike and most of its energy. It is also clear that the usage of the nine packets filters most of the noise in which the neural spike is embedded.

[4] 106w Various signal processing methods have been used to extract the features of neural spikes in the process of detection and sorting (see Lewicki, 1998;Wheeler, 1999 for reviews). One family of methods, is to describe each spike by a set of features. Simple features that are known to be used for describing spikes are the maximal and minimal peaks, the width of the spike etc. Fig. 15 present a histogram of the peak amplitude of events. Those are the same events that comprise the matrix plot in Fig. 13. It is clear that there is no clear separation of the five spike types from the noise group.

[5] 120w Another family of methods, is to filter the data utilizing a series of templates of spike waveforms (Lewicki, 1998). Variants of this method based on the principal components analysis (PCA) have been developed (Lewicki, 1998). The concept of PCA is to transform the data into a new set of coordinates that signify the direction of largest variations in the data. These directions (named principal components) are known to be the eigenvectors of the covariance matrix of the data (Glaser, 1971). The eigenvectors form an orthogonal basis, and the magnitude of the relative eigenvalues designates the importance of the components in representing the data. The top principal components are expected to enhance the differences among the shapes of different spike types.

[6] 137w There are a vast variety of methods used by researchers for feature extraction of neural spikes. We compare our sorting results of our method with the widely used PCA (presented by Lewicki, 1998), and to a method somewhat related to our own, recently presented by Zouridakis and Tam (1997) that uses the WT coefficients rather than WPD coefficients. As will be shown further in this section, our method is superior to these methods both in detecting and in sorting overlapping spikes. Moreover, as stated in the introduction, in this manuscript we do not compare different clustering methods, but rather demonstrate the efficiency of our method merely using the simple k -means classifier. The effect of different clustering methods on the spike Table 1 Performance of detection and sorting of the WPD method compared with WT and PCA