CN105718425A - 一种非线性数据全局互相关性的并行定量计算方法 - Google Patents

Abstract

本发明公开了一种非线性数据全局互相关性的并行定量计算方法，采用高性能并行计算中的CPU多线程技术和基于CUDA的GPGPU多线程方法开发并行NLI算法，将并行计算思想融入到算法中去，改进并提高算法的性能，辅助多通道信号的分析，将数据的处理从双通道扩展到多通道，并结合S估计器将数据的分析从局部到整体再到局部这样一个过程，量化信号一定区域范围内的同步强度，从原始数据中挖掘出更多更准确更有用的信息来分析多通道信号的同步问题，在保留其算法效果的同时大大提高其执行效率。实验证明，本方法在实际的多通道非线性数据互相关性分析中具有更好的效率和可用性。

Description

一种非线性数据全局互相关性的并行定量计算方法

技术领域

本发明属于数据分析技术领域，涉及一种数据相关性分析方法，尤其涉及一种同时揭示相关强度与相关方向的、全局的、非线性数据全局互相关性的并行定量计算方法。

背景技术

科学与实际工程应用中，复杂系统的监测数据(如大脑神经信号)具有明显的非线性特性，复杂系统多个组成部分之间同步现象的研究是当今复杂系统科学领域中最活跃的课题之一。以神经信号为例，同步测量是对大脑活动时产生的两个或更多连续时间序列的同步估计，独立的时间序列同步值低，相关的时间序列同步值高。利用多通道脑电信号分析方法研究神经信号的同步特征有助于深入理解生理病理情况下大脑各区域间的协作和调制机制，为脑部神经疾病的诊断提供依，从而进行正确的治疗，通过研究理解导致脑部神经疾病的机理，也有助于对疾病的提前预防。([文献1])

现有的一些二元信号同步性分析方法，比如交叉相关性([文献1])，频谱相干性([文献2])，相位同步性([文献3])和非线性互相关性(Non-linearInterdependency，NLI)([文献4])等方法。非线性互相关性(NLI)是一种非对称性的相关性测量方法，用来描述信号的相位同步强度，可用于衡量广义同步。NLI算法不进能够计算成对二元信号的强度，还能计算二元非线性数据同步的方向([文献5])，并且具有很好的抗噪声干扰能力。

而多通道信号比二元信号所包含的信息更多。为了进一步分析多通道信号之间的同步性关系，(如，为了探索不同大脑区域的工作机理我们需要分析其不同通道的神经信号的相关性)，多通道信号分析需要一个计算全局相关性的方法。因此，相比较于二元信号相关性分析，由于多通道信号的普适性，多通道信号分析收到了越来越多的重视。现有的一些多通道信号互相关性分析方法，S估计器([文献6])，关联性矩阵分析([文献7])，互信息，局部定向相关性([文献8])。文献9对这些方法做了一个比较，非线性互相关性方法(NLI)对噪声的鲁棒性明显强于其它的方法。因此，相比较而言，NLI算法更适合于多通道数据全局互相关性的定量计算。

然而，在多通道信号分析上应用NLI算法是高度计算密集型的，NLI算法主要应用于研究两通道间信号关系，虽然该算法在处理脑电信号上有其独特的优势，但当数据量较大或通道数较多时，该算法的执行效率明显下降，从而制约了该算法的发展和使用。因此亟待发明一种快速有效的多通道非线性数据相关性分析算法。

[文献1]C.Carmeli，M.G.Knyazeva，G.M.Innocenti，andO.D.Feo，“AssessmentofEEGsynchronizationbasedonstate-spaceanalysis，”Neuroimage，vol.25，no.2，pp.339–354，2005.

[文献2]A.Kraskov，“Synchronizationandinterdependencemeasuresandtheirapplicationstotheelectroencephalogramofepilepsypatientsandclusteringofdata，NIC-Directors，”Ph.D.，Jülich，Germany，2004.

[文献3]F.Mormann，K.Lehnertz，P.David，andC.E.Elger，“MeanphasecoherenceasameasureforphasesynchronizationanditsapplicationtotheEEGofepilepsypatients，”PhysicaD，vol.144，no.3–4，pp.358–369，2000

[文献4]M.LeVanQuyen，J.Soss，V.Navarro，R.Robertson，M.Chavez，M.Baulac，andJ.Martinerie，“PreictalstateidentificationbysynchronizationchangesinlongtermintracranialEEGrecordings，”Clin.Neurophysiol.，vol.116，no.3，pp.559–568，2005.

[文献5]M.Breakspear，J.R.Terry，K.J.Friston，A.W.F.Harris，L.M.Williams，K.Brown，J.Brennan，andE.Gordon，“AdisturbanceofnonlinearinterdependenceinscalpEEGofsubjectswithfirstepisodeschizophrenia，”NeuroImage，vol.20，no.1，pp.466–478，2003.

[文献6]M.G.Knyazeva，G.M.Innocenti，C.Carmeli，andO.D.Feo，“AssessmentofEEGsynchronizationbasedonstate-spaceanalysis，”Neuroimage，vol.25，no.2，pp.339–354，2005.

[文献7]X.Li，D.Cui，P.Jiruska，J.E.Fox，X.Yao，andJ.G.R.Jefferys，“Synchronizationmeasurementofmultipleneuronalpopulations，”J.Neurophysiol.，vol.98，pp.3341–3348，2007.

[文献8]K.SameshimaandL.A.Baccala，“Usingpartialdirectedcoherencetodescribeneuronalensembleinteractions，”J.Neurosci.Methods，vol.94，pp.93–103，1999.

[文献9]T.Kreuz，F.Mormann，R.Andrzejak，A.Kraskov，K.Lehnertz，andP.Grassberger，“Measuringsynchronizationincoupledmodelsystems:Acomparisonofdifferentapproaches，”PhysicaD，vol.225，no.1，pp.29–42，2007.

发明内容

为了解决上述技术问题，本发明提供了一种同时揭示相关强度与相关方向的、全局的、非线性数据全局互相关性的并行定量计算方法。

本发明所采用的技术方案是：一种非线性数据全局互相关性的并行定量计算方法，其特征在于，包括以下步骤：

步骤1：对M通道原始脑电数据shuru进行直接分解，M＞2，以两两通道数据最为原始数据进入下层循环；设定初始通道i＝1，j＝1；

步骤2：取第i通道原始数据shurui；

步骤3：取第j通道原始数据shuruj；

步骤4：设置数据时间窗Epoch和滑动距离Overlap；

步骤5：设定初始循环次数为t＝1，计算循环次数N：

N＝length(shuru-Epoch)/(Epoch-Overlap)

步骤4：按照时间窗从数据输入shurui和shuruj里取数据，计算当前窗口数据偏移量(t-1)*(Epoch-Overlap)，得到当前数据X_t(k)和Y_t(k)；

步骤5：运用NLI算法计算两个数据的相关性，得到独立性测量S、H、N、M；

步骤6：判断t＞N是否成立；

若是，则顺序执行下述步骤6；

若否，则t＝t+1，并回转执行所述步骤4；

步骤7：判断j＞M是否成立；

若是，则顺序执行下述步骤8；

若否，则j＝j+1，并回转执行所述步骤3；

步骤8：判断i＞M是否成立；

若是，则顺序执行下述步骤9；

若否，则i＝i+1，并回转执行所述步骤2；

步骤9：对独立性测量S、H、N、M组成的相关性矩阵进行特征值特征向量分解；

对数据矩阵X_M×T＝{x_i(k)}，i＝1,…,M，k＝1,…,T的协方差矩阵C进行特征值分解，其中M表示通道数量，T表示时间窗内的数据点数，获得特征值λ_i，并对其进行归一化，得到归一化特征值：

${{\mathrm{\λ}}_i}^{\′}=\frac{{\mathrm{\λ}}_i}{\underset{i=1}{\overset{M}{\Σ}}{\mathrm{\λ}}_i}=\frac{{\mathrm{\λ}}_i}{tr\left(C\right)};$

式中tr(C)为协方差矩阵的迹；

步骤10：基于特征值特征向量计算S估计器，S估计器的定义为：

$S=1+\frac{\underset{i=1}{\overset{M}{\Σ}}{{\mathrm{\λ}}_i}^{\′}log\left({{\mathrm{\λ}}_i}^{\′}\right)}{1og\left(M\right)}\∈\[0,1\].$

作为优选，步骤5中所述NLI算法，其具体实现过程包括以下子步骤：

步骤5.1：给定两个输入的信号数据X和Y，重构两个输入信号的延迟向量x_n＝(x_n,…,x_n-(m-1)τ)和y_n＝(y_n,…,y_n-(m-1)τ)，其中n＝1,…,N，m为嵌入维数，τ为预定义的延迟时间；N＝T-(m-1)τ，1≤n≤N；所有延迟向量组成了新的矩阵X＝(x₁,x₂,…,x_N)，Y＝(y₁,y₂,…,y_N)；

步骤5.2：分别计算X，Y中两两延迟向量之间的距离distXX,distXY，即distxx_n＝||x_n-x_j||，distxy_n＝||y_n-y_j||，其中j＝1,2,…,n；对计算结果distXX＝(distxx₁,distxx₂,…,distxx_n)，distXY＝(distxy₁,distxy₂,…,distxy_n)按由小到大的顺序进行排序，取距离最小的k个计算结果所对应的时间索引序列r,s；

步骤5.3：令r_n,j和s_n,j分别表示与x_n和y_n两个信号序列的最近的延迟向量的时间索引，j＝1,…,k，对于每个x_n，其最近的k个延迟向量的均方欧几里德距离定义为：

$R_n^{\left(k\right)}\left(X\right)=\frac{1}{k}\underset{j=1}{\overset{k}{\Σ}}{(x_n-x_{r_{n,j}})}^2---\left(1\right);$

同理，在时间序列Y上，x_n与最近的k个索引向量之间的均方欧几里德距离定义为：

$R_n^{\left(k\right)}\left(X\right|Y)=\frac{1}{k}\underset{j=1}{\overset{k}{\Σ}}{(x_n-x_{s_{n,j}})}^2---\left(2\right);$

在x_n中除x_n以外任意点的均方距离定义为：

$R_n\left(X\right)={(N-1)}^{-1}\underset{j\≠n}{\Σ}{(x_n-x_j)}^2---\left(3\right);$

如果X的活动独立于Y，则r_n,j和s_n,j就没有特定的关系，而且满足：

$R_n\left(X\right)\≈R_n^k\left(X\right|Y)>>R_n^k\left(X\right)---\left(4\right);$

反之，如果X和Y两个序列信号之间存在相互影响的关系，则满足：

$R_n\left(X\right)>>R_n^k\left(X\right|Y)\≈R_n^k\left(X\right)---\left(5\right);$

步骤5.4：定义独立性测量S为S^k(X|Y)，其计算方式如下：

$S^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\Σ}}\frac{R_n^{\left(k\right)}\left(X\right)}{R_n^{\left(k\right)}\left(X\right|Y)}---\left(6\right);$

由式4、式5可得0<S^k(X|Y)≤1，较低的值表示X和Y是独立的，S值越高则表示X和Y之间存在影响关系；

步骤5.5：利用几何平均的方法计算另一个独立性测量指数H，

$H^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}log\frac{R_n\left(X\right)}{R_n^{\left(k\right)}\left(X\right|Y)}---\left(7\right);$

步骤5.6：计算另一个独立性测量指数N，

$N^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\frac{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right|Y)}{R_n\left(X\right)}---\left(8\right);$

步骤5.7：计算独立性测量指数M，

$M^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\frac{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right|Y)}{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right)}---\left(9\right);$

以上几个测量中是以X序列为基准进行定义的，类似地，以Y序列为基准可以得到独立性测量S^(k)(Y|X)，H^(k)(Y|X)，N^(k)(Y|X)和M^(k)(Y|X)的定义。

作为优选，步骤5中所述NLI算法，是一种基于CPU并行化算法，其具体实现过程是每两个时间序列为一组，将每组时间窗长度数据作为一个待处理任务，放到任务队列中，当应用程序启动并执行创建线程池的时候，一定数量的任务被调用到线程池进行NLI运算，当有任务处理完成，线程空闲时，下一任务被调用，继续处理，如此执行直到所有任务执行完毕，即完成多通道数据的每两个通道的NLI算法，线程池才被销毁。

作为优选，步骤5中所述NLI算法，是一种基于CUDA的GPU并行算法，其具体实现包括以下子步骤：

第一步：对待处理数据进行相空间重构，此部分串行计算，每个通道待处理的N个数据经过重构变成M×D的矩阵，此步结束后得到N×M×D个数据，即X_n＝(x₁,x₂,…,x_N)，n＝N；其中，当嵌入为数和延迟时间不同时，嵌入维数取D；

第二步：对N个M×D的矩阵分别进行向量元素间差的平方运算，每个M×D的矩阵经过M次向量间计算化为M个M×D的矩阵，即X_(i,n)＝(x_n-x_j)²，i＝1,2,…,N，n,j＝1,2,…,M，此部分并行计算，依据通道数N和每个矩阵向量间计算次数M共分配M×N个线程块，每个线程块里完成M×D个数据的差的平方运算，所以为CUDA分配多个线程，每个线程处理一部分数据，此步结束后得到N×M×M×D个数据；

第三步：

3.1：依据第二步的结果，对差的平方求和，得到向量间的距离distXX，根据向量间距离公式，对N×M×M×D个数据每D个数求一个和，仍然使用M×N个线程块，每个线程块里D个线程，规约求和，循环M次，完成M×D个数据，最后得到N×M×M个数据，即N通道的distXX＝(distxx₁,distxx₂,…,distxx_M)；

3.2：依据第二步的结果，对差的平方求和，与3.1不同的是对N×M×M×D个数据每M×D个数据求一个和，并乘以一定的系数，得到的结果为x_n序列除x_n以外的任意点的均方距离R_n(X)，为后面的计算作准备；仍然使用10×2020个线程块，采用规约求和，由于规约求和要求待求数据为2的幂次，对于N×M个非2的幂次数据，分为两部分，把其中2的幂次的数采用规约求和的方法，剩下的采用普通求和方法，最后得到N×M个数据，乘以一定系数后的结果即是N通道数据的时间序列中任意点的均方距离R_n(X)；

第四步：对3.1的结果distXX采用双调排序法进行并行排序，由于次排序法要求待排数据为2的幂次，在排序前首先对N×M×M个数据进行补0，将N个M×M的方阵通过每行补0变成M×(M+D)的矩阵，N×M×(M+D)个数据分成两次完成并行排序，每次(N/2)×M×(M+D)个元素；每行待排数据为M+D，分配N×M个线程块，每个block里(M+D)/2个线程，以(M+D)为单位，每个线程块里负责(M+D)/2个元素的排序，每两个线程块完成(M+D)个元素的最终排序，其他各线程块负责(N/2)×M×(M+D)个数据其他元素的处理；此步结束后得到N×M×(M+D)数据，每M+D个数据按由大到小的顺序排列；

第五步：从排序后得到的矩阵中，找出除0之外的最小的k个值对应的索引，利用索引从重构后矩阵中找出与基向量最近的k个向量，得到N×M×D×k×N个元素，每个矩阵内部向量间有N个M×D×k，其中基向量为N个M×D的矩阵中的每个向量；计算矩阵内部和矩阵间每个向量与其最近的k个向量间的距离，即均方欧几里德距离；分配M×D个线程块，每个block里64个线程，每个线程处理N个数据，完成N×M×D×k×N个元素的差的平方和运算；得到N×N×M个元素，分别为和

本发明采用高性能并行计算中的CPU多线程技术和基于CUDA的GPGPU多线程方法开发并行NLI算法，将并行计算思想融入到算法中去，改进并提高算法的性能，辅助多通道信号的分析，将数据的处理从双通道扩展到多通道，并结合S估计器将数据的分析从局部到整体再到局部这样一个过程，量化信号一定区域范围内的同步强度，从原始数据中挖掘出更多更准确更有用的信息来分析多通道信号的同步问题。

与现有的传统的基于非线性互相关性方法(NLI)相比，本发明具有以下优点和有益效果：

(1)本发明通过对传统的NLI算法进行改进，使得算法能够快速的处理大规模高维非线性数据定量计算问题；

(2)从处理双通道数据变成能够处理多通道数据，从处理两两数据同步变成全局同步，并且支持定量测量；支持同步方向；且具有抗强噪声的特点；

(3)运用并行计算大大提升算法效率。应用基于CPU多线程的线程池技术和基于GPU多线程的CUDA技术，利用两种技术都达到了NLI算法并行化的目的，而且并行化后的算法比之前在效率和精度上都有了提高；

(4)同时本发明结合S估计器对多通道非线性数据进行分析，提取出数据更多的物理特性，为数据相关性分析提供了高效的解决方案。

附图说明

附图1：本发明实施例的流程图。

附图2：本发明实施例的NLI算法流程图。

附图3：本发明实施例的NLI算法CPU并行化流程图。

附图4：本发明实施例的NLI程序并行设计流程图。

附图5:本发明实施例的CUDA线程分配示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图1，本发明提供的一种非线性数据全局互相关性的并行定量计算方，包括以下步骤：

步骤一：初始化。多输入多通道信号shuru进行直接分解，两两通道数据最为原始数据进入下层循环，初始通道数据i＝1，j＝1。本实施例以输入第i通道原始数据shurui(图1中shuru1)和第j通道原始数据shuruj(图1中shuru2)为例，设置数据时间窗Epoch和滑动距离Overlap。

步骤二：根据初始输入的通道数据和算法参数Epoch和Overlap计算循环次数N，并设定初始循环次数为t＝1；

N＝length(shuru-Epoch)/(Epoch-Overlap)；

步骤三：按照时间窗从数据输入shuru1和shuru2里取数据，计算当前窗口数据偏移量(t-1)*(Epoch-Overlap)，得到当前数据X_t(k)和Y_t(k)；

步骤四：运用NLI算法计算两个信号数据的相关性，以X(t)和Y(t)为例，得

到独立性测量S、H、N、M，其主要流程图如图2所示：

步骤4.1：给定两个输入的信号数据X和Y，首先重构两个输入信号的延迟向量x_n＝(x_n,…,x_n-(m-1)τ)和y_n＝(y_n,…,y_n-(m-1)τ)，其中n＝1,…,N，m为嵌入维数，τ为预定义的延迟时间。N＝T-(m-1)τ，1≤n≤N。所有延迟向量组成了新的矩阵X＝(x₁,x₂,…,x_N)，Y＝(y₁,y₂,…,y_N)

步骤4.2：分别计算X，Y中两两延迟向量之间的距离distXX,distXY，即distxx_n＝||x_n-x_j||，distxy_n＝||y_n-y_j||，其中j＝1,2,…,n。对计算结果distXX＝(distxx₁,distxx₂,…,distxx_n)，distXY＝(distxy₁,distxy₂,…,distxy_n)按由小到大的顺序进行排序，取距离最小的k个计算结果所对应的时间索引序列r,s。

步骤4.3：令r_n,j和s_n,j，j＝1,…,k分别表示与x_n和y_n两个信号序列的最近的延迟向量的时间索引，对于每个x_n，其最近的k个延迟向量的均方欧几里德距离定义为：

$R_n^{\left(k\right)}\left(X\right)=\frac{1}{k}\underset{j=1}{\overset{k}{\&Sigma;}}{(x_n-x_{r_{n,j}})}^2---\left(1\right);$

同理，在时间序列Y上，x_n与最近的k个索引向量之间的均方欧几里德距离定义为：

$R_n^{\left(k\right)}\left(X\right|Y)=\frac{1}{k}\underset{j=1}{\overset{k}{\&Sigma;}}{(x_n-x_{s_{n,j}})}^2---\left(2\right);$

在x_n中除x_n以外任意点的均方距离定义为：

$R_n\left(X\right)={(N-1)}^{-1}\underset{j\&NotEqual;n}{\&Sigma;}{(x_n-x_j)}^2---\left(3\right);$

如果X的活动独立于Y，则r_n,j和s_n,j就没有特定的关系。而且满足：

$R_n\left(X\right)\&ap;R_n^k\left(X\right|Y)>>R_n^k\left(X\right)---\left(4\right);$

反之，如果X和Y两个序列信号之间存在相互影响的关系，则满足：

$R_n\left(X\right)>>R_n^k\left(X\right|Y)\&ap;R_n^k\left(X\right)---\left(5\right);$

步骤4.4：计算独立性测量S：

定义独立性测量S为S^k(X|Y)，其计算方式如下：

$S^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\frac{R_n^{\left(k\right)}\left(X\right)}{R_n^{\left(k\right)}\left(X\right|Y)}---\left(6\right);$

由式4、式5可得0<S^k(XY)≤1，较低的值表示X和Y是独立的，S值越高则表示X和Y之间存在影响关系。

步骤4.5：计算另一个独立性测量指数H，H独立性是利用几何平均的方法得到的。H的计算方式如下：

$H^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}log\frac{R_n\left(X\right)}{R_n^{\left(k\right)}\left(X\right|Y)}---\left(7\right);$

H在S基础上进行改进，但其值域没有进行归一化处理。

步骤4.6：计算另一个独立性测量指数N，N的计算采用了一种新的平均方法，其计算方式如下：

$N^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\frac{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right|Y)}{R_n\left(X\right)}---\left(8\right);$

N进行了归一化的处理，而且比S稳定。

步骤4.7：计算独立性测量指数M，根据定义只有当时，N＝1，这种情况不可能发生，独立性测量指数M定义为：

$M^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\frac{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right|Y)}{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right)}---\left(9\right);$

以上几个测量中是以X序列为基准进行定义的，类似地，以Y序列为基准可以得到独立性测量S^(k)(Y|X)，H^(k)(Y|X)，N^(k)(Y|X)和M^(k)(Y|X)的定义，由于NLI是一种不对称测量，所以可以根据二者的差值判断两个序列间驱动和响应的关系，这也是其他非线性测量方法所不具备的。

步骤五：对独立性测量S，H，N，M组成的相关性矩阵进行特征值特征向量分解。对数据矩阵X_M×T＝{x_i(k)}，i＝1,…,M，k＝1,…,T的协方差矩阵C进行特征值分解，其中M表示通道数量，T表示时间窗内的数据点数，获得特征值λ_i，并对其进行归一化，得到归一化特征值：

${{\mathrm{\&lambda;}}_i}^{\&prime;}=\frac{{\mathrm{\&lambda;}}_i}{\underset{i=1}{\overset{M}{\&Sigma;}}{\mathrm{\&lambda;}}_i}=\frac{{\mathrm{\&lambda;}}_i}{tr\left(C\right)}---\left(10\right);$

式中tr(C)——协方差矩阵的迹。

步骤六：S估计器同步性计算。S估计器的定义为：

$S=1+\frac{\underset{i=1}{\overset{M}{\&Sigma;}}{{\mathrm{\&lambda;}}_i}^{\&prime;}log\left({{\mathrm{\&lambda;}}_i}^{\&prime;}\right)}{1og\left(M\right)}\&Element;\&lsqb;0,1\&rsqb;---\left(11\right);$

S估计器是一种多变量的同步估计器，与同步强度成正比，非常适合用于分析高密度的多通道信号同步分析。

为提高NLI算法运行效率，本发明公开两种并行化设计方案。

请见图3，本实施例步骤四中的NLI算法基于CPU多线程的NLI并行设计子步骤如下：

根据算法和线程池技术的特点，每两个时间序列为一组，将每组时间窗长度数据作为一个待处理任务，放到任务队列中，当应用程序启动并执行创建线程池的时候，一定数量的任务被调用到线程池进行NLI运算，当有任务处理完成，线程空闲时，下一任务被调用，继续处理，如此执行直到所有任务执行完毕，即完成多通道数据的每两个通道的NLI算法，线程池才被销毁。程序执行示意图如图3所示。CPU多线程并行从外部粗粒度水平对NLI算法进行并行化。

请见图4，本实施例步骤四中的NLI算法基于CUDA的GPU并行设计子步骤如下：

第一步：对待处理数据进行相空间重构，此部分串行计算，每个通道待处理的N个数据经过重构变成M×D(嵌入为数和延迟时间不同，嵌入维数取D)的矩阵，此步结束后得到N×M×D个数据，即X_n＝(x₁,x₂,…,x_N)，n＝N。

第二步：对N个M×D的矩阵分别进行向量元素间差的平方运算，每个M×D的矩阵经过M次向量间计算化为M个M×D的矩阵，即X_(i,n)＝(x_n-x_j)²，i＝1,2,…,N，n,j＝1,2,…,M，此部分并行计算，依据通道数N和每个矩阵向量间计算次数M共分配M×N个线程块，每个线程块里完成M×D个数据的差的平方运算，所以为CUDA分配多个线程，每个线程处理一部分数据，此步结束后得到N×M×M×D个数据.

第三步：

(1)依据第二步的结果，对差的平方求和，得到向量间的距离distXX，根据向量间距离公式，对N×M×M×D个数据每D个数求一个和，仍然使用M×N个线程块，每个线程块里D个线程，规约求和，循环M次，完成M×D个数据，最后得到N×M×M个数据，即N通道的distXX＝(distxx₁,distxx₂,…,distxx_M)；

(2)依据第二步的结果，对差的平方求和，与(1)不同的是对N×M×M×D个数据每M×D个数据求一个和，并乘以一定的系数，得到的结果为x_n序列除x_n以外的任意点的均方距离R_n(X)，为后面的计算作准备。仍然使用10×2020个线程块，采用规约求和，由于规约求和要求待求数据为2的幂次，对于N×M个非2的幂次数据，可以分为两部分，把其中2的幂次的数采用规约求和的方法，剩下的采用普通求和方法，最后得到N×M个数据，乘以一定系数后的结果即是N通道数据的时间序列中任意点的均方距离R_n(X)；

第四步：对第三步(1)的结果distXX采用双调排序法进行并行排序，由于次排序法要求待排数据为2的幂次，在排序前首先对N×M×M个数据进行补0，将N个M×M的方阵通过每行补0变成M×(M+D)的矩阵，N×M×(M+D)个数据分成两次完成并行排序，每次(N/2)×M×(M+D)个元素。每行待排数据为M+D，分配N×M个线程块，每个block里(M+D)/2个线程，以(M+D)为单位，每个线程块里负责(M+D)/2个元素的排序，每两个线程块完成(M+D)个元素的最终排序，其他各线程块负责(N/2)×M×(M+D)个数据其他元素的处理。此步结束后得到N×M×(M+D)数据，每M+D个数据按由大到小的顺序排列；

第五步：从排序后得到的矩阵中，找出除0之外的最小的k个值对应的索引，利用索引从重构后矩阵中找出与基向量(N个M×D的矩阵中的每个向量)最近的k个向量，得到N×M×D×k×N个元素：每个矩阵内部向量间有N个M×D×k。计算矩阵内部和矩阵间每个向量与其最近的k个向量间的距离，即均方欧几里德距离。分配M×D个线程块，每个block里64个线程，每个线程处理N个数据，完成N×M×D×k×N个元素的差的平方和运算。得到N×N×M个元素，分别为和

本实施例采用的是CUDA4.2(统一计算设备架构)，而传统的小波相干性方法是基于C语言实现的。实验的硬件环境是NVIDIAGeForceGTX480显卡，处理器是Intel(R)Core(TM)i7-2600。待处理的原始数据为10通道数据，采样率1024Hz，数据长度2000s，对原始数据进行降采样，降频后的采样频率为256Hz，采用移动窗的方法，窗口长度为8s，滑动5s，即每次待处理的数据为10×2048个，NLI并行设计如下：

第一步：对待处理数据进行相空间重构，此部分串行计算，每个通道待处理的2048个数据经过重构变成2020×8(嵌入为数和延迟时间不同，得到的矩阵不同，本实施例中延迟时间取4，嵌入维数取8)的矩阵(实际程序中将矩阵化为向量参加运算，后面几步涉及矩阵的运算都是做同样处理)，此步结束后得到10×2020×8个数据，即X_n＝(x₁,x₂,…,x_N)，n＝10。

第二步：对10个2020×8的矩阵分别进行向量元素间差的平方运算，每个2020×8的矩阵经过2020次向量间计算化为2020个2020×8的矩阵，即X_(i,n)＝(x_n-x_j)²，i＝1,2,…,10，n,j＝1,2,…,2020，此部分并行计算，依据通道数10和每个矩阵向量间计算次数2020共分配2020×10个线程块，每个线程块里完成2020×8个数据的差的平方运算，所以分配160个线程，每个线程处理101个数据，此步结束后得到10×2020×2020×8个数据，图5是CUDA的线程分配示意图；

第三步：

(1)依据第二步的结果，对差的平方求和，得到向量间的距离distXX，根据向量间距离公式，对10×2020×2020×8个数据每8个数求一个和，仍然使用10×2020个线程块，每个线程块里8个线程，规约求和，循环2020次，完成2020×8个数据，最后得到10×2020×2020个数据，即10通道的distXX＝(distxx₁,distxx₂,…,distxx₂₀₂₀)；

(2)依据第二步的结果，对差的平方求和，与(1)不同的是对10×2020×2020×8个数据每2020×8个数据求一个和，并乘以一定的系数，得到的结果为x_n序列除x_n以外的任意点的均方距离R_n(X)，为后面的计算作准备。仍然使用10×2020个线程块，采用规约求和，由于规约求和要求待求数据为2的幂次，对于2020×8个非2的幂次数据，可以分为两部分，把其中2的幂次的数采用规约求和的方法，剩下的采用普通求和方法，具体线程分配为，首先每个block里分配32个线程，规约求和，循环505次，得到10×2020×505个数据，然后重新分配线程，每个线程块里1个线程，每个线程完成505个数据的求和，最后得到10×2020个数据，乘以一定系数后的结果即是10通道数据的时间序列中任意点的均方距离R_n(X)；

第四步：对第三步(1)的结果distXX采用双调排序法进行并行排序，由于次排序法要求待排数据为2的幂次，在排序前首先对10×2020×2020个数据进行补0，将10个2020×2020的方阵通过每行补0变成2020×2048的矩阵，10×2020×2048个数据分成两次完成并行排序，每次5×2020×2048个元素。N＝5×2020×2048，每行待排数据为2048，分配N/1024个线程块，每个block里1024/2个线程，以2048为单位，每个线程块里负责1024个元素的排序，每两个线程块完成2048个元素的最终排序，其他各线程块负责5×2020×2048个数据其他元素的处理。此步结束后得到10×2020×2048数据，每2048个数据按由大到小的顺序排列；

第五步：从排序后得到的矩阵中，找出除0之外的最小的k个值对应的索引，利用索引从重构后矩阵中找出与基向量(10个2020×8的矩阵中的每个向量)最近的k个向量，本实施例中k＝8，得到10×2020×8×8×10个元素：每个矩阵内部向量间有10个2020×8×8；各矩阵间向量有90个2020×8×8。计算矩阵内部和矩阵间每个向量与其最近的8个向量间的距离，即均方欧几里德距离。分配2020×10个线程块，每个block里64个线程，每个线程处理10个数据，完成10×2020×8×8×10个元素的差的平方和运算。得到10×10×2020个元素，分别为和

第六步：计算S^(k)(Y|X)，H^(k)(Y|X)，N^(k)(Y|X)和M^(k)(Y|X)并评估其同步性。

(1)原始数据在MATLAB上执行，记作M-NLI；

(2)基于C语言的单线程NLI执行，记作SC-NLI；

(3)基于线程池的多线程实现，线程数为4，记作DC-NLI；

(4)基于CUDA的NLI并行实现，记作G-NLI。

对每个通道512000个数据点的10通道数据进行处理分析，应用滑动窗方法，数据序列被划分为399个时间窗，各版本测试时间如表1所示。

表1程序执行时间比较

从表1可以看出，应用不同方法在样计算机操作环境下对相同数据处理所用的时间从G-NLI到M-NLI近似于指数增长，NLI算法四种实现中MATLAB程序耗时约833300秒，相当于9.6天，对于算法应用来说这个时间是不能容忍的；C语言程序耗时约127800秒，虽然比MATLAB语言执行效率提高了，但对于NLI算法的研究者和使用者来说，这个时间相对来说还是比较长；而经过改进后的CPU多线程和GPU多线程两种并行实现方法在执行效率上有了很明显的提高，特别是CUDA版本的NLI算法，耗时仅735秒，明显的提高了算法的执行效率。G-NLI相比于M-NLI提高了三个数量级，约1133倍，是SC-NLI的173倍，是DC-NLI的47倍。NLI算法的并行化设计明显提高了该算法的执行效率，缩短了研究周期，同时将NLI算法直接应用于多通道，扩大了研究尺度，使得算法应用的更加广泛。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (4)

1.一种非线性数据全局互相关性的并行定量计算方法，其特征在于，包括以下步骤：

步骤1：对M通道原始脑电数据shuru进行直接分解，M＞2，以两两通道数据最为原始数据进入下层循环；设定初始通道i＝1，j＝1；

步骤2：取第i通道原始数据shurui；

步骤3：取第j通道原始数据shuruj；

步骤4：设置数据时间窗Epoch和滑动距离Overlap；

步骤5：设定初始循环次数为t＝1，计算循环次数N：

N＝length(shuru-Epoch)/(Epoch-Overlap)

步骤4：按照时间窗从数据输入shurui和shuruj里取数据，计算当前窗口数据偏移量(t-1)*(Epoch-Overlap)，得到当前数据X_t(k)和Y_t(k)；

步骤5：运用NLI算法计算两个数据的相关性，得到独立性测量S、H、N、M；

步骤6：判断t＞N是否成立；

若是，则顺序执行下述步骤6；

若否，则t＝t+1，并回转执行所述步骤4；

步骤7：判断j＞M是否成立；

若是，则顺序执行下述步骤8；

若否，则j＝j+1，并回转执行所述步骤3；

步骤8：判断i＞M是否成立；

若是，则顺序执行下述步骤9；

若否，则i＝i+1，并回转执行所述步骤2；

步骤9：对独立性测量S、H、N、M组成的相关性矩阵进行特征值特征向量分解；

对数据矩阵X_M×T＝{x_i(k)}，i＝1,…,M，k＝1,…,T的协方差矩阵C进行特征值分解，其中M表示通道数量，T表示时间窗内的数据点数，获得特征值λ_i，并对其进行归一化，得到归一化特征值：

${{\mathrm{\&lambda;}}_i}^{\&prime;}=\frac{{\mathrm{\&lambda;}}_i}{\underset{i=1}{\overset{M}{\&Sigma;}}{\mathrm{\&lambda;}}_i}=\frac{{\mathrm{\&lambda;}}_i}{tr\left(C\right)};$

式中tr(C)为协方差矩阵的迹；

步骤10：基于特征值特征向量计算S估计器，S估计器的定义为：

$S=1+\frac{\underset{i=1}{\overset{M}{\&Sigma;}}{{\mathrm{\&lambda;}}_i}^{\&prime;}log\left({{\mathrm{\&lambda;}}_i}^{\&prime;}\right)}{\log \left(M\right)}\&Element;\&lsqb;0,1\&rsqb;.$

2.根据权利要求1所述的非线性数据全局互相关性的并行定量计算方法，其特征在于，步骤5中所述NLI算法，其具体实现过程包括以下子步骤：

步骤5.1：给定两个输入的信号数据X和Y，重构两个输入信号的延迟向量x_n＝(x_n,…,x_n-(m-1)τ)和y_n＝(y_n,…,y_n-(m-1)τ)，其中n＝1,…,N，m为嵌入维数，τ为预定义的延迟时间；N＝T-(m-1)τ，1≤n≤N；所有延迟向量组成了新的矩阵X＝(x₁,x₂,…,x_N)，Y＝(y₁,y₂,…,y_N)；

步骤5.2：分别计算X，Y中两两延迟向量之间的距离distXX,distXY，即distxx_n＝||x_n-x_j||，distxy_n＝||y_n-y_j||，其中j＝1,2,…,n；对计算结果distXX＝(distxx₁,distxx₂,…,distxx_n)，distXY＝(distxy₁,distxy₂,…,distxy_n)按由小到大的顺序进行排序，取距离最小的k个计算结果所对应的时间索引序列r,s；

步骤5.3：令r_n,j和s_n,j分别表示与x_n和y_n两个信号序列的最近的延迟向量的时间索引，j＝1,…,k，对于每个x_n，其最近的k个延迟向量的均方欧几里德距离定义为：

$R_n^{\left(k\right)}\left(X\right)=\frac{1}{k}\underset{j=1}{\overset{k}{\&Sigma;}}{(x_n-x_{r_{n,j}})}^2---\left(1\right);$

同理，在时间序列Y上，x_n与最近的k个索引向量之间的均方欧几里德距离定义为：

$R_n^{\left(k\right)}\left(X\right|Y)=\frac{1}{k}\underset{j=1}{\overset{k}{\&Sigma;}}{(x_n-x_{s_{n,j}})}^2---\left(2\right);$

在x_n中除x_n以外任意点的均方距离定义为：

$R_n\left(X\right)={(N-1)}^{-1}\underset{j\&NotEqual;n}{\&Sigma;}{(x_n-x_j)}^2---\left(3\right);$

如果X的活动独立于Y，则r_n,j和s_n,j就没有特定的关系，而且满足：

$R_n\left(X\right)\&ap;R_n^k\left(X\right|Y)>>R_n^k\left(X\right)---\left(4\right);$

反之，如果X和Y两个序列信号之间存在相互影响的关系，则满足：

$R_n\left(X\right)>>R_n^k\left(X\right|Y)\&ap;R_n^k\left(X\right)---\left(5\right);$

步骤5.4：定义独立性测量S为S^k(X|Y)，其计算方式如下：

$S^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\frac{R_n^{\left(k\right)}\left(X\right)}{R_n^{\left(k\right)}\left(X\right|Y)}---\left(6\right);$

由式4、式5可得0<S^k(X|Y)≤1，较低的值表示X和Y是独立的，S值越高则表示X和Y之间存在影响关系；

步骤5.5：利用几何平均的方法计算另一个独立性测量指数H，

$H^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}log\frac{R_n\left(X\right)}{R_n^{\left(k\right)}\left(X\right|Y)}---\left(7\right);$

步骤5.6：计算另一个独立性测量指数N，

$N^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\frac{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right|Y)}{R_n\left(X\right)}---\left(8\right);$

步骤5.7：计算独立性测量指数M，

$M^{\left(k\right)}\left(X\right|Y)=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\frac{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right|Y)}{R_n\left(X\right)-R_n^{\left(k\right)}\left(X\right)}---\left(9\right);$

以上几个测量中是以X序列为基准进行定义的，类似地，以Y序列为基准可以得到独立性测量S^(k)(Y|X)，H^(k)(Y|X)，N^(k)(Y|X)和M^(k)(Y|X)的定义。

3.根据权利要求1或2所述的非线性数据全局互相关性的并行定量计算方法，其特征在于，步骤5中所述NLI算法，是一种基于CPU并行化算法，其具体实现过程是每两个时间序列为一组，将每组时间窗长度数据作为一个待处理任务，放到任务队列中，当应用程序启动并执行创建线程池的时候，一定数量的任务被调用到线程池进行NLI运算，当有任务处理完成，线程空闲时，下一任务被调用，继续处理，如此执行直到所有任务执行完毕，即完成多通道数据的每两个通道的NLI算法，线程池才被销毁。

4.根据权利要求1或2所述的非线性数据全局互相关性的并行定量计算方法，其特征在于，步骤5中所述NLI算法，是一种基于CUDA的GPU并行算法，其具体实现包括以下子步骤：

第一步：对待处理数据进行相空间重构，此部分串行计算，每个通道待处理的N个数据经过重构变成M×D的矩阵，此步结束后得到N×M×D个数据，即X_n＝(x₁,x₂,…,x_N)，n＝N；其中，当嵌入为数和延迟时间不同时，嵌入维数取D；

第二步：对N个M×D的矩阵分别进行向量元素间差的平方运算，每个M×D的矩阵经过M次向量间计算化为M个M×D的矩阵，即X_(i,n)＝(x_n-x_j)²，i＝1,2,…,N，n,j＝1,2,…,M，此部分并行计算，依据通道数N和每个矩阵向量间计算次数M共分配M×N个线程块，每个线程块里完成M×D个数据的差的平方运算，所以为CUDA分配多个线程，每个线程处理一部分数据，此步结束后得到N×M×M×D个数据；

第三步：

3.1：依据第二步的结果，对差的平方求和，得到向量间的距离distXX，根据向量间距离公式，对N×M×M×D个数据每D个数求一个和，仍然使用M×N个线程块，每个线程块里D个线程，规约求和，循环M次，完成M×D个数据，最后得到N×M×M个数据，即N通道的distXX＝(distxx₁,distxx₂,…,distxx_M)；

3.2：依据第二步的结果，对差的平方求和，与3.1不同的是对N×M×M×D个数据每M×D个数据求一个和，并乘以一定的系数，得到的结果为x_n序列除x_n以外的任意点的均方距离R_n(X)，为后面的计算作准备；仍然使用10×2020个线程块，采用规约求和，由于规约求和要求待求数据为2的幂次，对于N×M个非2的幂次数据，分为两部分，把其中2的幂次的数采用规约求和的方法，剩下的采用普通求和方法，最后得到N×M个数据，乘以一定系数后的结果即是N通道数据的时间序列中任意点的均方距离R_n(X)；

第四步：对3.1的结果distXX采用双调排序法进行并行排序，由于次排序法要求待排数据为2的幂次，在排序前首先对N×M×M个数据进行补0，将N个M×M的方阵通过每行补0变成M×(M+D)的矩阵，N×M×(M+D)个数据分成两次完成并行排序，每次(N/2)×M×(M+D)个元素；每行待排数据为M+D，分配N×M个线程块，每个block里(M+D)/2个线程，以(M+D)为单位，每个线程块里负责(M+D)/2个元素的排序，每两个线程块完成(M+D)个元素的最终排序，其他各线程块负责(N/2)×M×(M+D)个数据其他元素的处理；此步结束后得到N×M×(M+D)数据，每M+D个数据按由大到小的顺序排列；

第五步：从排序后得到的矩阵中，找出除0之外的最小的k个值对应的索引，利用索引从重构后矩阵中找出与基向量最近的k个向量，得到N×M×D×k×N个元素，每个矩阵内部向量间有N个M×D×k，其中基向量为N个M×D的矩阵中的每个向量；计算矩阵内部和矩阵间每个向量与其最近的k个向量间的距离，即均方欧几里德距离；分配M×D个线程块，每个block里64个线程，每个线程处理N个数据，完成N×M×D×k×N个元素的差的平方和运算；得到N×N×M个元素，分别为和