基于概率神经网络的锋电位实时分类算法*

祝晓平 韩业强 郝耀耀 王东 陈耀武†

(1.浙江大学生物医学工程与仪器科学学院，浙江杭州310027;2.浙江大学求是高等研究院，浙江杭州310027)

脑－机接口是一种将脑电信息转化为外围设备控制命令的技术，它在生物体与外部环境之间建立了一种新型的信息交流与控制通道，从而为残障人士的运动控制提供了一种新方法［1］.脑－机接口技术现已成为神经科学与信息工程技术交叉学科的热点课题［2－4］，其中运动信息的神经解码是一个重要的研究方向［5］.

锋电位分类是运动信息神经解码研究中的关键问题之一.大脑中的神经元通过动作电位即锋电位进行信号的传递、交流和处理.每个电极通常会把周边若干个神经元的锋电位记录下来［6］，因此，要准确获取生物神经网络所传达的信息，就需要将各个神经元的锋电位区分开来［7］.在实际的脑－机接口实验和应用中，锋电位实时在线分类不但可以减轻动物脑电信号采集部分和神经信号解码部分之间传输数据的负担，而且是实现实时解码的先决条件，因此，实时、高效、准确地完成神经元锋电位的分类是实用脑－机接口设备的必要条件.

传统的锋电位分类算法包括特征提取和聚类两个步骤.特征提取是指从锋电位信号中提取特征值的过程，其常用算法主要有主成分分析法(PCA)［8－9］、模式匹配［7，9－10］和小波分析［10－11］等.聚类是将锋电位按特征值进行分类的过程，通常采用k均值聚类算法［8］.然而，k均值聚类算法需要知道分类的数量，并且对初始值的选择很敏感［7］.Yang等［7］使用数字集成电路完成锋电位的检测和特征提取，但聚类算法仍然在个人计算机上实现，因而该算法没有被广泛地应用;Zumsteg等［8］采用噪声白化的主成分分析算法和最大后验分类算法来实现锋电位的实时分类，但需要较长的训练时间;Dai等［9］实现了模式匹配、主成分分析法、径向基函数神经网络(RBFNN)、反向传播神经网络(BPNN)和小波分析5种锋电位分类方法并对其性能进行了比较，结果显示模式匹配方法具有最佳的性能，而反向传播神经网络的性能最差;Vollgraf等［10］针对封闭形式的约束二次优化问题，提出了一种基于时域优化的最佳线性滤波算法，该算法比基于频域的滤波算法具有更好的分类效果;Hulata等［11］采用小波包分解(WPD)方法来实现锋电位的分类，能够更好地处理锋电位重叠的情况，但在小波包选取上没有统一的方法.以上算法均注重算法的分类性能，而忽略了实时在线分类的性能指标.商用Cerebus 128通道数据采集系统(美国Cyberkinetics Inc.公司生产)采用模式匹配方法来实现锋电位的实时在线分类，但在模式的选择上需要人为设定，且设备体积过于庞大，限制了实验中动物的活动范围.为了实现便携式实时处理的脑－机接口，需要提出一种新模式来完成锋电位分类的实时计算，并达到低功耗、便携式实时处理的设计要求.

有鉴于此，文中提出了一种新型的基于概率神经网络(PNN)的锋电位实时分类算法.该算法将特征提取和聚类问题统一到PNN中，快速训练PNN后，根据训练数据求得贝叶斯准则下的最大后验概率值，从而获得输入锋电位分类信息.文中还完成了该算法的FPGA实现，并结合嵌入式系统设计技术以达到低功耗、便携式实时处理的设计要求.

1 概率神经网络

Specht［12］在1990 年首次详细介绍了 PNN，将贝叶斯估计放置于一个前馈神经网络中，依据最小风险的贝叶斯决策实现最优分类，其优点是训练简单、迅速，只要有足够的训练数据，就可以保证获得贝叶斯准则下的最优解.PNN的应用领域包括分类问题、映射问题、归纳概率估计问题等.Zhou等［13］成功地将PNN应用于大鼠压杆试验中的神经解码.Vinitha等［14］实现了基于PNN架构的人脸识别系统.

文中采用的PNN结构如图1所示，共分为4层:输入层、模式层、求和层、输出层.图中，X为测试向量(即一个未分类的神经锋电位向量)，X=(X1，X2，…，Xn)，Wi为分类 Ci的训练向量(即具有分类号的神经锋电位向量)，Di为分类Ci中的训练样本数，Si为分类Ci对应的分布密度.模式层是由m类训练向量样本构成，第i类包含Di个训练向量.将分类好的训练向量填充模式层，即可完成该PNN的构建.

测试时，输入层把测试向量X送至所有的模式层单元.每个模式层单元计算测试向量X和该单元中的训练向量Wi之间的距离，然后进行非线性运算，对应的激活函数为

图1 概率神经网络的结构Fig.1 Architecture of PNN

式中，σ 为平滑系数，(Wi－X)(Wi－X )T为测试向量X和训练向量Wi之间的距离.对输入向量和训练向量做归一化处理，则激活函数可以进一步简化成

其中Zi=X·Wi.每个模式层单元的结构见图2.

图2 模式层单元的结构Fig.2 Architecture of unit in pattern layer

求和层的单元个数等于类别数，各单元把模式层单元的输出按照类别号相加求和，获得各分类对应的分布密度.在文中的锋电位分类应用中，由于微电极在采集神经元信号时，单个微电极一般能够采集到4个左右的神经元信号，因此将求和层的单元个数设定为4.输出层从求和层中求得输入锋电位相对于各类别号的分布密度最大值，并将其对应的类别号作为PNN的输出.这样就完成了该测试向量即一个神经元锋电位的分类.

2 算法的性能分析与FPGA实现

2.1 生物手术及实验

文中实验系统的组成如图3所示，包括生物部分和信号处理部分.生物部分包括动物手术和行为学训练;信号处理部分包括数据采集、锋电位检测和锋电位分类.

图3 实验系统示意图Fig.3 Schematic diagram of experimental system

动物实验的所有操作均遵守国家实验动物管理和使用的有关规定.实验采用体重为280g左右的成年雄性Sprague Dawley大鼠.手术前实验动物用戊巴比妥钠(300mg/kg)腹腔注射麻醉后，置于脑立体定位仪(美国Stoelting Inc.公司生产)上，参考大鼠脑图谱，在左侧初级运动皮质层(M1)植入2×8的慢性阵列电极(排间距为0.4mm，排内间距为0.3mm，电极丝材料为包裹Tefloen的镍镉合金，电极尖端横截面裸露，美国Micro Probes Inc.公司生产)，电极尖端的位置在前肢区的第5层.电极埋置完后利用牙科水泥固定电极并封住伤口.

动物训练中，老鼠需要通过压杆并超过一定的阈值来获得水的奖赏.为了保证实验的顺利进行，老鼠在不训练时仅给少量水，使其处于口渴状态.当它能够较为熟练地压杆喝水时，对其进行手术.手术后大鼠恢复一周，再将其置于实验箱中继续进行压杆喝水实验.共有5只老鼠参加通过压杆来获取水的训练，其中3只老鼠各自进行了手术，植入16通道的微电极阵列.由于植入电极会对神经元造成损伤，因此个别电极通道存在没有数据的情况，文中实验只获得36个微电极通道(按顺序标注为1－36)的数据，并作为实验数据源.实验期间，使用Cerebus 128通道数据采集系统对大鼠的神经数据以30 kHz的采样率进行等间隔采集，压杆的压力数据通过压力传感器也以500 Hz的采样率被记录下来.由于Cerebus 128通道数据采集系统内置锋电位的实时检测和分类，因此不仅可以直接记录分类前的神经数据，而且可以记录分类后的神经数据.分类后的神经数据被作为标准分类参考来进行实验效果的对比，分类后的锋电位信号可用于神经解码以获取大鼠的压杆意图.

2.2 算法性能分析

选用实验数据对文中算法、k均值聚类算法、基于PCA的锋电位分类算法进行测试并比较了它们的性能.3种算法的锋电位分类结果如图4所示.

图4 3种算法的分类结果Fig.4 Sorting results of three algorithms

由实验结果可知:文中算法、k均值聚类算法和基于PCA的分类算法的平均分类准确率分别为93.86%、93.33%和79.96%，文中算法和k均值聚类算法的分类准确率高于基于PCA的分类算法，且文中算法的分类准确率略高于k均值聚类算法.

基于PCA的分类算法是通过PCA获得贡献率总和达到95%的前几个主成分，并用于锋电位的聚类分析，然后采用最大最小距离法进行聚类.在采用最大最小距离法确定聚类中心时，判定是否存在第3个聚类中心的条件为［15］

式中:dj1和dj2分别为第j个样本的锋电位与第1和第2个聚类中心的距离，d12为第1、第2个聚类中心之间的距离;m为检验参数，文中取值为1/2.依此类推来确定是否存在新的聚类中心.

从算法自身来看，k均值聚类算法需要明确分类的类别数k，同时需要通过标准测度函数的收敛来确定聚类的中心，这需要反复迭代，且k均值的初值选择是随机的，可能因聚类中心的初值选择不同而出现标准测度函数不收敛的情况，因此存在不确定因素.在锋电位分类的应用中，确定聚类中心的过程可称为训练过程.从k均值聚类中心的确定步骤来看，其训练时间并不可控，且实验中k均值聚类算法需要进行多次聚类才能获得最佳的分类效果.相比之下，文中算法仅需要将训练数据进行存储即可完成训练过程，简单且迅速，需要极短的训练时间.一旦完成训练，之后的分类将和训练样本的分类保持一致，这样便于对算法的分类效果进行评估.但k均值聚类算法因每次的初始聚类中心为随机选择，所以会出现k均值的分类类别号与用于评估的标准类别号不一致的情况，这虽然不影响最终的分类，但会给算法性能的评估带来不便.对于基于PCA的分类算法，在用PCA进行降维操作后，采用最大最小距离聚类时，要求有两种以上的分类类别，否则会将一类的类别分为两类，形成错误分类，而文中算法不会出现这种情况.

文中算法是针对单个通道的锋电位进行分类的，而在植入式脑－机接口中，微电极阵列中的电极数量较大，需要并行完成多通道锋电位的并行分类工作.为此，文中还给出了所提算法的基于FPGA的实现，以便在完成多通道锋电位的并行分类工作时，大幅提升文中算法的计算速度.

2.3 FPGA的设计和实现

文中使用Xilinx公司的ML506评估板作为实现锋电位分类算法的硬件平台，其FPGA型号为Virtex－5 系列的芯片 XC5VSX50T－1FFG1136，适于处理复杂的数字信号，使用的语言是Verilog HDL.文中算法的FPGA实现结构如图5所示，共有5个模块:控制模块、训练模块、距离计算模块、指数计算模块、求和比较模块.

在实验中，每个电极的锋电位被分为1－4类，当信号超过锋电位检测设定的阈值时，则认为有锋电位产生，就将前12个采样点、当前采样点及后面的11个采样点数据作为该锋电位的波形数据.这些参数的值可以根据实际情况进行调整.

2.3.1 控制模块

控制模块包括逻辑控制和总线接口.逻辑控制负责内部有限状态机的改变以及数据流的走向;总线接口负责和外部主机之间的交互，在交互部分使用了若干个握手协议.外部主机首先检查FPGA的内部操作状态，再对特定的地址发送命令，然后等待FPGA的中断信号，确认操作是否可以执行.这种交互机制使得系统具有可配置性，且能够提高交互的可靠性.

2.3.2 训练模块

训练模块把从控制模块接收到的数据按其类型进行存储.为了实现并行的距离计算，训练模块用将24个采样点的训练向量分别存放到24个RAM块中，这种处理方式可以保证在一个时钟周期内完成对24点数据的读取.此外，还需要将和训练向量对应的分类号数据存放到单独的RAM块中，为最终的分类输出做好准备.训练向量为100个锋电位信号，每个信号用24个采样点来表示，每个采样点用单浮点精度来表示，同样对应的分类数据为100个单浮点精度的分类号数据，因此，训练模块占用的RAM块大小为10 kB.每个锋电位和与之对应的分类号存放将消耗25个时钟周期，因而整个训练模块将占用2500个时钟周期.

图5 概率神经网络的FPGA实现结构Fig.5 Architecture of FPGA implementation of PNN

2.3.3 距离计算模块

距离计算模块完成输入测试量和训练量之间的距离计算:

式中，Xi为测试量，Wij为训练量，平滑系数σ=15.56.由于该模块中涉及了大量的浮点计算，因此，需要采用一种并行计算结构(见图6)来加快算法的计算速度.首先，使用24个浮点减法器来完成测试量和训练量对应点的并行减法计算;接着，使用24个浮点乘法器完成减法结果的并行乘法计算;然后，使用23个浮点加法器完成乘法结果的累加计算;最后，使用1个浮点乘法器完成累加结果和平滑系数的乘法计算.该模块的浮点计算器是调用浮点计算IP［16］，每个浮点计算器需要调用2个FPGA片上的DSP48Es资源，每次浮点操作需要消耗一个时钟周期，由于采用了流水线的操作模式，那么需要8个时钟周期的延时，因此完成输入测试量与所有训练量(100个)的距离计算将需要108个时钟周期.

2.3.4 指数计算模块

指数计算模块完成PNN激活函数(式(1))的计算，通过将自变量分解为整数和小数两部分再分别进行计算来求解:

式中，z为整数部分，f为小数部分.对指数函数而言，当自变量x大到一定值时，exp(－x)的值就接近于0，文中设定该值为31，因此，在计算整数部分expz时，只需从含有32个元素的查找表中进行查找即可.

对于小数部分exp f，将其表示成双曲正弦函数和双曲余弦函数的和:

则 exp f的值可通过 CORDIC 方法［17－18］来计算.文中使用的CORDIC IP核具有两个局限性:(1)只能对定点数据进行处理，因此需要在浮点格式与定点格式之间进行转换;(2)输入变量定义域为，因此需要把小数部分的值分为两部分来计算，即

图6 距离并行计算结构Fig.6 Architecture of parallel calculation of distance

图7 指数运算流程图Fig.7 Flowchart of exponent calculation

指数部分的整体运算流程图如图7所示.输入自变量来自前一级距离计算的输出dj，指数部分的输出结果为exp(－dj).在获取自变量后，首先完成自变量的分解，获得自变量的整数部分和小数部分.对于整数部分，判断是否大于31，若是则将溢出标志位置位，整个指数运算结果输出为0;否则通过查找表获得整数部分对应的函数值，输出为该函数值.对于小数部分，通过标志位来选择第1部分小数的值为exp0或者为exp0.5，将第2部分小数的值(不大于0.5)通过浮点转定点获得CORDIC IP核可以接受的数据格式，并通过IP核计算求得双曲正弦函数和双曲余弦函数的值，再将定点格式转换为浮点格式，相加即可获得第2部分小数的值.

然后，将整数部分、第1部分小数值和第2部分小数值相乘并取倒数，可获得计算结果.由于此模块涉及了一系列的转换、乘、加操作，所以需要通过缓存将3部分的值对齐.

此模块的计算占用了2个浮点乘法器、1个浮点加法器、1个浮点除法器、1个浮点转定点和2个定点转浮点，共占用了6个FPGA片上的DSP48Es.其中CORDIC的计算延迟为28个时钟周期，其它计算延迟均为1个时钟周期，相关逻辑操作的延迟为6个时钟周期，故此模块的计算延迟为41个时钟周期.

2.3.5 求和比较模块

求和比较模块负责将指数运算结果按类求和，并从所有的和中求得最大值及其对应的分类号.如图5所示，先将指数运算结果存放在指数结果存储器中，对不同的分类分别进行求和，最后将4个求和中的最大值及其对应的分类号作为最终的分类结果输出.此模块需要占用400 B的RAM资源，用于存放100个指数运算结果，并且需要1个浮点加法器和3个浮点比较器.加法操作需要1个时钟周期，取数操作和时序控制需要4个时钟周期，因此，对100个指数运算结果按类进行累加需要500个时钟周期.3个比较器的延迟为6个时钟周期，所以该模块共需要506个时钟周期.

综上所述，文中PNN的训练需要2500个时钟周期，输入的每个锋电位分类计算需要655个时钟周期.

3 实验结果与分析

FPGA的设计软件为Xilinx ISE 11.4，FPGA型号为 Virtex－5vsx50，Xilinx Virtex－5 系列［19］的 FPGA最大时钟频率是200 MHz，然而考虑到总线接口和外部主机的交互，文中采用了100 MHz的参考时钟作为全局时钟，即完成PNN的训练需要25μs，FPGA计算一个锋电位分类需要6.55 μs;同时在Matlab 7.8.0.347上实现了PNN算法完成锋电位信号的分类.计算机配置为Intel®coreTMDuo CPU E6550@2.3GHz和 2.0GB DDR2－800.

微电极采集获得的36通道的神经元信号数据的分类结果如图8所示.从图8可知:基于FPGA的PNN进行锋电位分类的时间是固定不变的，平均分类准确率和基于Matlab的PNN是一致的;基于Matlab的PNN完成锋电位分类的平均时间(36通道数据，每通道执行1000次的平均计算时间)为317.22μs，36通道实现锋电位分类的平均准确率为93.82%.

图8 基于Matlab和FPGA的PNN的分类结果Fig.8 Sorting results of PNNs respectively based on Matlab and FPGA

由Cerebus 128通道数据采集系统、基于Matlab的PNN和基于FPGA的PNN完成第12通道的锋电位的分类结果如图9所示.基于Matlab和基于FPGA的PNN的分类准确率一样，相对标准的Cerebus只有6.55%的差别.

基于Matlab和基于FPGA的PNN的平均分类准确率和平均分类时间如表1所示.从表1可知，在不同的实验平台上，同样的数据输入位宽，分类准确率相同，说明FPGA中激活函数的逼近计算并没有影响分类精度，但基于FPGA的PNN的分类时间明显少于基于Matlab的PNN.

图9 3种方法的分类结果比较Fig.9 Comparison of sorting results of three methods

表1 基于Matlab和FPGA的PNN的平均分类准确率和分类时间Table 1 Average sorting accuracy and sorting time of PNNs respectively based on Matlab and FPGA

从实验结果可以看出，基于FPGA的PNN的平均分类准确率为93.82%，训练时间为25μs，分类时间为6.55μs，足以完成锋电位的实时在线分类.在相同的分类准确率下，基于FPGA的PNN由于采用了并行计算结构，因而在分类速度上比基于Matlab的PNN快47.43倍.FPGA实现PNN模型时的资源使用情况如表2所示.

表2 FPGA的资源使用情况Table 2 Resource utilization of FPGA

从表2可知，文中算法的FPGA实现只占用了较少的资源，只有Slice LUTs和DSP48Es占用了50%左右的资源，这是因为它们采用了浮点精度的计算.

4 结语

为了实现锋电位信号的实时分类，文中提出了一种新型的基于FPGA的概率神经网络算法.实验结果表明，该算法的平均分类准确率可以高达93.82%，其分类速度比基于Matlab的PNN快47.43倍.今后将重点研究与前端采集装置的实时连接，同时对FPGA资源的利用率进行改进，以期为实现便携式实时处理的脑－机接口设备提供实验基础.

［1］ Sanchez Justin C，Principe Jose C.Brain－machine interface engineering［M］.San Francisco:Morgan＆ Claypool Publishers，2007:1－20.

［2］ Lebedev M A，Nicolelis M A L.Brain－machine interface:past，present and future ［J］.Trends in Neurosciences，2006，29(9):536－546.

［3］ Ye X S，Wang P，Liu J，et al.A portable telemetry system for brain stimulation and neuronal activity recording in freely behaving small animals ［J］.Journal of Neuroscience Methods，2008，174(2):186－193.

［4］ 代建华，章怀坚，张韶岷，等.大鼠运动皮层神经元集群锋电位时空模式分析［J］.中国科学C辑:生命科学，2009，39(8):736－745.Dai Jian－hua，Zhang Huai－jian，Zhang Shao－min，et al.Spatiotemporal analysis of the rat motor cortex neurons clustered spike［J］.Science in China Series C:Life Sciences，2009，39(8):736－745.

［5］ Nicolelis M A L.Brain－machine interfaces to restore motor function and probe neural circuits［J］.Nature Reviews Neuroscience，2003，4(5):417－422.

［6］ 祝晓平，王东，陈耀武.基于提升小波的神经元锋电位并行检测方法［J］.华南理工大学学报:自然科学版，2011，39(10):19－25，31.Zhu Xiao－ping Wang Dong，Chen Yao－wu.Parallel detection method of neuronal spikes based on lifting wavelet［J］.Journal of South China University of Technology:Natural Science Edition，2011，39(10):19－25，31.

［7］ Yang Z，Zhao Q，Liu W T.Neural signal classification using a simplified feature set with nonparametric clustering［J］.Neurocomputing，2009，73(1/2/3):412－422.

［8］ Zumsteg Z S，Kemere C，Driscoll S，et al.Power feasibility of implantable digital spike sorting circuits for neural prosthetic systems［J］.IEEE Transactions on Neural Systems and Rehabilitation Engineering，2005，13(3):272－279.

［9］ Dai J H，Liu X C，Yu Y，et al.Experimental study on neural spike sorting methods［C］∥Proceedings of the Second International Conference on Future Generation Communication and Networking.Hainan:IEEE，2008:230－233.

［10］ Vollgraf R，Obermayer K.Improved optimal linear filters for the discrimination of multichannel waveform templates for spike－sorting applications ［J］.IEEE Signal Processing Letters，2006，13(3):121－124.

［11］ Hulata E，Segev R，Ben－Jacob E.A method for spike sorting and detection based on wavelet packets and Shannon’s mutual information ［J］.Journal of Neuroscience Methods，2002，117(1):1－12.

［12］ Specht D F.Probabilistic neural networks ［J］.Neural Networks，1990，3(1):109－118.

［13］ Zhou F，Liu J，Yu Y，et al.Field－programmable gate array implementation of a probabilistic neural network for motor cortical decoding in rats［J］.Journal of Neuroscience Methods，2010，185(2):299－306.

［14］ Vinitha K V，Santhosh Kumar G.Face recognition using probabilistic neural networks［C］∥Proceedings of World Congress on Nature and Biologically Inspired Computing.Coimbatore:IEEE，2009:1388－1393.

［15］ 王静，封洲燕.多通道神经元锋电位检测和分类的新方法［J］.生物化学与生物物理进展，2009，36(5):641－647.Wang Jing，Feng Zhou－yan.A novel method for multichannel neuronal spike detection and classification［J］.Progress in Biochemistry and Biophysics，2009，36(5):641－647.

［16］ Xilinx.LogiCORE IP floating－pointer operator v5.0 data sheet［EB/OL］.(2011－03－01)［2011－09－10］.http:∥www.xilinx.com/support/documentation/ipdsp_arithmetic_floating－pt.htm.

［17］ Volder Jack E.The CORDIC trigonometric computing technique［J］.IRE Transactions on Electronic Computers，1959，8(3):330－334.

［18］ Xilinx.CORDIC v4.0 data sheet［EB/OL］.(2011－03－01)［2011－09－10］.http:∥www.xilinx.com/support/documentation/ipdsp_trigfunction_cordic.htm.

［19］ Xilinx.Virtex－5 FPGA family data sheet ［EB/OL］.(2009－02－06)［2011－09－10］.http:∥www.xilinx.com/products/virtex5/index.htm.