CN103093085B - 基于典型相关分析的稳态诱发电位的分析方法 - Google Patents

Abstract

本发明公开了一种应用于嵌入式处理器中的脑电稳态诱发电位的分析方法。所述方法包括：S1、选取一个待分析的频率，生成参考信号；S2、利用典型相关分析方法，将所述参考信号与待分析的稳态诱发电位进行定点运算，得到所述参考信号与所述稳态诱发电位的最大相关系数，从而得到所述稳态诱发电位在所述频率处的响应强度。该方法能够使定点嵌入式处理器实现对脑电稳态诱发电位的在线分析；基于该方法构建的嵌入式系统，对比传统的傅里叶分析方法具有更加高效和准确的优点。

Description

基于典型相关分析的稳态诱发电位的分析方法

技术领域

本发明涉及数字信号处理技术领域，具体涉及嵌入式处理器技术，特别涉及一种基于典型相关分析的稳态诱发电位的分析方法。

背景技术

如果以一定频率的闪烁刺激人眼，则可以在枕区的初级视觉皮层处，通过头皮脑电（EEG，Electroencephalogram）检测到与刺激频率相应的基频以及各次谐波成分，这种脑电信号被称为稳态视觉诱发电位（SSVEP，SteadyStateVisualEvokedPotential）。与之类似，如果给予受试者一定调制频率下的听觉刺激，通过头皮脑电，可以检测到听觉稳态响应（ASSR，AuditorySteadyStateResponse）。

传统的稳态诱发电位的分析方法主要针对单个电极的频域功率谱，采用F检验，判断在刺激频率以及各次谐波处是否出现了显著的频率成分。2007年，Lin等人首先将典型相关分析（CCA，CanonicalCorrelationAnalysis）的方法引入到SSVEP信号的检测中，随后，2009年，Bin等人成功构造出了基于CCA方法的SSVEP范式下脑机接口系统。

CCA方法是Hotelling于1936年提出的。它的核心思想是寻找两组信号之间的最大线性相关性。这是一种成熟的方法，已经广泛应用于气象、医学、环境、农业等领域。将这种方法应用到对脑电稳态诱发电位的分析中，给出了一种针对特定几个频率点的高效的识别方法，可以快速准确地判断出特定频率处的响应是否存在以及响应的大小。

过去对于脑电信号的分析和处理，都是基于计算机的浮点运算。为了将对稳态诱发电位的分析方法移植到嵌入式处理器中，本发明提出了一种基于定点CCA的稳态诱发电位分析方法。通过这种方法，可以在任意嵌入式处理器中实现对稳态诱发电位的在线处理，判断出稳态响应是否存在以及推断出刺激频率。

发明内容

（一）所要解决的技术问题

本发明的目的在于提供一种基于典型相关分析的稳态诱发电位的分析方法，从而使定点嵌入式处理器能够实现对脑电稳态诱发电位的在线分析。

（二）技术方案

为了解决上述技术问题，本发明提出了一种基于典型相关分析的稳态诱发电位的分析方法，所述方法包括以下步骤：

S1、选取一个待分析的频率，生成参考信号；

S2、利用典型相关分析方法，将所述参考信号与待分析的稳态诱发电位进行定点运算，得到所述参考信号与所述稳态诱发电位的最大相关系数，从而得到所述稳态诱发电位在所述频率处的响应强度。

可选的，步骤S2之后进一步包括步骤：

S3、选取另一个待分析的频率，重复步骤S1-S2，得到所述稳态诱发电位在另一个频率处的响应强度；

S4、重复步骤S3，直至得到所述稳态诱发电位在全部待分析的频率处的响应强度，从而判断出所述稳态诱发电位响应强度最高的频率。

可选的，步骤S2中，所述定点运算采用的数据格式为2补码表示的小数。

可选的，步骤S2具体包括：

S2-1、将所述参考信号及所述稳态诱发电位进行零均值化处理，得到向量X和向量Y；

S2-2、将向量X和向量Y进行QR分解，得到矩阵Q1和矩阵Q2；

S2-3、将矩阵Q1和矩阵Q2相乘，得到矩阵Q；

S2-4、将矩阵Q进行SVD分解，得到两个空域滤波器矩阵以及最大相关系数，所述最大相关系数即为所述参考信号与所述稳态诱发电位的最大相关系数。

可选的，步骤S2-1中，所述稳态诱发电位的零均值化处理是通过FIR滤波器实现的。

可选的，步骤S2-2具体包括：

通过Schimidt正交化实现向量X和向量Y的QR分解。

可选的，步骤S2-4中，矩阵Q的SVD分解是通过Jacobi迭代实现的。

可选的，步骤S1具体包括：

选取一个待分析的频率，生成基频加倍频的参考信号。

（三）有益效果

本发明主要针对脑电稳态诱发电位给出了一种可在定点嵌入式处理器中实现的分析方法。该方法对原始输入数据的要求不高，且可根据具体的数据长度调整实际使用的数据位宽，平衡功耗、计算时间与计算精度的要求。对稳态诱发电位进行分析的目的是判断某些特定频率点处的响应，该方法对比传统的傅里叶分析方法具有更加高效和准确的优点。传统的傅里叶分析在数据点数达到一定的时候，无法获得很高的精度，且速度较慢。本发明提出的方法在数据点数增加时，只需适当调整整体数据的位宽，而计算的速度不会受到特别大的影响。

附图说明

图1是本发明提出的技术方案的基本流程图。

图2是本发明一个实施例中CCA方法的基本流程图。

图3是本发明一个实施例中定点运算相对于浮点运算的QR分解结果的精度误差示意图。

图4是本发明一个实施例中定点运算相对于浮点运算的SVD分解左边矩阵的误差示意图。

具体实施方式

本发明提出了一种基于典型相关分析的稳态诱发电位的分析方法，如图1所示，所述方法包括以下步骤：

S1、选取一个待分析的频率，生成参考信号；

S2、利用典型相关分析方法，将所述参考信号与待分析的稳态诱发电位进行定点运算，得到所述参考信号与所述稳态诱发电位的最大相关系数，从而得到所述稳态诱发电位在所述频率处的响应强度。

优选的，步骤S2之后进一步包括步骤：

S3、选取另一个待分析的频率，重复步骤S1-S2，得到所述稳态诱发电位在另一个频率处的响应强度；

S4、重复步骤S3，直至得到所述稳态诱发电位在全部待分析的频率处的响应强度，从而判断出所述稳态诱发电位响应强度最高的频率。

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。

图2表示了一种CCA方法的基本流程图。能够证明，通过图2所示的流程，实现CCA方法，在对输入的数据进行合理的限制之后，整个计算过程所产生的数据不会超过[-1,1）这个范围。

1、零均值化过程不会溢出，只要待零均值化的向量，每个元素的值都不大于0.5。在本发明的一个实施例中，选择的原始脑电信号来自于24位模数转换器的采样结果，选择的数据位宽为32位，从而这个步骤的要求一定可以满足。

2、Schimidt正交化方法实现的QR分解过程不会溢出，只要保证整个过程中原矩阵的列向量的模均不会超过1。类似于步骤1，因为数据只有24位，而选择的位宽为32位，一次分析用数据长度为400点，从而这个步骤对于向量模的要求也能满足。

3、矩阵相乘的过程不会发生溢出，因为对于实际的数据，采样时的噪声等问题，不可能出现两个待相乘的矩阵完全相同的情况。

4、通过Jacobi迭代的方法实现的SVD分解过程不会溢出，因为待分解的矩阵来自于两个正交矩阵相乘的结果，整个分解过程矩阵的结构保持不变。

优选的，在进行定点CCA计算时，算法采用的数据格式为利用2补码表示的小数。具体来说即，对于位数为n+1的数据，有效部分为n位，最高位为符号位。首先将数据看成2补码的整数M，然后该数据实际表示的实数大小为M/2ⁿ。本实施例中选择n=31。

根据图2，步骤S2中的典型相关分析方法具体包括：

S2-1、将所述参考信号及所述稳态诱发电位进行零均值化处理，得到向量X和向量Y。在本实施例中，通过FIR滤波器替代零均值化过程。脑电数据的采样率为200Hz，需要进行50Hz限波，以及滤除低频和高频成分。本实施例中选用平均滤波器，采用[1,0,0,0,-1]的滤波器系数。

S2-2、将向量X和向量Y进行QR分解，得到矩阵Q1和矩阵Q2。在本实施例中，通过Schimidt正交化方法实现QR分解。

S2-3、将矩阵Q1和矩阵Q2相乘，得到矩阵Q；

S2-4、将矩阵Q进行SVD分解，得到两个空域滤波器矩阵以及最大相关系数，所述最大相关系数即为所述参考信号与所述稳态诱发电位的最大相关系数。在本实施例中，通过Jacobi迭代的方法实现SVD分解。

本实施例中的嵌入式设备选择为TI公司的C5000系列数字信号处理器（DSP，DigitalSignalProcessor），TMS320C5515。编程语言采用C语言调用汇编语言函数的形式。最底层的数学计算采用汇编语言完成，CCA方法的整合利用C语言调用汇编语言函数的方式完成。

利用汇编语言完成的基本运算包括：

1、加法和减法，基本的二补码数据的加减法。

2、乘法，利用两倍于数据位宽的累加器和乘法器完成乘法的计算，如果嵌入式处理器的数据位宽不足，采用乘数与被乘数分别拆分成两部分再进行4次相乘的手段。

3、除法，首先求取除数的倒数，然后用被除数乘以除数的倒数。本实施例中，倒数的求取采用逐位比较的方法。

4、三角函数的计算，采用基于Taylor展开的方法。本实施例中，Taylor展开到7阶。

本实施例分析的SSVEP信号的离线数据来自于一项具有6个目标的SSVEP实验。6个目标的闪烁频率分别为7、8、9、10、11、12。分析的目的是，从受试者的脑电信号中，判断出这6个频率点中哪个频率点的响应最强。CCA方法分析的两组输入信号分别来自于原始的脑电数据，以及构造的特定频率的参考信号。参考信号的构造方法如下式所示，选定一个基频ft，然后构造最高位N_h的倍频：

$Y_f=\left[\begin{array}{c}\sin \left(2\mathrm{\πft}\right)\\ \cos \left(2\mathrm{\πft}\right)\\ \·\\ \·\\ \·\\ \sin \left(2\mathrm{\π}{N}_h\mathrm{ft}\right)\\ \cos \left(2\mathrm{\π}{N}_h\mathrm{ft}\right)\end{array}\right]$

参考信号因为是选定的，所以可以无需执行完整的在线CCA流程。具体地说，参考信号首先在计算机上通过浮点运算，执行完QR分解，再将数据格式转化为与定点算法相同的数据格式，直接存储到嵌入式处理器中，然后与处理器在线完成计算的脑电数据的QR分解结果一起，完成后续的计算。

图3和图4显示了通过上述方法计算得到的结果与计算机浮点运算的结果之间的误差。其中图3表示脑电数据QR分解的精度误差，具体来说是，144段400点长的脑电数据，进行QR分解，结果的正交性误差。图4展示了SVD分解过程的精度误差，具体来说，是864个8*6的矩阵进行SVD分解后，左边矩阵的正交性的误差。与此同时，整个定点CCA过程，造成的最终最大相关系数误差不超过10^-5。

通过比较6个不同频率点处的相关系数的大小，可以判断出，在哪个频率点上获得了最大的响应，从而判断出受试者在实验中所注视的目标的频率。

以上所述仅是本发明的优选实施方式，应当指出，对于本领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和替换，这些改进和替换也应视为本发明的保护范围。

Claims (6)

1.一种基于典型相关分析的稳态诱发电位的分析方法，其特征在于，所述方法包括以下步骤：

S1、选取一个待分析的频率，生成参考信号；

S2、利用典型相关分析方法，将所述参考信号与待分析的稳态诱发电位进行定点运算，得到所述参考信号与所述稳态诱发电位的最大相关系数，从而得到所述稳态诱发电位在所述频率处的响应强度；

所述步骤S2中，所述定点运算采用的数据格式为2补码表示的小数，具体包括：

S2-1、将所述参考信号及所述稳态诱发电位进行零均值化处理，得到向量X和向量Y；

S2-2、将向量X和向量Y进行QR分解，得到矩阵Q1和矩阵Q2；

S2-3、将矩阵Q1和矩阵Q2相乘，得到矩阵Q；

S2-4、将矩阵Q进行SVD分解，得到两个空域滤波器矩阵以及最大相关系数，所述最大相关系数即为所述参考信号与所述稳态诱发电位的最大相关系数。

2.根据权利要求1所述的方法，其特征在于，步骤S2之后进一步包括步骤：

S3、选取另一个待分析的频率，重复步骤S1-S2，得到所述稳态诱发电位在另一个频率处的响应强度；

S4、重复步骤S3，直至得到所述稳态诱发电位在全部待分析的频率处的响应强度，从而判断出所述稳态诱发电位响应强度最高的频率。

3.根据权利要求1所述的方法，其特征在于，步骤S2-1中，所述稳态诱发电位的零均值化处理是通过FIR滤波器实现的。

4.根据权利要求1所述的方法，其特征在于，步骤S2-2具体包括：

通过Schimidt正交化实现向量X和向量Y的QR分解。

5.根据权利要求1所述的方法，其特征在于，步骤S2-4中，矩阵Q的SVD分解是通过Jacobi迭代实现的。

6.根据权利要求1所述的方法，其特征在于，步骤S1具体包括：

选取一个待分析的频率，生成基频加倍频的参考信号。