SPKtool
A tool for neurophysiological data processing
锋电位分析工具
版本1.0 2008-8 1 简介及安装
1.1简介
SPKtool是一个在MATLAB平台上编写的神经元锋电位分析工具软件，提供了从锋电位检测、特征提取、聚类分析到放电序列分析的各种功能：
锋电位检测：
-
显示记录信号中锋电位峰值的分布曲线；
-
通过阈值检测锋电位；
-
通过原始信号、能量、非线性能量进行检测；
-
三种锋电位波形对齐方式：阈值、波谷、波峰。
-
GUI交互式界面。
特征提取：
-
多种常用特征：PCA，峰值、能量等；
-
小波方法；
-
特征分别图，特征密度分布图等。
聚类分析：
-
K-means；
-
模板匹配；
-
基于高斯混合模型的EM算法；
-
Valley seeking；
-
通过绘制轮廓线，手动聚类；
-
对聚类结果进行合并，删除等。
动作电位串分析：
-
相邻动作电位间隔分布图(ISI)；
-
Poincare map；
-
相关函数图；
-
锋电位分布光栅图；
-
放电频次直方图；
-
刺激前后放电直方图（PSTH）；
-
刺激前后放电光栅图。
1.2安装及界面
本程序在MATLAB 7.1 SP3上开发完成，可能需要该版本或更新的版本的MATLAB才能运行。
将程序压缩包拷贝到MATLAB中的work目录下，解压缩，产生SPKtool目录。在MATLAB中，选择菜单File->Set Path … ，在“Set Path”对话框中选择“Add with Subfolders…”，把该目录及子目录添加到MATLAB的路径中去，选择“Save”按钮保存。
在命令窗口输入：spktool，打开GUI界面如下图。
程序界面中包含三个部分，菜单、信息表格和底部的状态栏。软件的所有功能均通过菜单命令来操作。信息表格列出了分析数据的主要信息，如采样频率，记录长度，文件名称、锋电位数目、特征名称、各个分类包含的锋电位数目等。对于多通道数据，通过信息表格右上角的Channel下拉菜单选择当前通道。状态栏给出程序运行状态的简明信息。
2 数据文件
SPKtool没有和其它软件进行数据交换的接口，不能直接读入记录系统产生的数据。用户需要首先将数据（连续记录信号、锋电位波形矩阵等）导入MATLAB，并按照特定的变量名保存到.mat文件（MATLAB中默认的数据文件格式），才能在SPKtool中打开，进行分析。
在菜单中选择File->Load…，打开文件选择对话框。加载数据后，根据文件中的变量，自动设置各个菜单的状态。
SPKtool通过变量名区分数据，因此，对数据文件(.mat)中的变量名有严格的要求。
2.1原始记录数据
原始记录数据文件中至少应包含两个变量：
dataflow：连续记录数据，二维数组，行对应通道，列对应采样点；
fs： 采样频率，数值。
默认情况下，SPKtool对读入的数据进行低通滤波，使用6阶Butterworth滤波器，截至频率4500Hz。可以打开MLoadRaw.m，找到butter函数对应的位置，修改滤波参数。
2.2锋电位波形数据
从原始信号中检测出的锋电位波形数据，至少应该包含以下三个变量：
posind： 锋电位波形在连续数据中的采样位置标记。行向量。第n个波形发生的时刻为posind(n)/fs。
waveforms： 锋电位波形数据，二维数组。每行对应一个锋电位的采样数据。例如，1000个波形，每个波形含40个数据点， 则waveforms尺度为1000×40；
fs： 采样频率。
数据文件中还可以包含如下一些变量：
featurespace： 锋电位特征空间，二维数组，每行包含一个波形的特征。如1000个波形的PC1和PC2特征，则featurespace为1000×2。特征的数目不限。
featnames：特征名称，包含字符串的单元数组，每个字符串对应一个特征的名称。例如: 两个特征的名称分别为PC1和PC2， 则 featnames 为{‘PC1’,’PC2’};
unitcell： 包含分类信息的单元数组。每个单元对应一个类别中的锋电位序号列表。例如，unitcell{1}对应没有归类的锋电位编号。 unitcell{3} 是第二个类别中的锋电位编号数组。
对于多通道数据，上面这些变量需要放在一个名为ResultSet的结构体数组内，每个结构体对应一个通道。
刺激信息
stimmark: 刺激的开始和结束时间。行向量。例如[20 30 60 70]，表示两次刺激，第一次刺激从第20s开始，到30s结束，第二次刺激60s开始，70s结束。
上述变量可以共存于一个文件中。在程序中，如果找不到dataflow这个变量，将关闭锋电位检测功能（不能激活Detect菜单），找不到waveforms，将关闭特征提取功能（不能激活Features菜单）。
另外，程序中用到了两个全局变量，采样频率fs和定义各个类别色彩的矩阵colorary，所以在程序运行过程中，请勿执行clear all等能够清除全局变量的命令。
2.3数据保存和导出
选择菜单File->Save …，保存分析结果。文件格式同样是.mat文件，变量名称与读入数据相同。（注意，连续记录数据不包含在该数据文件内。）可以用SPKtool重新打开保存的文件，也可以直接用读入MATLAB，将特定的数据导出。
选择菜单File->Save Current Channel…，仅保存当前通道的分析数据。
菜单File-> Experimental Info …，设置试验信息，包含记录位置、时间、坐标、数据文件名称等信息。
菜单File-> Export to TXT file …，导出数据到文本文件。这将产生两个文本格式的文件，后缀名为.inf的文件包含上面的试验信息，以及数据信息（通道数、单元数等）。后缀名为.txt的文件是一个数据表格，第一列对应锋电位时间（即变量posind），第二列为类别编号，后面的各列为波形数据。
3 锋电位检测
如果读入的数据中包含原始连续记录信号（即dataflow变量），菜单Detect可用，其中包含以下两个命令：
3.1峰值密度分布曲线
菜单Detect->Plot Peak Histogram，绘出原始记录数据的峰值密度分布曲线。幅度为对数坐标。
3.2锋电位检测界面
菜单Detect->Load Detector…，加载锋电位检测，界面如下：
参数设置：
Detect On：该下拉菜单中包含了三种波形检测方式，即分别在原始信号（Raw Signal），原始信号的能量(Energy)和非线性能量(Nonlinear Energy)1曲线上进行检测；
Threshold：阈值。输入阈值后左图中的红色水平线也相应调整。信号中与阈值线有交点的波形将根据后面的参数进行截取；
Spike length：每个锋电位波形的长度，单位为ms。如果采样频率为40KHz，1ms长度的锋电位波形包含40个数据点；
Prethreshold ratio：对齐点以前数据长度比例。对齐点前保留的数据长度与波形数据总长度的比值。例如包含40个点的锋电位波形，对齐方式为valley时，在波谷位置前保留10个采样点，则该参数需要设置为10/40=0.25；
Dead time：不应期长度，单位为ms。例如，如果采样频率为40KHz，该参数设置为2时，程序在检测到一个波形后，将跳过该波形后紧邻的80个点，然后才开始寻找下一个波形。
Alignment：波形对齐方式，阈值位置（gate）、波谷 (valley)、和波峰 (peak) 对齐三种形式。阈值和波谷的对比如下：
a) gate b) valley
使用方法：
中间的坐标轴根据Detect On的设置，显示原始信号或其（非线性）能量曲线。底部坐标轴在（非线性）能量检测方式下显示原始信号波形。左侧坐标轴显示中间坐标轴中数据的峰值分布曲线（横坐标为对数坐标）。
可以利用顶部的工具栏，对中间坐标轴进行缩放和平移，左侧坐标轴的纵坐标和底部坐标轴的横坐标也将进行相应的调整。
参数设置完毕后，选择Detect进行波形检测，检测完成后，将打开一个单独的绘图窗口显示锋电位波形。如果不符合要求，可以修改参数，重新检测。最后，选择Done完成，并关闭该界面。
需要说明的是，由于原始记录信号中的数据量通常很大，为了加快运行速度，坐标轴中只绘出前0.4e6个采样点。当然，锋电位波形检测还是针对所有的采样数据。
右侧的参数中除Detect On外，具有自动保存功能，再次运行Detector时将自动加载前一次设置的参数，即使你前一次选择Cancel退出。
在切换不同的检测方式时，阈值将根据中坐标轴中的数据自动调整为：均值+/-3倍标准差。
4 锋电位特征提取
通过阈值检测出锋电位波形后，或者读入的数据中保护锋电位数据，这时，Feature菜单将可用。在该菜单中，可以进行锋电位特征的选择、计算、修改、绘图等功能。通常需要定义特征后，才能打开Sort菜单进行聚类操作。
4.1特征定义
菜单 Features->Feature Define…，打开如下的特征选择对话框：
左侧All Features 列表中列出了所有可用的特征。中间有两个特征列表：Features Vs Features 中的特征用于绘制特征比较图（菜单 Features-> Features Vs Features）；Features For Cluster 才是用于聚类分析的特征。本说明中，除非特别指明，特征均指用于聚类的特征。在右侧可以修改一些特征相关参数。
可用特征：
peak：锋电位波峰处的幅值；
peak.pos：峰值对应的位置；
valley：锋电位波谷处的幅值；
valley.pos：波谷对应的位置；
peak-valley：峰峰值，波峰到波谷的振幅。
width：波形宽度，即波峰位置到波谷位置的距离：peak.pos-valley.pos；
energy：波形的能量，即采样点的平方和。
nonlinear.energy：所有点的非线性能量之和，幅度平方与频率平方的乘积，其计算方法为：NEO(t) = v(t)2- v(t-1)*v(t+1)，（v(t)为记录信号）；
Slice：锋电位波形中某个时刻的幅值，对应参数Slice Parameters。例如每个锋电位含32个采样点，Slice Parameters为0.5，则Slice特征为锋电位中第32×0.5=16个采样点的数值。Slice Parameters可以包含多个值，提取多个采样点，例如 0.5，0.75 则第16和第24个采样点的数值都作为特征。
pca：主成份分析方法，可以在右侧PCA Parameters中设置要提取的主成份序号；
ridge.pca：基于小波脊的PCA特征提取方法，首先对波形逐个进行小波变换（使用mexh小波），然后提取出小波脊，再对小波脊上的小波系数进行PCA，提取特征。需要三个参数：小波分析的频率范围Frequency limits，小波脊的数目Number of Ridges，以及需要提取的主成份序号PCA Parameters。
完成特征选取后，选择Done，开始提取用于聚类分析的特征，并给出特征分布散点图。
4.2特征修改
特征比例缩放
Features->Scale Features，打开如下图的对话框，在其中，可以调整特征的相对范围，可能有利于改善K-Means的聚类结果。
删除特征
Features->Remove Features，删除部分使用中的特征。一般特征的计算不需要花很长时间，重新选择并计算即可。在同时使用Ridge.pca和PCA作为特征时，由于共用PCA参数，可以通过该工具重新选择。
4.3特征分布散点图
Features->Plot Features 绘制特征分布散点图，如果定义了两个特征，散点图是二维的，如果有三个或更多的特征，则只使用前三个，以三维形式显示，并同时画出特征向三个平面上的投影分布，如下图。
4.4特征密度图
Features-> Features Density：特征分布密度图。如下图，该图只能绘出前两个特征的分布密度。
4.5特征对比图
Features-> Features Vs Features，绘制特征分布的对比图，注意该功能中使用的特征与聚类分析的特征是独立的。并且，每次绘图都需要重新计算特征。
如果已经划分类别，不同类别的特征将以不同颜色区别显示。
4.6主成份
Features-> Principal Components：给出特征提取中主成份曲线。例如特征定义中使用第一个和第二个主成份，则给出PC1和PC2的形状曲线。
5 聚类分析
聚类功能集中于菜单Sort中。在选择具体的聚类方法之前，应该先设置聚类的参数。
5.1参数设置
菜单Sort->Options…, 打开参数设置对话框：
其中共有四个参数：
Unit Number：划分类别的数目，K-Means，期望最大化(EM)方法，以及模板匹配方法使用该参数；
K-Means Radius：K-Means方法的聚类半径，值越小，聚类的半径越小。模板匹配方法也受该参数影响；
Parzen Mult：寻谷聚类算法的初始半径大小，初始半径越小，获得的类别数目越多。
Sort On：设置聚类的锋电位范围，选择All Spikes时对所有锋电位的特征进行聚类；选择Unsorted Spikes时仅对未分类的锋电位进行聚类。
5.2手动聚类
菜单Sort->Contour…，对应手动聚类，选择该菜单后，弹出如下图的手动聚类对话框。
对话框中间给出了未分类的特征的分布图。底部有四个控制按钮。首先单击New Contour按钮，光标将变成图中所示的十字形，在特征分布图中单击左键，将出现一个红色的十字形标记点，这样依次点击，就可以确定一个类别的大致轮廓，之后，鼠标右击，这些点将自动连接成一条封闭的曲线。然后选择底部的Update按钮，曲线内的点就被划分到一个类别。再选择New Contour绘制下一个类别的边界，依此操作，直到聚类完成。
也可以连续绘制多条边界线，然后再选择Update，显示聚类结果。最后选择Done将显示的聚类结果传回主程序，选择Cancel放弃。
手工聚类的类别数目不受前面聚类参数对话框中的类别数目的影响。
5.3自动聚类方法
SPKtool中包含了几种常见的聚类算法，简单说明如下：
K-Means
根据特征之间的距离来进行聚类。用户指定聚类的数目（聚类参数对话框中），该算法首先从样本（特征）中随机抽取确定聚类中心，然后计算每个样本到各个聚类中心的距离，并将其归到距离最近的类中，第一次聚类完成后，计算各个样本的均值作为新的聚类中心，重新计算每个样本到新的聚类中心的距离，并将其划分到聚类最近的类别中，如此循环，直到聚类中心的位置稳定[2]。
Template 模板匹配法
在SPKtool中，模板方法的与K-Means基本上是一样的，不同之处仅在于直接用锋电位的波形数据作为特征，而K-Means使用锋电位降维后的特征数据。
Gaussian EM 基于高斯混合模型的期望最大化算法
该算法首先用K-Means方法确定初始聚类，然后将锋电位特征的分布看成符合混合高斯分布，用期望最大化算法确定高斯混合模型的均值、协方差等参数，然后计算每个特征在各个高斯分布中的概率，并划分到概率最大的一个类别（高斯分布）中[2]。
Valley Seeking 寻谷聚类算法
在特征空间中寻找特征分布密度的极小位置，作为聚类的边界。算法的流程可以简述如下3:
1. 确定初始半径，可以取为样本标准差的某个倍数；初始半径的设置影响聚类的数目；
2. 在该半径范围内，确定每个特征的近邻数目，以及近邻的集合。如果特征的数目为N，需要初始化一个N×N的数组来存放每个特征的近邻编号，计算完成后，可以根据最大近邻集合中元素的数目，删去多余的元素。
3. 如果一个特征的近邻数目，超过它每一个近邻各自的近邻数目，则该特征确定了一个特征密度峰。如此，设置一个最小近邻标准，可以找出全部峰。
4. 五个峰确定5个类别，给每个峰的近邻特征设置类别标记；
5. 将所有特征按照近邻数目从大到小排序；
6. 从近邻数目最大的特征开始，看其近邻中属于哪一类别的最多，就将该特征也归入这一类；如此，完成全部特征的聚类。
5.4聚类结果调整
移除未分类的锋电位
菜单Sort->Edit units …中，选择复选框Remove unsorted spikes。未分类的锋电位将从数据中删去。
合并分类
同样在菜单Sort->Edit units …中，在Merge Units框中输入欲合并的单元编号；例如，输入“2 3”将单元2和单元3合并为单元2，输入“0 3”删除单元3（其中的波形转为未归类状态）。
移除特征离散度很大的锋电位
可以首先通过轮廓线聚类方法，将要保留的特征划分到单元1中，然后按照上面的方法，先移除没有归类的边缘特征，再将单元1与0（未分类）合并。
6 绘图功能
6.1锋电位波形
菜单：Plot->Spikes：绘出所有锋电位波形；
菜单：Plot-> Sorted spikes：分布绘出各个类别的波形。每个波形图左下的数字标明了该分类中波形的数目。
6.2特征分布图
菜单：Plot-> Features：绘出特征分布散点图；
菜单：Plot->Sorted Features：画出分类后的特征散点图；如下图。
6.3时间间隔(ISI)分布
菜单：Plot->ISI …，打开如下的ISI参数设置窗口：
Min. Interval: 最小间隔时间；
Max. Interval: 最大间隔时间；
Number of Bins: 阶梯图的分段数目，数值越大，ISI曲线越平滑；
Units: 选择单元，只画出指定单元的ISI分布图。
6.4 间隔-时间分布图
菜单：Plot->ISI versus time…，绘制相邻放电间隔的时间分布图。参数设置对话框：
X Min. :横坐标（时间）最小值；
X Max.：横坐标最大值；
Units：选择需要进行绘图的单元编号。
6.5庞加莱图
菜单：Plot-> Poincare maps：该图以与前一个放电的时间间隔为纵坐标，与后一个放电的时间间隔为横坐标，标出每一个放电的位置。坐标轴的范围（即最大时间间隔）由ISI中的Min. Interval和Max Interval 指定。如下图。
6.6放电频次阶梯图
菜单：Plot->Fire Rate…：绘制放电的频次直方图。参数设置对话框：
X Min. :横坐标（时间）最小值；
X Max.：横坐标最大值；
Time per Bin：每个阶梯表示的时间长度；
Units：选择需要进行绘图的单元编号。
Add stimulation markers：选中该项后，可以在下面的Markers中输入刺激标记（对应数据文件中的simmark变量）。例如，刺激开始的时刻为49.89885s，结束的时刻为262.4599s，在Markers中填写如下：
放电频次分布直方图如下，刺激时间将以一个黑框标示出来。
6.7放电光栅图
菜单：Plot->Rasters …：绘制放电的光栅分布图。参数设置对话框：
X Min. :横坐标（时间）最小值；
X Max.：横坐标最大值；
Units：选择需要进行绘图的单元编号。
Add stimulation markers：与频次直方图中相同。
6.8相关函数图
菜单：Plot->CrossCorrelograms…：绘制单元间的互相关和自相关函数。参数设置对话框：
Pre-ref. Time：横坐标（时间）最小值，计算相关函数时信号的偏移范围；
Post-ref. Time：横坐标最大值；
Number of Bins：相关函数分段数目；
Units：参与计算的单元编号；
Ref Units：作为参考的单元编号。
6.9刺激前后放电直方图PSTH
菜单：Plot->Peristimulus Histogram … ：打开如下的参数设置对话框：
T. pre-Event：刺激事件前保留的时间长度；
T. post-Event：刺激事件后保留的时间长度；
Time/Bin：直方图的每段包含的时间长度；
Units：选择需要进行绘图的单元编号。
Stimulus Phase：选择刺激事件的相位，“on”表示刺激开始；“off”表示刺激结束；
Stimulus No.：当有多个刺激周期时，可以选择其中的一个刺激周期。设置为0时，将所有刺激周期的的PSTH叠加。
6.10刺激前后放电光栅图
菜单：Plot->Peristimulus Raster … ：打开如下的参数设置对话框：
T. pre-Event：刺激事件前时间长度；
T. post-Event：刺激事件后时间长度；
Units：选择需要进行绘图的单元编号。
Stimulus Phase：选择刺激事件的相位，“on”表示刺激开始；“off”表示刺激结束；
Stimulus No.：当有多个刺激周期时，可以选择其中的一个刺激周期。设置为0时，将所有刺激前后放电都画出。
在该组信号中，共有三个刺激事件，Stimulus No.设置为0时，三个刺激前后的放电都绘出，如下图。纵坐标名称中，E*表示事件的编号，U*表示单元编号：
部分参考文献：
1Mukhopadhyay S, Ray GC. A new interpretation of nonlinear energy operator and its efficacy in spike detection. IEEE Transactions on Bio-Medical Engineering. 1998 Feb ;45(2):180-7.
2 （美）迪达等 著，李宏东等 译. 模式分类（原书第2版）[M]. 北京: 机械工业出版社，2003.
3 Zhang C, Zhang X, Zhang MQ, Li Y. Neighbor number, valley seeking and clustering. Pattern Recogn. Lett. 2007 ;28(2):173-180.
其它文献：
梁培基，陈爱华. 神经元活动的多电极同步记录及神经信息处理 [M]. 北京：北京工业出版社, 2003.
Lewicki MS. A review of methods for spike sorting: the detection and classification of neural action potential [J]. Net Comput Neural Syst, 1998, 9(4): 53-78.
Answer.com: Expectation-maximization algorithm
