生物神经网络如人脑由称为神经元的特化细胞组成。有各种类型的神经元，但它们都基于相同的概念。称为神经递质的信号分子在突触处释放，突触是两个神经元之间的连接点。神经递质通过与细胞膜内的离子通道相互作用改变突触后细胞的膜电位。如果突触后细胞的去极化足够强，则在轴突小丘处产生动作电位。动作电位将沿着轴突行进并触发神经递质释放到突触间隙中，这将影响下一个神经元的膜电位。这样，信号可以通过（整个）网络从一个细胞传递到下一个细胞，其中动作电位是在突触处释放神经递质的触发器。神经沟通本质上是电化学的，并且知道何时以及在何种条件下产生动作电位可以为大脑的运作提供有价值的见解。

但研究活脑内个体神经元的活动是一项具有挑战性的任务。没有可用的非侵入性方法，通过该方法可以在单个细胞水平上实时监测神经活动。通常将电极插入大脑以记录其附近的电活动。在这些电生理记录中，动作电位表现为快速高振幅尖峰。但由于神经元密集地堆积在大脑内部，因此记录电极通常一次只能从一个神经元中获取峰值。因此，如果我们想要了解单个神经元的行为方式，我们需要从记录中提取峰值并检查它们是由一个还是可能是几个神经元生成的。从数据中提取和聚类峰值称为峰值排序。

开始

首先我们需要数据。Matlab中有一种很流行的spike 排序算法叫做Wave Clus。在他们web页面的data部分（https://www2.le.ac.uk/centres/csn/software），他们提供了一个测试数据集，我们将在这里使用它。根据网页上提供的信息，这段录音大约30分钟，来自一名癫痫患者。数据存储在.ncs文件中，这是制造记录系统的公司的数据格式。因此，如果我们想将记录读入Python，我们需要了解数据是如何存储的。在公司的web页面上有一个关于.ncs文件格式的详细描述，我们可以使用它来导入文件，Python代码如下：

>>> # Define data path

>>> data_file = './UCLA_data/CSC4.Ncs'

>>> # Header has 16 kilobytes length

>>> header_size = 16 * 1024

>>> # Open file

>>> fid = open(data_file, 'rb')

>>> # Skip header by shifting position by header size

>>> fid.seek(header_size)

>>> # Read data according to Neuralynx information

>>> data_format = np.dtype([('TimeStamp', np.uint64),

>>> ('ChannelNumber', np.uint32),

>>> ('SampleFreq', np.uint32),

>>> ('NumValidSamples', np.uint32),

>>> ('Samples', np.int16, 512)])

>>> raw = np.fromfile(fid, dtype=data_format)

>>> # Close file

>>> fid.close()

>>> # Get sampling frequency

>>> sf = raw['SampleFreq'][0]

>>> # Create data vector

>>> data = raw['Samples'].ravel()

正如我们在上面的代码中看到的那样， .ncs文件还包含一些有关录制的附加信息。最重要的是，采样频率告诉我们每秒记录的数据点数。通过采样频率和数据中的采样数，我们现在可以创建一个时间向量，允许我们随时间绘制信号。以下Python代码将执行此操作并绘制信号的第一秒。

>>> # Determine duration of recording in seconds

>>> dur_sec = data.shape[0]/sf

>>> # Create time vector

>>> time = np.linspace(0, dur_sec,data.shape[0])

>>> # Plot first second of data

>>> fig, ax = plt.subplots(figsize=(15, 5))

>>> ax.plot(time[0:sf], data[0:sf])

>>> ax.set_title('Broadband; Sampling Frequency: {}Hz'.format(sf), >>> fontsize=23)

>>> ax.set_xlim(0, time[sf])

>>> ax.set_xlabel('time [s]', fontsize=20)

>>> ax.set_ylabel('amplitude [uV]', fontsize=20)

>>> plt.show()

我们看到尖峰了吗?……不,我们没有。但我们确实在数据中看到了一些有节奏的活动。这是一个有意义的信号吗?答案是否定的。如果我们在这第一秒的数据中计算峰值的数量我们最终会得到60个峰值，这意味着我们看到的是60hz的振荡。我们可以通过绘制信号的功率谱来进一步证实这一点，显示在60hz处有一个清晰的峰值。

那么我们在看什么呢？在我们从中获取数据的网页上，该录音来自美国加州大学洛杉矶分校的Itzhak Fried实验室。美国的电源频率为60 Hz。因此，我们所关注的实际上是数据采集过程中房间内电子设备的电噪声。

找到信号中的峰值

即使数据中有60hz的噪声，我们仍然可以处理它。幸运的是，动作电位是只持续1到2毫秒的快速事件。所以我们能做的就是把原始的宽带信号过滤到一个不含60赫兹噪声的范围内。一个典型的过滤器设置是500到9000赫兹，我们的Python实现如下:

>>> # Import libarys

>>> from scipy.signal import butter, lfilter

>>> def filter_data(data, low=500, high=9000, sf, order=2):

>>> # Determine Nyquist frequency

>>> nyq = sf/2

>>> # Set bands

>>> low = low/nyq

>>> high = high/nyq

>>> # Calculate coefficients

>>> b, a = butter(order, [low, high], btype='band')

>>> # Filter signal

>>> filtered_data = lfilter(b, a, data)

>>>

>>> return filtered_data

使用上述函数处理数据将为我们提供信号的高频带或峰值信道。我们的期望是这个尖峰通道包含动作电位，不再有60 Hz的噪声。因此，让我们看一下滤波后的峰值信道，并将其与原始宽带信号进行比较。

正如预期的那样，峰值通道不再显示60 Hz振荡。最重要的是，我们终于可以看到录音中的第一个峰值。录制大约0.5秒后，在过滤后的数据中清晰可见。此外，现在我们现在可以查看，我们可以在未经过滤的数据中看到它。然而，由于60赫兹的噪音，很难发现。

从信号中提取峰值

现在我们把高频信号带分离出来我们可以提取出单个的峰值。为此，我们将编写一个执行以下操作的简单函数:

• 找到信号中超过某个阈值的数据点
• 在这些事件周围定义一个窗口并“删除它们”
• 将它们对准它们的峰值振幅

此外，我们还将定义一个上限。超过此阈值的数据点将被拒绝，因为它们可能是高频工件。这样的人为现象可能是由于病人的运动引起的，也可能是由于房间里的电灯泡的开关等电气事件引起的。您可以在Jupyter笔记本（https://github.com/akcarsten/akcarsten.github.io/blob/master/spike_sorting/Spike_sorting%20.ipynb）中详细查看脉冲提取函数。这里我们只看一下100个随机尖峰它们是用函数从信号中提取出来的。

在上面的图中，我们可以看到数据中至少有两种波形。一组高振幅的峰值，另一组具有更大的初始峰值。所以这些尖峰很可能是由多个神经元产生的。