CN114218986A - 基于eeg脑电信号数据的状态分类方法 - Google Patents
基于eeg脑电信号数据的状态分类方法 Download PDFInfo
- Publication number
- CN114218986A CN114218986A CN202111502840.4A CN202111502840A CN114218986A CN 114218986 A CN114218986 A CN 114218986A CN 202111502840 A CN202111502840 A CN 202111502840A CN 114218986 A CN114218986 A CN 114218986A
- Authority
- CN
- China
- Prior art keywords
- frequency
- matrix
- signal data
- data
- electroencephalogram
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 83
- 239000013598 vector Substances 0.000 claims abstract description 70
- 238000000605 extraction Methods 0.000 claims abstract description 22
- 230000008569 process Effects 0.000 claims abstract description 16
- 230000004927 fusion Effects 0.000 claims abstract description 14
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 11
- 238000013145 classification model Methods 0.000 claims abstract description 11
- 238000012360 testing method Methods 0.000 claims abstract description 8
- 238000012549 training Methods 0.000 claims abstract description 7
- 238000007781 pre-processing Methods 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 113
- 230000003595 spectral effect Effects 0.000 claims description 30
- 238000005070 sampling Methods 0.000 claims description 25
- 238000012880 independent component analysis Methods 0.000 claims description 18
- 238000002156 mixing Methods 0.000 claims description 18
- 230000014509 gene expression Effects 0.000 claims description 16
- 108010076504 Protein Sorting Signals Proteins 0.000 claims description 12
- 238000001914 filtration Methods 0.000 claims description 11
- 230000009467 reduction Effects 0.000 claims description 11
- 238000012545 processing Methods 0.000 claims description 10
- 238000012706 support-vector machine Methods 0.000 claims description 10
- 238000005457 optimization Methods 0.000 claims description 8
- 230000002087 whitening effect Effects 0.000 claims description 6
- 230000004424 eye movement Effects 0.000 claims description 5
- 238000004458 analytical method Methods 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 3
- 230000003340 mental effect Effects 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 230000000875 corresponding effect Effects 0.000 description 19
- 230000006870 function Effects 0.000 description 11
- 239000000203 mixture Substances 0.000 description 6
- 238000011160 research Methods 0.000 description 5
- 210000004556 brain Anatomy 0.000 description 4
- 238000000354 decomposition reaction Methods 0.000 description 3
- 210000001508 eye Anatomy 0.000 description 3
- 210000003128 head Anatomy 0.000 description 3
- 230000010365 information processing Effects 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000002790 cross-validation Methods 0.000 description 2
- 238000004880 explosion Methods 0.000 description 2
- 230000015654 memory Effects 0.000 description 2
- 238000003909 pattern recognition Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000004397 blinking Effects 0.000 description 1
- 230000007177 brain activity Effects 0.000 description 1
- 210000005252 bulbus oculi Anatomy 0.000 description 1
- 230000036992 cognitive tasks Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000000241 respiratory effect Effects 0.000 description 1
- 210000004761 scalp Anatomy 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000003936 working memory Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/12—Classification; Matching
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Mathematical Physics (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Artificial Intelligence (AREA)
- Mathematical Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Signal Processing (AREA)
- Computing Systems (AREA)
- Evolutionary Computation (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明提供一种基于EEG脑电信号数据的状态分类方法,包括以下步骤:S1,获取不同分类的EEG脑电信号数据、以及相应的垂直和水平眼电信号数据;S2,EEG脑电信号数据预处理;S3,基于过程的特征提取;S4,使用频‑空特征向量训练分类模型;S5,对样本空间的不同测试集数据进行分类。本发明采用多域特征提取的方法有效并最大程度上保证了提取的特征中所包含的信息量,将空间和频率两域的融合体现在了方法过程中,减少了多域特征融合后降维的步骤,提高了算法效率。
Description
技术领域
本发明涉及信号模式识别领域,尤其涉及一种基于EEG脑电信号数据的状态分类方法。
背景技术
当前随着以信息录入为主的信息化作业任务量的提高,装备或产品作业人员的功能定位由传统的机械操作类的任务角色逐步转换为以监督、决策为主的认知类任务角色,劳动密集型程度减少而带来的单个作业人员信息处理任务量和信息处理通道维度的增加,进而导致作业人员作业压力的增加。而作业压力增加很容易导致作业人员出现注意力资源不足,情景意识下降的问题,最终导致装备人机系统的绩效水平无法满足设计要求。
生理参数测量法通过对生理参数信号的采集、处理等过程实现对信号的分类,具有客观、实时、以及对被试者较少干扰的特点,在近期的研究中被越来越多的应用。在脑电信号、眼电指标、心电指标以及呼吸信号等诸多生理指标中,脑电信号因其的敏感性,较高的时间分辨率以及良好的适用性被广泛应用于对状态模式识别研究中。
在之前的研究中,存在基于脑电的功率谱密度(Power Spectral Density,PSD),使用支持向量机的方法对脑电信号进行分类的研究;此外,在基于脑电的乘员信息处理作业脑力负荷状态识别模型的研究中通过融合小波包分解(Wavelet PacketDecomposition,WPD)算法完成脑电参数的预处理,建立了负荷状态分类的脑电指标特征输入空间,随后通过粒子群优化算法(Particle Swarm Optimization,PSO)对支持向量机(Support Vector Machine,SVM)的参数进行寻优处理,最终构建了脑电信号状态识别模型。上述研究中,只使用了脑电数据的频域信息及特征,而在真实采集到的EEG(Electroencephalogram脑电信号)数据具有多导联空间时序信息特征,同时频域维度下存在良好的脑电信号分类信息特征,在此基础上,本专利提出保留脑电信号的多维度特性,提升对脑电信号分类的性能。
发明内容
为了克服现有技术的缺陷,本发明的目的在于通过提取脑电信号的多域特征,对脑电信号进行分类。本发明提供一种基于EEG脑电信号数据的状态分类方法,其包括以下步骤:
S1,获取不同分类的EEG脑电信号数据、以及相应的垂直和水平眼电信号数据;
获取两类EEG脑电信号数据、以及相应的垂直和水平眼电信号数据;通过从数据库中获取,或者通过测试仪中的电极获取,获取的EEG脑电信号数据都已知属于哪个分类;脑电信号数据由多通道的连续的采样数据组成,对每次获取的EEG脑电信号数据进行切片,构成一个样本,每个样本选用A个采样点作为样本的长度样本的步进长度为STEP个采样点,每个样本包含了多个通道的采样数据;
S2,EEG脑电信号数据预处理;
首先对脑电信号数据进行带通滤波,所述带通滤波频率为1-30Hz;
随后使用独立成分分析的方法对脑电信号数据中的眨眼和眼球运动伪迹进行去除;得到经预处理后的脑电信号数据;
S3,基于过程的特征提取,包括以下步骤:
S31,频域特征提取;
对于单个样本的单个通道根据快速傅里叶变换得到样本的功率谱密度;
将频率进行分段,得到B个频率段,根据不同频率段下的功率谱密度特征向量集合得到频域脑电信号数据矩阵Vfr,Vfr为N×Mfr维的矩阵,其中N为通道数量,Mfr是单个通道在fr频率段下的特征维度,fr为B个频率段中的一个,矩阵元素为相应频率段下离散频率点k处的功率谱密度值;
S32,空域特征提取;
S321,生成降维后的空间滤波器;
得到在fr频率段下降维后空间滤波器Wfr,c,矩阵Wfr,c为2m*N的矩阵,其中m为设定值,N为通道数量;
S322,生成单个样本的频-空特征向量;
将步骤S31获得的频域脑电信号数据矩阵Yfr与降维后空间滤波器Wfr,c相乘得到特征矩阵Yfr:
Yfr=Wfr,cVfr (19)
其中,特征矩阵Yfr为2m*Mfr的矩阵;
求取特征矩阵Yfr中的每行数据的方差,将特征矩阵Yfr的行数据替换为相应行数据的方差,特征矩阵Yfr变为维度为2m的特征向量;将所有频率段的特征向量直接连接组合后,成为单个样本的频-空特征向量,单个样本的频-空特征向量的维度为2m*B;
S4,训练分类模型;
重复步骤S3得到脑电信号数据的多个样本的频-空特征向量,使用频-空特征向量以及依据步骤S1获取的脑电信号数据分类对每个频-空特征向量设置相应的标签,采用支持向量机方法根据频-空特征向量和标签对分类模型进行训练;
优选的,步骤S31,频域特征提取步骤具体为:
脑电信号数据由多通道的连续的采样数据组成,对于单个样本选用1024个采样点作为样本的长度取值,从脑电信号数据中选取多个样本,样本的步进长度为512个采样点,对于单个样本的单个通道使用式(3)得到单个样本单个通道的功率谱密度:
其中,X(n)为有限长度的随机信号序列,n为随机信号序列X(n)长度,即时域信号点的数量,Pk为每个离散频率点k处的功率谱密度,Pk一起组成了随机信号序列X(n)的功率谱密度,N为快速傅里叶变换后得到的离散频率点的点数,N=n;FFT[X(n)]是长度为n的随机信号序列X(n)的快速傅里叶变换过程,FFT*[X(n)]为FFT[X(n)]的共轭表达式,fs为采样频率;
将频率进行分段,得到4个频率段,分别是delta为1-4Hz、theta为4-8Hz、alpha为8-12Hz和beta为13-30Hz;不同频率段下的功率谱密度特征向量集合表示为:
Pdelta={Pfreq},freq∈delta (4)
Ptheta={Pfreq},freq∈theta (5)
Palpha={Pfreq},freq∈alpha (6)
Pbeta={Pfreq},freq∈beta (7)
其中,Pfreq为计算得到的离散功率谱密度,Pdelta、Ptheta、Palpha、Pbeta为不同频率段下功率谱密度向量的集合,单个通道下4个频率段的特征维度分别为Mfr,Mfr的值为相应频率段fr下离散频率点k的数量;fr表示频率段,取值为集合{delta,theta,alpha,beta}中的一个;
根据不同频率段下的功率谱密度特征向量集合得到频域脑电信号数据矩阵Vfr,Vfr为N×Mfr维的矩阵,其中N为通道数量30,Mfr是单个通道在fr频率段下的特征维度,fr为4个频率段中的一个,矩阵元素为相应频率段下离散频率点k处的功率谱密度值。
优选的,步骤S321,生成降维后的空间滤波器步骤具体为:
对于同一个对象有两类脑电信号数据,如果脑电信号数据属于第一类,则将Vfr表示为Vfr,1,如果脑电信号数据属于第二类,则将Vfr表示为Vfr,2,分别得到同一个对象在fr频率段下的频域脑电信号数据矩阵Vfr,1和Vfr,2,在此基础上计算出两类脑电信号数据归一化后的协方差矩阵C1和C2:
式中,Vfr,1 T和Vfr,2 T分别用以表示Vfr,1和Vfr,2的转置,trace(V)用于表示矩阵的迹,将同一对象的每一类脑电信号数据在fr频率段下所有样本的归一化的协方差矩阵进行均值处理,之后通过矩阵求和的运算方法得到融合协方差矩阵C:
由于融合协方差矩阵C为正定矩阵,其转化成为特征向量和对角矩阵表达式的形式:
C=UΛUT (11)
其中:Λ为特征值组成的对角矩阵,同时特征值按照降序输入对角阵中;U为特征向量矩阵,对其进行白化转换表示为:
P=Λ-0.5UT (12)
S1=PC1PT (13)
S2=PC2PT (14)
对于S1和S2,它们具有公共的特征向量,且存在如下的条件,如果:
S1=HΛ1HT (15)
则有:
S2=HΛ2HT (16)
I=Λ1+Λ2 (17)
其中:S1和S2两矩阵的特征值之和为1,I为单位矩阵;因此对于特征向量的矩阵H,当对应一个类别S1具有最大的特征值时,另一类S2具有最小的特征值;将脑电信号数据白化变换为特征值最大的特征向量的方法为判别两类脑力信号的最优特征值向量;通过上述方法得到的投影矩阵,即空间滤波器W表示为:
W=HTP (18)
矩阵H的首尾行向量能够体现出两类信号的最大化差异信息,故提取空间滤波器W矩阵前后各m行构成在fr频率段下降维后空间滤波器Wfr,c,矩阵Wfr,c为2m*N的矩阵,m为设定值。
优选的,步骤S2中使用独立成分分析的方法对脑电信号数据中的眨眼和眼球运动伪迹进行去除,得到经预处理后的脑电信号数据;具体为:
独立成分分析方法的基础为假设某单一信号输入端收集到的数据是由n个独立成分混合得到的,使用数学公式表示为:
yj=wj1x1+wj2x2+…+wjnxn,j=1,2,3…,m (1)
式(1)转化为以下的形式:
YCOL=WMIXXINDEP (2)
该方法通过线性变换将多信号源融合信号转换为多个独立源成分信号的线性组合;式(2)即为独立成分分析假设的数学表达式,式中XINDEP=[x1,x2,…,xn]T表示统计学含义的n维独立源的向量表达式,而YCOL=[y1,y2,…,ym]T表示采集到的m维信号数据的向量表达式;WMIX为混合矩阵;上述独立源信号通过m*n的混合矩阵WMIX进行混合;在分析过程中独立信号源XINDEP和混合矩阵WMIX均是未知的,独立成分分析使用采集到m维信号的随机向量作为数据输入,通过优化迭代算法估计n维独立源信号XINDEP;
使用独立成分分析方法将一个脑电信号数据样本分解为与脑电信号数据通道数量一致的独立分量,根据计算分解后的各分量与水平及垂直眼电信号的相关性水平,确定与眼电信号高度相关的独立分量,并将其作为眼电伪迹进行去除;将XINDEP中判定为伪迹的独立源的向量内的值改为0,将修改后的X输入公式(2)得到重构后的m维脑电信号。
优选的,还包括步骤S5,对样本空间的不同测试集数据进行分类,具体步骤为:
对样本空间中的不同测试集数据进行分类时,执行步骤S2、S31后得到频域脑电信号数据矩阵Vfr,并直接乘以对象已有的Wfr,c,得到单个样本的频-空特征向量,将脑电信号数据的样本的频-空特征向量输入到训练好的分类模型中,对脑电信号数据进行分类。
与现有技术相比,本发明具有以下有益效果:
(1)针对脑电信号的单域特征提取特征信息量有限、分类性能较差的问题,本发明采用多域特征提取的方法有效并最大程度上保证了提取的特征中所包含的信息量,相较于单域特征提取的方法,提高了脑电信号分类的性能。
(2)针对脑电信号典型多域特征提取因其是分别提取不同域特征在后期组合构建的特征矩阵,而造成维度爆炸需要后续降维的问题。本发明将空间和频率两域的融合体现在了方法过程中,减少了多域特征融合后降维的步骤,大大提高了算法效率。
(3)针对多模态信号组合特征提取导致的多模态信号采集难度大、信号同步要求高的问题,在降低操作难度、时间复杂度的基础上,保证了脑电信号数据分类的性能。
附图说明
图1是本发明的流程结构示意图;
图2是本实施例中可穿戴式脑电采集电极点具体位置。
具体实施方式
本专利的附图用来进一步理解本发明的主要技术内容,同时作为本发明的框架性指导部分,与本发明的具体实施例文字部分结合,一起解释本发明,但并不够成对本发明的限制。
本实施例公开了一种基于EEG脑电信号数据的状态分类方法如图1所示,具体包括以下步骤:
S1,获取不同分类的EEG脑电信号数据、以及相应的垂直和水平眼电信号数据;
不同分类的EEG脑电信号数据、以及相应的垂直和水平眼电信号数据可以通过从数据库中获取,或者通过测试仪中的电极获取。获取的EEG脑电信号数据都已知属于哪个分类。
本实施例中将EEG脑电信号数据分为两类,对于每个对象分别获取属于不同分类的EEG脑电信号数据。脑电信号数据由多通道的连续的采样数据组成,对每次获取的EEG脑电信号数据进行切片,构成一个样本,每个样本选用1024个采样点作为样本的长度样本的步进长度为512个采样点,每个样本包含了多个通道的采样数据。
当从电极获取脑电信号时,结合现行装备或产品中多任务并行处理的任务属性特征,通过设置不同工作记忆任务,利用脑电信号测试仪中的电极获取脑电信号数据、以及相应的垂直和水平眼电信号。本实施例中不同的工作记忆任务就是指任务对记忆要求的难度不同,如输入上1秒看到的数字值和位置获取的脑电信号数据为第一类,输入10秒前看到的数字值和位置获取的脑电信号数据为第二类。脑电信号电极点位置使用“10-20国际标准导联系统”标准,具体导联位置如图2所示。
本步骤用于获取原始脑电信号数据和原始垂直和水平眼电信号数据,原始脑电信号数据为30通道,垂直和水平眼电信号用于对脑电信号的伪迹进行去除。
S2,EEG脑电信号数据预处理;
由于脑电EEG信号本质上为头皮表面各电极点的电位差,而外部噪声、头部或眼球运动均会导致采集到的脑电信号数据不准确。因此,需要对脑电信号数据进行伪迹的去除。首先对脑电信号数据进行带通滤波,本发明选择带通滤波频率为1-30Hz,该频率范围在之前的研究中被证明与脑部活动有较强的相关性,而带通滤波的方法可以去除属于高频段的伪迹信号。
随后使用现有的独立成分分析(Independent components analysis,ICA)的方法对脑电信号数据中的眨眼和眼球运动伪迹进行去除。独立成分分析隶属于盲源信号分离,其分析基础为假设某单一信号输入端收集到的数据是由n个独立成分混合得到的,使用数学公式表示为:
yj=wj1x1+wj2x2+…+wjnxn,j=1,2,3…,m (1)
式(1)可以转化为以下的形式:
YCOL=WMIXXINDEP (2)
该方法通过线性变换将多信号源融合信号转换为多个独立源成分信号的线性组合。式(2)即为独立成分分析假设的数学表达式,式中XINDEP=[x1,x2,…,xn]T表示统计学含义的n维独立源的向量表达式,而YCOL=[y1,y2,…,ym]T表示采集到的m维信号数据的向量表达式。WMIX为混合矩阵。上述独立源信号通过m*n的混合矩阵WMIX进行混合。在分析过程中独立信号源XINDEP和混合矩阵WMIX均是未知的,独立成分分析使用采集到m维信号的随机向量作为数据输入,通过优化迭代算法估计n维独立源信号XINDEP。
使用独立成分分析方法将一个脑电信号数据样本分解为与脑电信号数据通道数量一致的独立分量,本实施例中通道数量为30,因此将脑电信号数据样本分解为30个独立分量。之所以选取该数量的独立分量,主要考虑以下两个因素,首先,独立成分分析输出的独立成分维度不能超过输入的维度,即获取的脑电信号数据通道个数。其次,相关研究表明独立成分输出的维度越高,最终分类的准确率也就越高。
在完成独立成分分析后,根据计算分解后的各分量与水平及垂直眼电信号的相关性水平,确定与眼电信号高度相关的独立分量,并将其作为眼电伪迹进行去除。将XINDEP中判定为伪迹的独立源的向量内的值改为0,将修改后的X输入公式(2)得到重构后的m维脑电信号。
本步骤得到的结果为经过带通滤波及去除典型伪迹信号后提出的脑电信号数据,被称为经预处理后的脑电信号数据。
S3,基于过程的特征提取;
本步骤输入的数据为经预处理后的脑电信号数据。
本实施例使用基于过程的脑电信号数据特征提取方法,总体而言是将频域特征提取结果作为空域滤波特征提取算法的输入。进一步的,对于频域特征提取算法,本实施例选取脑电信号数据不同频率段的能量特征作为特征输入;而对于空域滤波提取算法,本实施例提出了基于频域特征的共空间模式(Based on frequency domain features-CommonSpatial Pattern,BFF-CSP)方法构建特征空间。下文将对所构建的特征提取方法进行详细说明。
S31,频域特征提取;
本发明采用直接法求取功率谱密度这一脑电信号频域能量特征,这种方法是基于快速傅里叶变换而实现的,随后提取不同功率下的功率谱密度。
该方法的理论基础为瑞利能量定理,该定理假设对于函数平方的积分等于傅里叶变换后的平方的积分。因此对于有限长度的随机信号序列X(n),通过快速傅里叶变换得到每个离散频率点k处的功率谱密度的表达式为:
式中,Pk为每个离散频率点k处的功率谱密度,Pk一起组成了随机信号序列X(n)的功率谱密度,n为随机信号序列X(n)长度,即时域信号点的数量,N为快速傅里叶变换后得到的离散频率点的点数,N=n;FFT[X(n)]是长度为n的随机信号序列X(n)的快速傅里叶变换过程,FFT*[X(n)]为FFT[X(n)]的共轭表达式,fs为采样频率。
在本发明中,脑电信号数据由多通道的连续的采样数据组成,对每次获取的EEG脑电信号数据进行切片,构成一个样本,每个样本选用1024个采样点作为样本的长度,样本的步进长度为512个采样点,每个样本包含了多个通道的采样数据。对于单个样本的单个通道使用式(3)得到单个样本单个通道的的功率谱密度,并根据脑电信号的频率分布情况提取了以下四类频率指标,分别是delta为1-4Hz、theta为4-8Hz、alpha为8-12Hz和beta为13-30Hz。不同频率段下的功率谱密度特征向量集合可以表示为:
Pdelta={Preq},freq∈delta (4)
Ptheta={Pfreq},freq∈theta (5)
Palpha={Pfreq},freq∈alpha (6)
Pbeta={Pfreq},freq∈beta (7)
上式中:Pfreq为计算得到的离散功率谱密度,Pdelta、Ptheta、Palpha、Pbeta为不同频率段下功率谱密度向量的集合,单个通道下四个频率段的特征维度分别为Mfr,值为相应频率段fr下离散频率点k的数量;fr表示频率段,取值为集合{delta,theta,alpha,beta}中的一个;所有通道中某个确定频率段fr的特征维度值都是相同的。
因此本步骤得到按频域特征提取后形成的样本空间,被称为频域脑电信号数据矩阵Vfr,Vfr为N×Mfr维的矩阵,其中N为通道数量,Mfr是单个通道在fr频率段下的特征维度,fr为集合{delta,theta,alpha,beta}中的一个,矩阵元素为相应频率段下离散频率点k处的功率谱密度值。
由于本实施例中采用的脑电信号数据的通道维度N为30个,因此样本对应于每一个频率段的离散功率谱特征向量维度分别为:Mdelta*30、Mdelta*30、Malpha*30、Mbeta*30。
S32,空域特征提取;
基础的共空间模式是一类时空序列信号的中的空域滤波特征的提取方法,该方法可以从多通道采集的高维时序脑电信号数据中得到特定的空间分布组成。共空间模式算法的基础思想是利用矩阵的对角化,获取一组最优空间滤波器进行投影,在经过最优空间滤波器的处理后可以最大化两类信号之间的差异,进而取得具有较高差异化的特征空间。原始的空间从时域信号中构建最优空间滤波器。本发明将上一步骤输出的频域特征样本空间按照频率段分别作为本步骤的输入,将频-空域特征融合体现在算法的过程中而不是后期的特征组合。
S321,生成降维后的空间滤波器;
对于同一个对象有两类脑电信号数据,如果脑电信号数据属于第一类,则将相应Vfr表示为Vfr,1,如果脑电信号数据属于第二类,则将相应Vfr表示为Vfr,2,分别得到同一个对象在fr频率段下的频域脑电信号数据矩阵Vfr,1和Vfr,2,在此基础上可以计算出两类脑电信号数据归一化后的协方差矩阵C1和C2。
式中,Vfr,1 T和Vfr,2 T分别用以表示Vfr,1和Vfr,2的转置,trace(V)用于表示矩阵的迹。将同一对象的每一类脑电信号数据在fr频率段下所有样本的归一化的协方差矩阵进行均值处理,之后通过矩阵求和的运算方法得到融合协方差矩阵C:
由于融合协方差矩阵C为正定矩阵,其可以转化成为特征向量和对角矩阵表达式的形式:
C=UΛUT (11)
其中:Λ为特征值组成的对角矩阵,同时特征值按照降序输入对角阵中。U为特征向量矩阵,对其进行白化转换可表示为:
P=Λ-0.5UT (12)
S1=PC1PT (13)
S2=PC2PT (14)
对于S1和S2,它们具有公共的特征向量,且存在如下的条件,如果:
S1=HΛ1HT (15)
则有:
S2=HΛ2HT (16)
I=Λ1+Λ2 (17)
其中:S1和S2两矩阵的特征值之和为1,I为单位矩阵。因此对于特征向量的矩阵H,当对应一个类别S1具有最大的特征值时,另一类S2具有最小的特征值。因此将脑电信号数据白化变换为特征值最大的特征向量的方法为判别两类脑力信号的最优特征值向量。通过上述方法得到的投影矩阵,即空间滤波器W可以表示为:
W=HTP (18)
随后进行特征提取,将降维的思路融入其中,由于前述步骤中H的首尾向量能够体现出两类信号的最大化差异信息,故提取空间滤波器W矩阵前后各m行构成在fr频率段下降维后空间滤波器Wfr,c,矩阵Wfr,c为2m*N的矩阵。
S322,生成单个样本的频-空特征向量;
将步骤S31获得的频域脑电信号数据矩阵Vfr与降维后空间滤波器Wfr,c相乘得到特征矩阵Yfr:
Yfr=Wfr,cVfr (19)
此时,特征矩阵Yfr为2m*Mfr的矩阵。
求取特征矩阵Yfr中的每行数据的方差,并进行对数化处理。经过本步骤,特征矩阵Yfr变为维度为2m的特征向量,即单个样本在fs频率段的特征向量维度为2m,将所有频率段的特征向量直接连接组合后,输出单个样本的频-空特征向量。本实施例中为四个频率段,因此单个样本的频-空特征向量的维度为2m*4。
S4,训练分类模型;
S41,特征空间的预处理;
重复步骤S3可得到多个样本的频-空特征向量,使用频-空特征向量以及依据步骤S1获取的脑电信号数据分类对每个频-空特征向量设置相应的标签,频-空特征向量所对应的原始EEG脑电信号数据为第一分类的标签设定为1,为第二分类的标签设定为-1。
S42,判别分类模型;
采用现有支持向量机方法根据频-空特征向量和标签对分类模型进行训练。
基础的支持向量机方法作为一种典型的判别分类方法在分类问题中被广泛使用,但是在脑电信号数据的判别分类问题中,样本数据往往是线性不可分的,本实施例中使用基于核函数的支持向量机方法来解决这一问题。
支持向量机的算法分类性能和具体泛化能力直接取决于核函数以及松弛变量的引入。松弛变量表征了离散点离散的距离,为解决这类离群点在分类过程误分类的问题,引入惩罚因子进行调节,惩罚因子的取值越大,离散点对于模型分类过程中的重要度越高。另外对于非线性样本数据,支持向量机可引入核函数将输入特征转换入更高维度的特征空间内,实现输入数据的线性可分性。在众多的核函数中,高斯核函数(Radial BasisFunction,RBF)具有判别性能良好,映射范围较大的特点,因此本专利引入高斯核函数方法构建模型。对于惩罚因子C和高斯核函数参数gamma的参数优化过程中选用网格寻优和交叉验证结合的方法,来获取最优化的参数和结果。
具体而言,惩罚因子C的取值范围是[0.001,0.01,0.1,1,100],高斯核函数参数gamma的取值范围是[0.1,1,10,100,1000]。
本实施例可选用4折交叉验证来实现参数的筛选和优化选择。
本实施例使用最优的参数和模型评估结果作为输出。
S5,对样本空间的测试集数据进行分类;
对样本空间的测试集数据进行分类时,只需要执行步骤S2、S31后得到频域脑电信号数据矩阵Vfr,并直接乘以现有的Wfr,c即可得到单个样本的频-空特征向量,将待测试的脑电信号数据的样本的频-空特征向量输入到训练好的分类模型中,即可对脑电信号数据进行分类。
本专利提出的方法有效地解决了单域特征提取特征信息量有限性能较差、多域特征提取造成维度爆炸需要后续降维操作难度较大的问题,在降低操作难度、时间复杂度的基础上,保证了分类的性能,这对后续针对装备或产品的对象压力状态分类结果实现人机功能再分配、人机系统自适应自动化提供了方法层面的支撑和依据。
最后应说明的是:以上所述的各实施例仅用于说明本发明的技术方案,而非对其限制;尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述实施例所记载的技术方案进行修改,或者对其中部分或全部技术特征进行等同替换;而这些修改或替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。
Claims (5)
1.一种基于EEG脑电信号数据的状态分类方法,其特征在于:其包括以下步骤:
S1,获取不同分类的EEG脑电信号数据、以及相应的垂直和水平眼电信号数据;
获取两类EEG脑电信号数据、以及相应的垂直和水平眼电信号数据;通过从数据库中获取,或者通过测试仪中的电极获取,获取的EEG脑电信号数据都已知属于哪个分类;获取的脑电信号数据由多通道的连续的采样数据组成,对每次获取的EEG脑电信号数据进行切片,构成一个样本,每个样本选用A个采样点作为样本的长度,样本的步进长度为STEP个采样点,每个样本包含了多个通道的采样数据;
S2,EEG脑电信号数据预处理;
首先对脑电信号数据进行带通滤波,所述带通滤波频率为1-30Hz;
随后使用独立成分分析的方法对脑电信号数据中的眨眼和眼球运动伪迹进行去除;得到经预处理后的样本;
S3,基于过程的特征提取,包括以下步骤:
S31,频域特征提取;
对于单个预处理后的样本的单个通道根据快速傅里叶变换得到功率谱密度;
将频率进行分段,得到B个频率段,根据不同频率段下的功率谱密度特征向量集合得到频域脑电信号数据矩阵Vfr,Vfr为N×Mfr维的矩阵,其中N为通道数量,Mfr是单个通道在fr频率段下的特征维度,fr为B个频率段中的一个,矩阵元素为相应频率段下离散频率点k处的功率谱密度值;
S32,空域特征提取;
S321,生成降维后的空间滤波器;
得到在fr频率段下降维后空间滤波器Wfr,c,矩阵Wfr,c为2m*N的矩阵,其中m为设定值,N为通道数量;
S322,生成单个样本的频-空特征向量;
将步骤S31获得的频域脑电信号数据矩阵Vfr与降维后空间滤波器Wfr,c相乘得到特征矩阵Yfr:
Yfr=Wfr,cVfr (19)
其中,特征矩阵Yfr为2m*Mfr的矩阵;
求取特征矩阵Yfr中的每行数据的方差,将特征矩阵Yfr的行数据替换为相应行数据的方差,特征矩阵Yfr变为维度为2m的特征向量;将所有频率段的特征向量直接连接组合后,成为单个样本的频-空特征向量,单个样本的频-空特征向量的维度为2m*B;
S4,训练分类模型;
重复步骤S3得到脑电信号数据的多个样本的频-空特征向量,使用频-空特征向量以及依据步骤S1获取的脑电信号数据分类对每个频-空特征向量设置相应的标签,采用支持向量机方法根据频-空特征向量和标签对分类模型进行训练。
2.根据权利要求1所述的基于EEG脑电信号数据的状态分类方法,其特征在于:所述步骤S31,频域特征提取步骤具体为:
脑电信号数据由多通道的连续的采样数据组成,对于单个样本选用1024个采样点作为样本的长度取值,从脑电信号数据中选取多个样本,样本的步进长度为512个采样点,对于单个样本的单个通道使用式(3)得到单个样本单个通道的功率谱密度:
其中,X(n)为有限长度的随机信号序列,n为随机信号序列X(n)长度,即时域信号点的数量,Pk为每个离散频率点k处的功率谱密度,Pk一起组成了随机信号序列X(n)的功率谱密度,N为快速傅里叶变换后得到的离散频率点的点数,N=n;FFT[X(n)]是长度为n的随机信号序列X(n)的快速傅里叶变换过程,FFT*[X(n)]为FFT[X(n)]的共轭表达式,fs为采样频率;
将频率进行分段,得到4个频率段,分别是delta为1-4Hz、theta为4-8Hz、alpha为8-12Hz和beta为13-30Hz;不同频率段下的功率谱密度特征向量集合表示为:
Pdelta={Pfreq},freq∈delta (4)
Ptheta={Pfreq},freq∈theta (5)
Palpha={Pfreq},freq∈alpha (6)
Pbeta={Pfreq},freq∈beta (7)
其中,Pfreq为计算得到的离散功率谱密度,Pdelta、Ptheta、Palpha、Pbeta为不同频率段下功率谱密度向量的集合,单个通道下4个频率段的特征维度分别为Mfr,Mfr的值为相应频率段fr下离散频率点k的数量;fr表示频率段,取值为集合{delta,theta,alpha,beta}中的一个;
根据不同频率段下的功率谱密度特征向量集合得到频域脑电信号数据矩阵Vfr,Vfr为N×Mfr维的矩阵,其中N为通道数量30,Mfr是单个通道在fr频率段下的特征维度,fr为4个频率段中的一个,矩阵元素为相应频率段下离散频率点k处的功率谱密度值。
3.根据权利要求1所述的基于EEG脑电信号数据的状态分类方法,其特征在于:所述步骤S321,生成降维后的空间滤波器步骤具体为:
对于同一个对象有两类脑电信号数据,如果脑电信号数据属于第一类,则将Vfr表示为Vfr,1,如果脑电信号数据属于第二类,则将Vfr表示为Vfr,2,分别得到同一个对象在fr频率段下的频域脑电信号数据矩阵Vfr,1和Vfr,2,在此基础上计算出两类脑电信号数据归一化后的协方差矩阵C1和C2:
式中,Vfr,1 T和Vfr,2 T分别用以表示Vfr,1和Vfr,2的转置,trace(V)用于表示矩阵的迹,将同一对象的每一类脑电信号数据在fr频率段下所有样本的归一化的协方差矩阵进行均值处理,之后通过矩阵求和的运算方法得到融合协方差矩阵C:
由于融合协方差矩阵C为正定矩阵,其转化成为特征向量和对角矩阵表达式的形式:
C=UΛUT (11)
其中:Λ为特征值组成的对角矩阵,同时特征值按照降序输入对角阵中;U为特征向量矩阵,对其进行白化转换表示为:
P=Λ-0.5UT (12)
S1=PC1PT (13)
S2=PC2PT (14)
对于S1和S2,它们具有公共的特征向量,且存在如下的条件,如果:
S1=HΛ1HT (15)
则有:
S2=HΛ2HT (16)
I=Λ1+Λ2 (17)
其中:S1和S2两矩阵的特征值之和为1,I为单位矩阵;因此对于特征向量的矩阵H,当对应一个类别S1具有最大的特征值时,另一类S2具有最小的特征值;将脑电信号数据白化变换为特征值最大的特征向量的方法为判别两类脑力信号的最优特征值向量;通过上述方法得到的投影矩阵,即空间滤波器W表示为:
W=HTP (18)
矩阵片的首尾行向量能够体现出两类信号的最大化差异信息,故提取空间滤波器W矩阵前后各m行构成在fr频率段下降维后空间滤波器Wfr,c,矩阵Wfr,c为2m*N的矩阵,m为设定值。
4.根据权利要求1所述的基于EEG脑电信号数据的状态分类方法,其特征在于:所述步骤S2中使用独立成分分析的方法对脑电信号数据中的眨眼和眼球运动伪迹进行去除,得到经预处理后的脑电信号数据;具体为:
独立成分分析方法的基础为假设某单一信号输入端收集到的数据是由n个独立成分混合得到的,使用数学公式表示为:
yj=Wj1x1+wj2x2+…+wjnxn,j=1,2,3…,m (1)
式(1)转化为以下的形式:
YCOL=WMIXXINDEP (2)
该方法通过线性变换将多信号源融合信号转换为多个独立源成分信号的线性组合;式(2)即为独立成分分析假设的数学表达式,式中XINDEP=[x1,x2,...,xn]T表示统计学含义的n维独立源的向量表达式,而YCOL=[y1,y2,...,ym]T表示采集到的m维信号数据的向量表达式;WMIX为混合矩阵;上述独立源信号通过m*n的混合矩阵WMIX进行混合;在分析过程中独立信号源XINDEP和混合矩阵WMIX均是未知的,独立成分分析使用采集到m维信号的随机向量作为数据输入,通过优化迭代算法估计n维独立源信号XINDEP;
使用独立成分分析方法将一个脑电信号数据样本分解为与脑电信号数据通道数量一致的独立分量,根据计算分解后的各分量与水平及垂直眼电信号的相关性水平,确定与眼电信号高度相关的独立分量,并将其作为眼电伪迹进行去除;将XINDEP中判定为伪迹的独立源的向量内的值改为0,将修改后的X输入公式(2)得到重构后的m维脑电信号。
5.根据权利要求1所述的基于EEG脑电信号数据的状态分类方法,其特征在于:还包括步骤S5,对样本空间的不同测试集数据进行分类,具体步骤为:
对样本空间中的不同测试集数据进行分类时,执行步骤S2、S31后得到频域脑电信号数据矩阵Vfr,并直接乘以对象已有的Wfr,c,得到单个样本的频-空特征向量,将脑电信号数据的样本的频-空特征向量输入到训练好的分类模型中,对脑电信号数据进行分类。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111502840.4A CN114218986B (zh) | 2021-12-10 | 2021-12-10 | 基于eeg脑电信号数据的状态分类方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111502840.4A CN114218986B (zh) | 2021-12-10 | 2021-12-10 | 基于eeg脑电信号数据的状态分类方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114218986A true CN114218986A (zh) | 2022-03-22 |
CN114218986B CN114218986B (zh) | 2024-05-07 |
Family
ID=80700657
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111502840.4A Active CN114218986B (zh) | 2021-12-10 | 2021-12-10 | 基于eeg脑电信号数据的状态分类方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114218986B (zh) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115227263A (zh) * | 2022-07-15 | 2022-10-25 | 无锡智蜂科技有限公司 | 基于EEG信号theta振荡进行调控的神经刺激系统 |
CN117958841A (zh) * | 2024-03-28 | 2024-05-03 | 北京理工大学 | 一种脑电信号头动伪影去除方法、装置、存储介质及产品 |
CN118220195A (zh) * | 2024-05-24 | 2024-06-21 | 广汽埃安新能源汽车股份有限公司 | 一种基于脑电的车辆控制方法及装置 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107239142A (zh) * | 2017-06-01 | 2017-10-10 | 南京邮电大学 | 一种结合公共空间模式算法和emd的脑电信号特征提取方法 |
WO2018014436A1 (zh) * | 2016-07-18 | 2018-01-25 | 天津大学 | 一种提高情绪识别模型时间鲁棒性的情绪脑电识别方法 |
CN108042132A (zh) * | 2017-12-27 | 2018-05-18 | 南京邮电大学 | 基于dwt和emd融合csp的脑电特征提取方法 |
CN110889082A (zh) * | 2019-12-03 | 2020-03-17 | 中国航空综合技术研究所 | 基于系统工程理论的人机工程设备综合评价方法 |
WO2020143300A1 (zh) * | 2019-01-07 | 2020-07-16 | 哈尔滨工业大学(深圳) | 听觉注意状态觉醒度识别方法、装置及存储介质 |
CN112515685A (zh) * | 2020-11-10 | 2021-03-19 | 上海大学 | 基于时频共融的多通道脑电信号通道选择方法 |
CN113158793A (zh) * | 2021-03-15 | 2021-07-23 | 东北电力大学 | 一种基于多特征融合的多类运动想象脑电信号识别方法 |
-
2021
- 2021-12-10 CN CN202111502840.4A patent/CN114218986B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2018014436A1 (zh) * | 2016-07-18 | 2018-01-25 | 天津大学 | 一种提高情绪识别模型时间鲁棒性的情绪脑电识别方法 |
CN107239142A (zh) * | 2017-06-01 | 2017-10-10 | 南京邮电大学 | 一种结合公共空间模式算法和emd的脑电信号特征提取方法 |
CN108042132A (zh) * | 2017-12-27 | 2018-05-18 | 南京邮电大学 | 基于dwt和emd融合csp的脑电特征提取方法 |
WO2020143300A1 (zh) * | 2019-01-07 | 2020-07-16 | 哈尔滨工业大学(深圳) | 听觉注意状态觉醒度识别方法、装置及存储介质 |
CN110889082A (zh) * | 2019-12-03 | 2020-03-17 | 中国航空综合技术研究所 | 基于系统工程理论的人机工程设备综合评价方法 |
CN112515685A (zh) * | 2020-11-10 | 2021-03-19 | 上海大学 | 基于时频共融的多通道脑电信号通道选择方法 |
CN113158793A (zh) * | 2021-03-15 | 2021-07-23 | 东北电力大学 | 一种基于多特征融合的多类运动想象脑电信号识别方法 |
Non-Patent Citations (3)
Title |
---|
姜月;邹任玲;: "基于多特征融合的运动想象脑电信号识别研究", 中国医学物理学杂志, no. 05, 25 May 2019 (2019-05-25) * |
安凯等: "装备维修人机工程设计要求标准设想", 航空标准化与质量, no. 01, 31 January 2021 (2021-01-31), pages 3 - 5 * |
李明爱;崔燕;杨金福;郝冬梅;: "基于HHT和CSSD的多域融合自适应脑电特征提取方法", 电子学报, no. 12, 15 December 2013 (2013-12-15) * |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115227263A (zh) * | 2022-07-15 | 2022-10-25 | 无锡智蜂科技有限公司 | 基于EEG信号theta振荡进行调控的神经刺激系统 |
CN117958841A (zh) * | 2024-03-28 | 2024-05-03 | 北京理工大学 | 一种脑电信号头动伪影去除方法、装置、存储介质及产品 |
CN118220195A (zh) * | 2024-05-24 | 2024-06-21 | 广汽埃安新能源汽车股份有限公司 | 一种基于脑电的车辆控制方法及装置 |
Also Published As
Publication number | Publication date |
---|---|
CN114218986B (zh) | 2024-05-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114218986B (zh) | 基于eeg脑电信号数据的状态分类方法 | |
CN111329474B (zh) | 基于深度学习的脑电身份识别方法、系统及信息更新方法 | |
Ranjan et al. | Ocular artifact elimination from electroencephalography signals: A systematic review | |
CN111616682B (zh) | 基于便携式脑电采集设备的癫痫发作预警系统及应用 | |
Mahajan et al. | Classification of EEG using PCA, ICA and Neural Network | |
CN110781945A (zh) | 一种融合多种特征的脑电信号情感识别方法及系统 | |
CN104586387A (zh) | 一种时、频、空域多参数脑电特征提取与融合方法 | |
CN114533086B (zh) | 一种基于空域特征时频变换的运动想象脑电解码方法 | |
CN109480834A (zh) | 一种基于快速多维经验模态分解的脑电信号分类方法 | |
CN110135285B (zh) | 一种使用单导设备的脑电静息态身份认证方法及装置 | |
CN109657646B (zh) | 生理时间序列的特征表示与提取方法、装置及存储介质 | |
CN111184509A (zh) | 一种基于传递熵的情绪诱导脑电信号分类方法 | |
CN113536882B (zh) | 一种多类运动想象脑电信号特征提取及分类方法 | |
CN111310656A (zh) | 基于多线性主成分分析的单次运动想象脑电信号识别方法 | |
Mousa et al. | A novel brain computer interface based on principle component analysis | |
CN113208613A (zh) | 基于fhls特征选择的多模态bci时序优化方法 | |
CN114081503A (zh) | 去除脑电信号中眼电伪迹的方法 | |
CN113180659A (zh) | 一种基于三维特征和空洞全卷积网络的脑电情感识别系统 | |
Nasehi et al. | A novel effective feature selection algorithm based on S-PCA and wavelet transform features in EEG signal classification | |
CN108470182B (zh) | 一种用于非对称脑电特征增强与识别的脑-机接口方法 | |
CN111067513B (zh) | 一种特征权重自学习的睡眠质量检测关键脑区判定方法 | |
CN115414051A (zh) | 一种脑电信号自适应窗口的情绪分类识别方法 | |
Jiang et al. | Picture-induced EEG signal classification based on CVC emotion recognition system | |
Hsu | Wavelet-coherence features for motor imagery EEG analysis posterior to EOG noise elimination | |
Helal et al. | Using autoencoders for feature enhancement in motor imagery Brain-Computer Interfaces |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |