CN109376651A - 一种基于cuda架构的gpu加速锋电位分类的系统 - Google Patents

Abstract

本发明提供一种基于CUDA架构的GPU加速锋电位分类系统，包括流式数据输入模块、滤波模块、锋电位提取模块、锋电位特征提取模块、锋电位聚类模块和锋电位增量聚类模块；本发明首先系统具有很好的抗噪性，可以提高分类的准确度，提高分类的速度，其次系统的全部模块基于CUDA架构实现，使得运行时间大幅缩短，在高维信道上可以达到实时分类，最后系统的模块化设计，具有良好的重用性和可扩展性，该系统极大的方便了脑神经科学方面的研究。

Description

一种基于CUDA架构的GPU加速锋电位分类的系统

技术领域

本发明涉及机器学习技术领域，尤其是一种基于CUDA架构的GPU加速锋电位分类系统。

背景技术

锋电位分类是用于分析电生理数据的一类技术，该技术是根据锋电位形状的相似性将锋电位分组为聚类。原则上，每个神经元倾向于激发特定形状的锋电位，所得到的聚类对应于不同的假定神经元的活动，锋电位分类的最终结果是确定哪个锋电位与这些神经元中的哪一个相对应，其中，CUDA架构(Compute Unified Device Architecture)是面向GPU(Graphic Processing Unit)图形处理器等设备的并行计算架构，是一套在GPU上完成高性能计算的解决方案，在CUDA环境上并行编程的接口API有CUDA C，C++，OpenCL等。CUDA内核程序被称为kernel函数，是在GPU上运行的并行计算函数，同时CUDA也包含很多高性能的函数库，如CUBLAS、CUSPARSE、CUSOLVER和THRUST等函数库，使开发者方便进行高性能并行程序的开发。

目前，锋电位分类系统主要有klusta系统、kilosort系统和mountainsort系统等。klusta系统聚类应用Masked EM算法，该算法基于高斯混合模型，抗噪性能较好，虽然达到了比较高的分类精度，但是该系统并没有进行GPU加速，系统运行时间无法接受，另外该系统无法支持在线数据的处理；kilosort系统采用了模板匹配算法，同时实现了在GPU上的加速，可以实现锋电位信号的实时处理，但是由于锋电位信号的不稳定性，该方法可能受到噪声信号的影响；mountainsort系统聚类采用ISO-SPLIT算法，在维度较低的信道上计算效率较高，但是在高维信道仍然无法达到实时处理。

发明内容

针对现有技术的不足，本发明提供一种基于CUDA架构的GPU加速锋电位分类系统。

本发明的技术方案为：一种基于CUDA架构的GPU加速锋电位分类系统，包括流式数据输入模块、滤波模块、锋电位提取模块、锋电位特征提取模块、锋电位聚类模块和锋电位增量聚类模块；

所述的流式数据输入模块用于获取脑电波信号采集设备采集的数据，并将数据传入 GPU内存中，并进行数据分块；

所述的滤波模块用于将分块的数据采用GPU加速的滤波器进行带通滤波，滤除背景噪声和局部场电位；

所述的锋电位提取模块用于将滤波后的数据按照每个信道计算得到一个高阈值和一个低阈值，根据阈值将波形块从原始数据中提取出来，并对提取出来的波形块平滑化和尖锋对齐；

所述的锋电位特征提取模块用于将提取到的每个波形块的每个信道进行特征分解，并选取每个信道前三维特征表示该信道的主要特征，用于之后的聚类工作；

所述的锋电位聚类模块用于根据提取的一定量初始波形块的特征，进行迭代聚类，并计算各个聚类的参数；

所述的锋电位增量聚类模块用于对之后的流式增量数据进行增量聚类，更新聚类参数，并将聚类结果传回到CPU端。

进一步的，所述的流式数据输入模块采用window-based的划分方式，window的时长为 1s，同时支持在线和离线数据输入。

进一步的，所述的滤波模块采用在GPU上加速实现的三阶巴特沃斯带通滤波器，低频为600Hz，高频为0.95*Nyquist。

进一步的，所述的锋电位提取模块利用GPU计算出每个信道的电位信号中位数V_median，那么阈值基数Θ＝V_median/0.6745，高阈值为4.5Θ，低阈值为2Θ，采用Flood fill算法得到波形块为32×32，波形块平滑化采用自然边界的三次样条插值算法。

进一步的，所述的锋电位特征提取模块对波形块特征分解采用GPU加速的SVD分解，每个信道取前三维特征表示该信道主要特征，将每个信道的主要特征拼接起组成96维的波形块特征。

进一步的，所述的锋电位聚类模首先采用Masked EM算法进行预聚类，其中预聚类数据时长为20s，然后对于增量数据采用增量Masked EM聚类，通过预聚类结果和增量锋电位数据特征进行快速聚类和参数更新，该算法可以使后来的数据可以利用之前数据的聚类信息进行增量计算，而不需要重复计算或者抛弃之前的计算信息，不仅可以增加计算的准确性，同时可以更好的处理流式数据，更加符合脑电波实时分析的场景要求。

进一步的，所述的Masked EM算法通过数据的特征根据重要程度计算出一个masked向量，从而掩蔽掉大部分的没有明显特征的维度，同时对有特征的维度进行加权，使得高维特征数据的聚类的计算量大大减少，提高聚类效率，算法的具体描述如下：

1、对原始数据点的特征生成掩蔽向量(mask向量)，即

其中，m_n，i为x_n，i的mask向量，x_n，i表示第n个数据点的第i个特征，sD_i表示第i个特征的阙值，a，β分别为阈值的上下边界系数；

2、通过掩蔽向量计算被masked的每个特征i的噪声的均值v_i和方差具体如下：

其中，表示特征i被掩蔽的样本数，表示所有数据点的特征i中 m_n，i＝0的特征的数目；

3、将原始数据点通过高斯混合的最大似然估计拟合为一个虚拟的分布，该分布可以表示为：

其中，为拟合之后的虚拟分布，表示均值为v_i，方差为的正态分布，prob m_m，i表示虚拟分布有m_n，i的概率被表示为x_n，i，同理，prob 1-m_n，i表示虚拟分布有 1-m_n，i的概率被表示为参数为(v_i，σ_i ²)的正态分布；

通过噪声分布的均值和方差，可以将原始数据拟合为输入数据，

y_n，i＝m_n，ix_n，i+(1-m_n，i)v_i

η_n，i＝z_n，i-(y_n，i)²；

其中，y_n，i表示被masked向量拟合后的输入数据的每个特征的期望，η_n，i表示被masked 向量拟合后的输入数据的每个特征的方差，z_n，i表示虚拟分布的第n个数据点的第i个特征的平方的期望，用于计算拟合后输入数据的方差，v_i为特征i的噪声的均值，σ_i表示所有数据点的特征i的噪声的标准差；

4、进行Mstep和Estep，通过汉明距离计算两个样本点掩蔽向量的距离，并将距离相近的样本点初始化为一个簇，计算簇的权重ω_k、均值μ_k与协方差∑_k，更新样本点对于每个簇的似然对数值，在Masked EM模型中，对于簇k更新权重、特征均值和协方差的公式如下：

其中，ω_k表示第k簇的权重，C_k表示第k簇的数据点的集合，N为数据点的总个数，(μ_k)_i表示第k个簇的第i个特征的均值，(μ_k)_j表示第k个簇的第j个特征的均值，M_k，i表示第k个簇中第i个特征被掩蔽的数据点的集合，(∑_k)_i，j表示第k个簇的第i个特征和第j个特征的协方差，δ_i，j为第i个特征和第j个特征的相关系数，表示簇k中特征i未被掩蔽的所有数据点的方差的和；π_n，k表示第n个数据点属于第i个簇的对数似然值，D为数据点的特征维度，T为矩阵或向量转置符号，det∑_k表示第k个簇的协方差矩阵∑_k的行列式的值，为对角线校正项，y_n为拟合后的输入数据，μ_k为簇k的均值，为第k 个簇的第i个特征的方差的逆；

5、根据得到的对数似然值π_n，k进行簇的重新分配，并通过计算簇的得分考虑删除和分裂簇的操作；

6、迭代第4，5步，直到满足一定的条件，退出迭代。

进一步的，对于增量数据采用锋电位增量聚类模块进行增量聚类，通过预聚类结果和增量锋电位数据特征进行快速聚类和参数更新，具体为：

通过Masked EM算法得到的聚类模型，对于第t+1次的增量数据样本点y_n+1，可以通过E步计算得到 表示t+1次的增量数据y_n+1的对数似然估计，假设样本点y_n+1的类别为K^*，那么有K^*＝arg max{π_n+1，i}，i＝1，2，....，K，其中，K为当前的簇的总个数，K^*为对数似然估计最大的簇，除了与样本点最相似的簇K^*，由于其余簇的估计参数变化都很轻微，故可近似为加入样本点y_n+1之前的参数，对于第K^*个簇的参数估计，通过上述的MStep和Estep更新和从而完成对样本点y_n+1的分类和参数更新。

本发明的有益效果为：首先系统具有很好的抗噪性，可以提高分类的准确度，提高分类的速度，其次系统的全部模块基于CUDA架构实现，使得运行时间大幅缩短，在高维信道上可以达到实时分类，最后系统的模块化设计，具有良好的重用性和可扩展性，该系统极大的方便了脑神经科学方面的研究。

附图说明

图1为本发明的结构框架图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明：

如图1所示，一种基于CUDA架构的GPU加速锋电位分类系统，包括流式数据输入模块、滤波模块、锋电位提取模块、锋电位特征提取模块、锋电位聚类模块和锋电位增量聚类模块；

所述的流式数据输入模块用于获取脑电波信号采集设备采集的数据，并将数据传入 GPU内存中，并进行数据分块，其中，所述的流式数据输入模块采用window-based的划分方式，window的时长为1s，同时支持在线和离线数据输入；

所述的滤波模块用于将方块的数据采用GPU加速的滤波器进行带通滤波，滤除背景噪声和局部场电位，其中，所述的滤波模块采用在GPU上加速实现的三阶巴特沃斯带通滤波器，低频为600Hz，高频为0.95*Nyquist；

所述的锋电位提取模块用于将滤波后的数据按照每个信道计算得到一个高阈值和一个低阈值，根据阈值将波形块从原始数据中提取出来，并对提取出来的波形块平滑化和尖锋对齐，其中，所述的锋电位提取模块利用GPU计算出每个信道的电位信号中位数V_median，那么阈值基数Θ＝V_median/0.6745，高阈值为4.5Θ，低阈值为2Θ，采用Flood fill算法得到波形块为32×32，波形块平滑化采用自然边界的三次样条插值算法。

所述的锋电位特征提取模块用于将提取到的每个波形块的每个信道进行特征分解，并选取每个信道前三维特征表示该信道的主要特征，用于之后的聚类工作，其中，所述的锋电位特征提取模块对波形块特征分解采用GPU加速的SVD分解，每个信道取前三维特征表示该信道主要特征，将每个信道的主要特征拼接起组成96维的波形块特征。

所述的锋电位聚类模块用于根据提取的一定量初始波形块的特征，进行迭代聚类，并计算各个聚类的参数，其中，所述的锋电位聚类模首先采用Masked EM算法进行预聚类，其中预聚类数据时长为20s，然后对于增量数据采用锋电位增量聚类模块进行增量聚类，通过预聚类结果和增量锋电位数据特征进行快速聚类和参数更新，该算法可以使后来的数据可以利用之前数据的聚类信息进行增量计算，而不需要重复计算或者抛弃之前的计算信息，不仅可以增加计算的准确性，同时可以更好的处理流式数据，更加符合脑电波实时分析的场景要求，

所述的Masked EM算法通过数据的特征根据重要程度计算出一个masked向量，从而掩蔽掉大部分的没有明显特征的维度，同时对有特征的维度进行加权，使得高维特征数据的聚类的计算量大大减少，提高聚类效率，算法的具体描述如下：

1、对原始数据点的特征生成掩蔽向量(mask向量)，即

其中，m_n，i为x_n，i的mask向量，x_n，i表示第n个数据点的第i个特征，sD_i表示第i个特征的标准差，a，β分别为阈值的上下边界系数；

2、通过掩蔽向量计算被masked的每个特征i的噪声的均值v_i和方差具体如下：

其中，表示特征i被掩蔽的样本数，表示所有数据点的特征i中 m_n，i＝0的特征的数目；

3、将原始数据点通过高斯混合的最大似然估计拟合为一个虚拟的分布，该分布可以表示为：

其中，为拟合之后的虚拟分布，表示均值为v_i，方差为的正态分布，prob m_m，i表示虚拟分布有m_n，i的概率被表示为x_n，i，同理，prob 1-m_n，i表示虚拟分布有 1-m_n，i的概率被表示为参数为(v_i，σ_i ²)的正态分布；

通过噪声分布的均值和方差，可以将原始数据拟合为输入数据，

y_n，i＝m_n，ix_n，i+(1-m_n，i)v_i

η_n，i＝z_n，i-(y_n，i)²；

其中，y_n，i表示被masked向量拟合后的输入数据的每个特征的期望，η_n，i表示被masked 向量拟合后的输入数据的每个特征的方差，z_n，i表示虚拟分布的第n个数据点的第i个特征的平方的期望，用于计算拟合后输入数据的方差，v_i为特征i的噪声的均值，σ_i表示所有数据点的特征i的噪声的标准差；

4、进行Mstep和Estep，通过汉明距离计算两个样本点掩蔽向量的距离，并将距离相近的样本点初始化为一个簇，计算簇的权重ω_k、均值μ_k与协方差∑_k，更新样本点对于每个簇的似然对数值，在Masked EM模型中，对于簇k更新权重、特征均值和协方差的公式如下：

其中，ω_k表示第k簇的权重，C_k表示第k簇的数据点的集合，N为数据点的总个数，(μ_k)_i表示第k个簇的第i个特征的均值，(μ_k)_j表示第k个簇的第j个特征的均值，M_k，i表示第k个簇中第i个特征被掩蔽的数据点的集合，(∑_k)_i，j表示第k个簇的第i个特征和第j个特征的协方差，δ_i，j为第i个特征和第j个特征的相关系数，表示簇k中特征i未被掩蔽的所有数据点的方差的和；π_n，k表示第n个数据点属于第k个簇的对数似然值，D为数据点的特征维度，T为矩阵或向量转置符号，det∑_k表示第k个簇的协方差矩阵∑_k的行列式的值，为对角线校正项，y_n为拟合后的输入数据，μ_k为簇k的均值，为第k 个簇的第i个特征的方差的逆；

5、根据得到的对数似然值π_n，k进行簇的重新分配，并通过计算簇的得分考虑删除和分裂簇的操作，具体为：通过π_n，k将每个样本点重新分配到似然函数最大的簇中，对于重新分配后的所有簇，其得分对于每个样本点，采用并行规约算法找到其最大对数似然的所属簇，更新所属簇，并通过并行规约算法求得更新S；

簇的删除和分裂，在初始时给定一个可能的最大的簇数量，之后在Mstep和Estep迭代过程中，通过判断当前簇是否满足分裂或删除的条件，动态的调整簇的数量；

簇的删除判断如下：对每个簇i，将其中的样本点重新分配到次大的似然函数值所对应的簇中，并更新所有簇中的每个样本点的似然函数的值和每个簇的惩罚项的值，更新之后重新计算得分记为S_i，对于每个簇都尝试删除，可以找到最小的S_new＝min(S₁，S₂，....，S_k)，如果 S_new＜S，则说明删除之后结果更优，所以对得分最小的簇进行删除，将该簇中的样本点重新分配；

簇分裂判断如下：针对每个簇，通过随机初始化为2个簇，并通过M步和E步进行更新，最终将分成的两个簇，一个簇为当前簇号，另一个簇新分配一个簇号，并通过M步和 E步更新之后，重新计算所有簇的得分，如果当前得分小于未分裂得分S，则表示当前簇分裂更优，所以将当前簇分裂为两个；

6、迭代第4，5步，直到满足一定的条件，退出迭代。

对于增量数据采用锋电位增量聚类模块进行增量聚类，通过预聚类结果和增量锋电位数据特征进行快速聚类和参数更新，具体为：

通过Masked EM算法得到的聚类模型，对于第t+1次的增量数据样本点y_n+1，可以通过 E步计算得到 表示t+1次的增量数据y_n+1的对数似然估计，假设样本点y_n+1的类别为K^*，那么有K^*＝arg max{π_n+1，i}，i＝1，2，....，K，其中，K为当前的簇的总个数，K^*为对数似然估计最大的簇，除了与样本点最相似的簇K^*，由于其余簇的估计参数变化都很轻微，故可近似为加入样本点y_n+1之前的参数，对于第K^*个簇的参数估计，通过上述的MStep和Estep 更新和从而完成对样本点y_n+1的分类和参数更新。

所述的锋电位增量聚类模块用于对之后的流式增量数据进行增量聚类，更新聚类参数，并将聚类结果传回到CPU端。

上述实施例和说明书中描述的只是说明本发明的原理和最佳实施例，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。

Claims (8)

1.一种基于CUDA架构的GPU加速锋电位分类系统，其特征在于：包括流式数据输入模块、滤波模块、锋电位提取模块、锋电位特征提取模块、锋电位聚类模块和锋电位增量聚类模块；

所述的流式数据输入模块用于获取脑电波信号采集设备采集的数据，并将数据传入

GPU内存中，并进行数据分块；

所述的滤波模块用于将分块的数据采用GPU加速的滤波器进行带通滤波，滤除背景噪声和局部场电位；

所述的锋电位提取模块用于将滤波后的数据按照每个信道计算得到一个高阈值和一个低阈值，根据阈值将波形块从原始数据中提取出来，并对提取出来的波形块平滑化和尖锋对齐；

所述的锋电位特征提取模块用于将提取到的每个波形块的每个信道进行特征分解，并选取每个信道前三维特征表示该信道的主要特征，用于之后的聚类工作；

所述的锋电位聚类模块用于根据提取的一定量初始波形块的特征，进行迭代聚类，并计算各个聚类的参数；

所述的锋电位增量聚类模块用于对之后的流式增量数据进行增量聚类，更新聚类参数，并将聚类结果传回到CPU端。

2.根据权利要求1所述的一种基于CUDA架构的GPU加速锋电位分类系统，其特征在于：所述的流式数据输入模块采用window-based的划分方式，window的时长为1s，同时支持在线和离线数据输入。

3.根据权利要求1所述的一种基于CUDA架构的GPU加速锋电位分类系统，其特征在于：所述的滤波模块采用在GPU上加速实现的三阶巴特沃斯带通滤波器，低频为600Hz，高频为0.95*Nyquist。

4.根据权利要求1所述的一种基于CUDA架构的GPU加速锋电位分类系统，其特征在于：所述的锋电位提取模块利用GPU计算出每个信道的电位信号中位数V_median，那么阈值基数Θ＝V_median/0.6745，高阈值为4.5Θ，低阈值为2Θ，采用Flood fill算法得到波形块为32×32，波形块平滑化采用自然边界的三次样条插值算法。

5.根据权利要求1所述的一种基于CUDA架构的GPU加速锋电位分类系统，其特征在于：所述的锋电位特征提取模块对波形块特征分解采用GPU加速的SVD分解，每个信道取前三维特征表示该信道主要特征，将每个信道的主要特征拼接起组成96维的波形块特征。

6.根据权利要求1所述的一种基于CUDA架构的GPU加速锋电位分类系统，其特征在于：所述的锋电位聚类模首先采用Masked EM算法进行预聚类，其中预聚类数据时长为20s，然后对于增量数据采用锋电位增量聚类模块进行增量聚类，通过预聚类结果和增量锋电位数据特征进行快速聚类和参数更新。

7.根据权利要求6所述的一种基于CUDA架构的GPU加速锋电位分类系统，其特征在于：所述的Masked EM算法通过数据的特征根据重要程度计算出一个masked向量，从而掩蔽掉大部分的没有明显特征的维度，同时对有特征的维度进行加权，使得高维特征数据的聚类的计算量大大减少，提高聚类效率，算法的具体描述如下：

1)、对原始数据点的特征生成掩蔽向量(mask向量)，即

其中，m_n,i为x_n，i的mask向量，x_n,i表示第n个数据点的第i个特征，sD_i表示第i个特征的阙值，a，β分别为阈值的上下边界系数；

2)、通过掩蔽向量计算被masked的每个特征i的噪声的均值v_i和方差具体如下：

其中，表示特征i被掩蔽的样本数，表示所有数据点的特征i中m_n，i＝0的特征的数目；

3)、将原始数据点通过高斯混合的最大似然估计拟合为一个虚拟的分布，该分布可以表示为：

其中，为拟合之后的虚拟分布，表示均值为v_i，方差为的正态分布，probm_m,i表示虚拟分布有m_n,i的概率被表示为x_n,i，同理，prob1-m_n，i表示虚拟分布有1-m_n,i的概率被表示为参数为(v_i,σ_i ²)的正态分布；

通过噪声分布的均值和方差，可以将原始数据拟合为输入数据，

y_n,i＝m_n,ix_n,i+(1-m_n,i)v_i

η_n,i＝z_n,i-(y_n,i)²；

其中，y_n,i表示被masked向量拟合后的输入数据的每个特征的期望，η_n,i表示被masked向量拟合后的输入数据的每个特征的方差，z_n,i表示虚拟分布的第n个数据点的第i个特征的平方的期望，用于计算拟合后输入数据的方差，v_i为特征i的噪声的均值，σ_i表示所有数据点的特征i的噪声的标准差；

4)、进行Mstep和Estep，通过汉明距离计算两个样本点掩蔽向量的距离，并将距离相近的样本点初始化为一个簇，计算簇的权重ω_k、均值μ_k与协方差∑_k，更新样本点对于每个簇的似然对数值，在Masked EM模型中，对于簇k更新权重、特征均值和协方差的公式如下：

其中，ω_k表示第k簇的权重，C_k表示第k簇的数据点的集合，N为数据点的总个数，(μ_k)_i表示第k个簇的第i个特征的均值，(μ_k)_j表示第k个簇的第j个特征的均值，M_k，i表示第k个簇中第i个特征被掩蔽的数据点的集合，(∑_k)_i，j表示第k个簇的第i个特征和第j个特征的协方差，δ_i,j为第i个特征和第j个特征的相关系数，表示簇k中特征i未被掩蔽的所有数据点的方差的和；π_n,k表示第n个数据点属于第i个簇的对数似然值，D为数据点的特征维度，T为矩阵或向量转置符号，det∑_k表示第k个簇的协方差矩阵∑_k的行列式的值，为对角线校正项，y_n为拟合后的输入数据，μ_k为簇k的均值，为第k个簇的第i个特征的方差的逆；

5)、根据得到的对数似然值π_n,k进行簇的重新分配，并通过计算簇的得分考虑删除和分裂簇的操作；

6)、迭代第4，5步，直到满足一定的条件，退出迭代。

8.根据权利要求6所述的一种基于CUDA架构的GPU加速锋电位分类系统，其特征在于：对于增量数据采用锋电位增量聚类模块进行增量聚类，通过预聚类结果和增量锋电位数据特征进行快速聚类和参数更新，具体为：

通过Masked EM算法得到的聚类模型，对于第t+1次的增量数据样本点y_n+1，可以通过E步计算得到 表示t+1次的增量数据y_n+1的对数似然估计，假设样本点y_n+1的类别为K^*，那么有K^*＝argmax{π_n+1,i}，i＝1,2，....,K,其中，K为当前的簇的总个数，K^*为对数似然估计最大的簇，除了与样本点最相似的簇K^*，由于其余簇的估计参数变化都很轻微，故可近似为加入样本点y_n+1之前的参数，对于第K^*个簇的参数估计，通过上述的MStep和Estep更新和从而完成对样本点y_n+1的分类和参数更新。