CN116662742A - 基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法 - Google Patents
基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法 Download PDFInfo
- Publication number
- CN116662742A CN116662742A CN202310776366.7A CN202310776366A CN116662742A CN 116662742 A CN116662742 A CN 116662742A CN 202310776366 A CN202310776366 A CN 202310776366A CN 116662742 A CN116662742 A CN 116662742A
- Authority
- CN
- China
- Prior art keywords
- frequency
- mask
- signal
- layer
- state
- 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
- 210000004556 brain Anatomy 0.000 title claims abstract description 43
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 41
- 238000000034 method Methods 0.000 title claims abstract description 36
- 238000005868 electrolysis reaction Methods 0.000 title claims abstract description 18
- 102100030511 Stanniocalcin-1 Human genes 0.000 claims abstract description 31
- 101710183756 Stanniocalcin Proteins 0.000 claims abstract description 29
- 238000013127 simulated treatment comparison Methods 0.000 claims abstract description 29
- 238000007781 pre-processing Methods 0.000 claims abstract description 7
- 102100024210 CD166 antigen Human genes 0.000 claims description 18
- 101000980840 Homo sapiens CD166 antigen Proteins 0.000 claims description 18
- 238000004422 calculation algorithm Methods 0.000 claims description 16
- 230000008569 process Effects 0.000 claims description 11
- 238000001914 filtration Methods 0.000 claims description 9
- 238000010276 construction Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 230000009466 transformation Effects 0.000 claims description 2
- 230000000694 effects Effects 0.000 description 12
- 230000005021 gait Effects 0.000 description 11
- 238000004458 analytical method Methods 0.000 description 10
- 230000006870 function Effects 0.000 description 9
- 230000033001 locomotion Effects 0.000 description 8
- 238000004364 calculation method Methods 0.000 description 7
- 206010017577 Gait disturbance Diseases 0.000 description 5
- 230000001537 neural effect Effects 0.000 description 5
- 210000005036 nerve Anatomy 0.000 description 4
- 230000009471 action Effects 0.000 description 3
- 238000001514 detection method Methods 0.000 description 3
- 210000003141 lower extremity Anatomy 0.000 description 3
- 101000701440 Homo sapiens Stanniocalcin-1 Proteins 0.000 description 2
- 101000701446 Homo sapiens Stanniocalcin-2 Proteins 0.000 description 2
- 208000018737 Parkinson disease Diseases 0.000 description 2
- 208000027089 Parkinsonian disease Diseases 0.000 description 2
- 206010034010 Parkinsonism Diseases 0.000 description 2
- 102100030510 Stanniocalcin-2 Human genes 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 210000000653 nervous system Anatomy 0.000 description 2
- 230000001575 pathological effect Effects 0.000 description 2
- 210000004281 subthalamic nucleus Anatomy 0.000 description 2
- 208000024891 symptom Diseases 0.000 description 2
- 210000001364 upper extremity Anatomy 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 208000012661 Dyskinesia Diseases 0.000 description 1
- 206010033799 Paralysis Diseases 0.000 description 1
- 208000006011 Stroke Diseases 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000004913 activation Effects 0.000 description 1
- 238000002266 amputation Methods 0.000 description 1
- 210000004227 basal ganglia Anatomy 0.000 description 1
- 230000006399 behavior Effects 0.000 description 1
- 239000000090 biomarker Substances 0.000 description 1
- 230000007177 brain activity Effects 0.000 description 1
- 230000003925 brain function Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 230000019771 cognition Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000006735 deficit Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 230000008451 emotion Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 210000003414 extremity Anatomy 0.000 description 1
- 239000012634 fragment Substances 0.000 description 1
- 230000005714 functional activity Effects 0.000 description 1
- 230000002267 hypothalamic effect Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000008449 language Effects 0.000 description 1
- 239000012528 membrane Substances 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 201000006417 multiple sclerosis Diseases 0.000 description 1
- 210000003205 muscle Anatomy 0.000 description 1
- 230000004770 neurodegeneration Effects 0.000 description 1
- 208000015122 neurodegenerative disease Diseases 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 238000010827 pathological analysis Methods 0.000 description 1
- 238000003909 pattern recognition Methods 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000011179 visual inspection Methods 0.000 description 1
- 210000000707 wrist Anatomy 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/369—Electroencephalography [EEG]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/369—Electroencephalography [EEG]
- A61B5/372—Analysis of electroencephalograms
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/7246—Details of waveform analysis using correlation, e.g. template matching or determination of similarity
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/10—Pre-processing; Data cleansing
-
- 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/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2123/00—Data types
- G06F2123/02—Data types in the time domain, e.g. time-series data
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
-
- 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
Landscapes
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Animal Behavior & Ethology (AREA)
- Theoretical Computer Science (AREA)
- Veterinary Medicine (AREA)
- General Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Heart & Thoracic Surgery (AREA)
- Surgery (AREA)
- Public Health (AREA)
- Biomedical Technology (AREA)
- Molecular Biology (AREA)
- Psychiatry (AREA)
- Medical Informatics (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Data Mining & Analysis (AREA)
- Evolutionary Biology (AREA)
- Signal Processing (AREA)
- Physiology (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Psychology (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开了基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，能够将隐马尔可夫模型与掩膜经验模态分解结合。具体方案为：采集脑电信号，并进行预处理。预处理后的脑电信号一方面进行掩膜经验模态分解、提取所需频率范围对应的本征模函数IMF分量，并计算所有IMF分量的希尔伯特包络；另一方面输入隐马尔科夫模型，通过设定状态的数量，得到每个状态的状态时间序列STC。将所有IMF分量的希尔伯特包络，与状态时间序列STC计算互相关系数，根据最大正相关系数进行STC与脑电频带的对应，如果不存在正相关系数则归于背景频带，最终得到五个频带。分别计算不同频带下分配的STC的时域特征，提取不同频带下状态占用率FO特征和状态生命时间LT特征。
Description
技术领域
本发明涉及医学信号处理技术领域，具体涉及基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法。
背景技术
大脑回路及内部网络能够产生复杂、丰富的动力状态，这些不断变化的状态表现出的功能激活和协同模式，反映了神经***在时间和空间上的特定的连接及相互作用。解码与大脑状态相关的信息，尤其是基于任务态的神经解码，对脑认知和脑功能网络的机制探索具有深远意义，如语言、情绪、运动等行为下的脑状态分析，大脑网络的时变重复模式识别和动态组织规律探索，上肢的手指、腕部、肘部等动作的分类以及下肢步态运动的周期性相位识别等。
运动是生命活动最基本的功能之一，运动的执行需要神经***对肢体和躯干各肌群进行精巧的调控来实现，因此针对特定动作下的大脑活动状态进行解码，其重要性不言而喻。相较于大量关于上肢运动意图的识别和分析，下肢运动，尤其是对于步态障碍的神经解码研究较少。目前，步态相关的神经控制机制尚不明确，难以进行有效的症状干预和治疗，步态调控也成为神经工程和医疗康复领域的一个关键挑战。对于下肢运动能力缺陷或严重丧失患者，如面向帕金森病、多发性硬化症等的神经退行性疾病以及中风、瘫痪或截肢患者，分析其步进过程的大脑状态，获取高精准的生物标志物，并开发高效神经解码算法，从而构建智能康复***辅助行走或实现神经调控，旨在改善病患生活质量。
过去针对步态的分析多使用传感器在腿部、足部或关节处记录加速度或角度信号来进行，从中提取步进周期的步幅、步速、步态变异性等运动学信息进行病理分析。该类特征固然能够直观反映步态变化，但难以揭示神经环路的生理或病理状态。而脑电信号在实时性、连续性方面具有较好特性，其中蕴含着丰富的时频和幅相信息，能够反映大脑在特定状态下的工作模式与功能活动，对于疾病表征十分有效。另一方面，以帕金森病为例，研究指出丘脑底核的神经活动与帕金森病的症状密切相关，即丘脑底核中低频活动的爆发限制了基底神经节***编码有关预期动作生理相关信息的能力，而应用局部场势动力学能够解码下丘脑活动，在运动执行之前准确预测预期运动，表明脑电信号用于运动解码具有可行性。另一方面，β频带神经活动的过度同步以及β功率增加与运动障碍有关，其中β爆发的时间特性，如振幅、持续时间等更能反映运动损伤。在过去基于阈值模型检测上述爆发活动的算法中，阈值多定义为归一化信号幅度的75％，爆发持续时间定义为超过设定阈值的时间，该类方法过度依赖阈值选择，且需要预先对信号频带加以限制，而隐马尔可夫模型HMM克服了阈值选择的任意性这一缺点，并能够同时检测多种爆发活动，十分适用于分析状态复杂的脑电信号。
此外，传统脑电分析使用直接滤波方法提取对应频带分量，并基于傅里叶变换进行频域分析。傅里叶变换的应用前提为信号平稳，且该方法仅能获得频率这一单维信息，难以准确表征高动态脑电信号。
在针对脑电信号的解码过程中，现有的方法使用隐马尔可夫模型HMM进行状态识别，适用于分析状态复杂的脑电信号，可以对多种爆发活动进行检测。当前已有应用脑电信号进行步态障碍的分析识别，使用隐马尔可夫模型HMM进行状态识别，但是该方法无法保留脑电信号的高动态变化和非线性特征，也无法自适应地解码处不同频带下的脑状态。
因此如何对脑电信号中关于步态的特征进行精细解码和检测，在保留脑电信号非线性的同时，自适应地实现不同频带下大脑状态的精细解码，是目前亟待解决的问题。
发明内容
有鉴于此，本发明提供了基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，能够将隐马尔可夫模型与掩膜经验模态分解结合，在保留脑电信号非线性的同时，自适应地实现不同频带下大脑状态的精细解码，计算出两个时域特征FO和LT。
为达到上述目的，本发明的技术方案包括如下步骤：
采集脑电信号，并进行预处理。
预处理后的脑电信号一方面进行掩膜经验模态分解、提取所需频率范围对应的本征模函数IMF分量，并计算所有IMF分量的希尔伯特包络。
预处理后的脑电信号另一方面输入隐马尔科夫模型，利用隐马尔科夫模型，通过设定状态的数量，得到每个状态的状态时间序列STC。
将所有IMF分量的希尔伯特包络，与状态时间序列STC计算互相关系数，根据最大正相关系数进行STC与脑电频带的对应，如果不存在正相关系数则归于背景频带，最终得到不同频带下分配的STC。
分别计算不同频带下分配的STC的时域特征，提取不同频带下状态占用率FO特征和状态生命时间LT特征。
进一步地，采集受试者步进过程中的脑电信号，并进行预处理，其中预处理的具体流程包括：对连续脑电数据进行50Hz陷波，然后进行8Hz的高通滤波，然后进行48Hz低通滤波，最后将原始数据降采样至100Hz，降采样后的采样率为SR。
进一步地，预处理后的脑电信号一方面进行掩膜经验模态分解、提取所需频率范围对应的本征模函数IMF分量，具体采用如下步骤：
MEMD使用正弦信号作为掩膜信号，振幅为输入信号的标准差，每一层掩膜信号的中心频率小于信号所包含的最高频率分量，每一层掩膜信号的相位均匀分布于整个周期；具体构建过程为：
预设所需频带的期望频率fd，期望频率的选择与脑电频带的划分相关，在每个脑电频带上选取一个期望频率；
对同一个输入信号以不同的期望频率完成MEMD算法，即S301-S303，每一个期望频率fd获得一个IMF分量；
S201：输入信号为预处理后的脑电信号，即s(t)，对于一个期望频率fd，输入信号进行q层分解，q计算为：q＝(SR/2)/(fd×p)，p是经验参数；SR降采样后的采样率；
每一层对应一个掩膜信号的中心频率，i的取值为1到q，q个中心频率构成一个向量fc，第i个中心频率
第i层分解的信号记为x(i)(t)，i初值取1时，第一层分解信号为x(t)(1)＝s(t)；
预先设置每层掩膜信号的相位数为l，相位在整个周期上均匀分布，则第i层第k个掩膜信号的相位k的取值为1到l，每一层共有l个掩膜信号；第i层第k个掩膜信号构造为/>的振幅A为输入信号s(t)的标准差，每一层的所有掩膜信号频率为中心频率/>
S202：构建中间信号后，对第i层第k个中间信号/>使用EMD经验模态分解得到第一级IMF分量/>循环该步骤S302，直至本层所有掩膜信号均被使用，即从/>到/>都做了EMD经验模态分解得到第一级IMF分量；
S203：将第i层的所有中间信号所得到的第一级IMF分量叠加后取均值记为d(i)，当前层的输入信号减去d(i)作为下一层的输入信号，本层循环结束，层数i加1，返回S201重复以上步骤进入下一层循环，直至最后一层中心频率得到的d(i)即为所需IMF分量，即期望频率fd对应脑电频带的时域信号分量。
进一步地，预设所需频带的期望频率fd，期望频率的选择与脑电频带的划分相关，在每个脑电频带上选取一个期望频率，具体为：
脑电频带划分为四个频带，分别为：α频带为[8Hz,12Hz]、低β频带为[13Hz,21Hz]、高β频带[22Hz,35Hz]和γ频带35Hz以上；
对应四个频带分别取对应期望频率fd分别为10Hz，18Hz，28Hz，40Hz。
进一步地，计算所有IMF分量的希尔伯特包络，具体为：
取希尔伯特包络的输入信号x(t)为分解得到的每个IMF分量，t为时刻；
首先对输入信号x(t)进行希尔伯特变换，得到然后将原输入信号作为实部，经过希尔伯特变换后的信号作为虚部，构造解析信号/>j为虚数单位；
令x(t)＝A(t)cos(ω0t+θ(t))，ω0指正弦信号的角频率，A(t)为瞬时振幅，即包络，θ(t)为瞬时相位；
带入解析信号整理得到解析信号的绝对值即为所求的希尔伯特包络。
进一步地，状态占用率FO特征为每个状态所有出现的时间在整个信号时间长度上占的比例；所述状态生命时间LT特征为一个状态未转移至其他状态前该状态的平均持续时间。
有益效果：
1、本发明提供的基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，能够将隐马尔可夫模型与掩膜经验模态分解结合，其中掩膜经验模态分解MEMD算法能够依据数据自身的时间尺度特征，将复杂的时域信号分解为有限个IMF，隐马尔可夫模型HMM用于脑电解码，则能够同时检测多种爆发活动，十分适用于分析状态复杂的脑电信号，本发明将二者结合，脑电信号非线性、非平稳特点使得上述方法适用性较低，经验模态分解EMD在处理非平稳、非线性信号上具有明显优势，能够依据数据自身时间尺度，将复杂时域信号自适应地分解为有限个本征模函数IMF，各IMF分量包含了原始信号在不同时间尺度上的局部特征信息，但该方法的两项缺点为模态混叠和端点效应。进一步改进的掩膜经验模态分解MEMD方法能够有效同时避免模态混叠和残余噪声，在脑电信号解构分析中取得了较为突出的效果。MEMD在信号分解时能够保留高动态变化的脑电信号的非线性特性，结合HMM对多种爆发活动的检测，自适应地解码出不同频带下的脑状态，实现了针对脑电信号的精细解码，计算出两个时域特征FO特征和LT特征。
2、本发明实施例提供的基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，针对MEMD算法进行改进，以设置中心频率的方式选取IMF分量。现有的MEMD算法中，通常的一个信号做一次MEMD会得到多个IMF，而本发明实施例中提出设置一个中心频率，自动选择最需要的那一个IMF，不需要在得到多个IMF后再通过肉眼观察进行选择，降低了人员肉眼观察带来的主观性差异，提高了分解的准确性。
附图说明
图1为本发明实施例提供基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法的流程图；
图2为本发明实施例中MEMD结果示例图；
图3为本发明实施例中状态时间特征计算示例图。
具体实施方式
下面结合附图并举实施例，对本发明进行详细描述。
本发明提供了基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，其技术流程如图1所示，即包括如下步骤：
步骤1)采集脑电信号，并进行预处理；所采集的脑电信号是受试者步进过程中的脑电信号。
对采集到的脑电信号进行预处理，具体操作流程包括：先对连续脑电数据进行50Hz陷波以避免工频噪声，然后进行8Hz的高通滤波；滤波后，对时域信号进行目测审视，剔除高通滤波效果不佳、动作伪迹依然明显的数据片段；接着进行48Hz低通滤波，将频率成分集中在脑电信号主要范围内；最后将原始数据降采样至100Hz，降采样后的采样率为SR。
步骤2)将预处理后的脑电信号进行MEMD分解，提取所需频率范围所对应的IMF分量。
MEMD算法能够依据数据自身的时间尺度特征，将复杂的时域信号分解为有限个IMF，分解得到的各IMF分量包含了原始信号在不同时间尺度上的局部特征信息。最终，输入信号可表示为多个IMF和一个残差项之和：
其中：n为时间点，m为IMF分量的序号，M为IMF分量的总数，I(n)为输入信号，IMFm(n)为时间点n上的第m个IMF分量，ResM(n)为时间点n对应的残差项。
本发明实施例中MEMD使用正弦信号作为掩膜信号，振幅为输入信号的标准差，每一层掩膜信号的中心频率小于信号所包含的最高频率分量，，每一层掩膜信号的相位均匀分布于整个周期。具体构建过程为：
预设所需频带的期望频率fd，期望频率的选择与脑电频带的划分相关，目前经典脑电频带通常划分为：α(8-12Hz)、低β(13-21Hz)、高β(22-35Hz)和γ(>35Hz)，这里对应不同频带取期望频率为10Hz，18Hz，28Hz，40Hz。对同一个输入信号以不同的期望频率完成MEMD算法，即S301-S303，每一个期望频率fd获得一个IMF分量。
S201：输入信号为s(t)，这里为预处理后的脑电信号。对于一个期望频率fd，输入信号需要进行q层分解，q计算为：q＝(SR/2)/(fd×p)，p是经验参数，需调试。每一层对应一个掩膜信号的中心频率，i的取值为1到q，q个中心频率构成一个向量fc，第i个中心频率现有的MEMD算法中，通常的一个信号做一次MEMD会得到多个IMF，而本发明实施例中在此处提出设置一个中心频率，自动选择最需要的那一个IMF，不需要在得到多个IMF后再通过肉眼观察进行选择，降低了人员肉眼观察带来的主观性差异，提高了分解的准确性。
第i层分解的信号记为x(i)(t)，第一层分解信号为x(t)(1)＝s(t)。
预先设置每层掩膜信号的相位数为l，相位在整个周期上均匀分布，则第i层第k个掩膜信号的相位k的取值为1到l，每一层共有l个掩膜信号。第i层第k个掩膜信号构造为/>的振幅A为输入信号s(t)的标准差，每一层的所有掩膜信号频率为中心频率/>
S202：构建中间信号后，对第i层第k个中间信号/>使用EMD经验模态分解得到第一级IMF分量/>循环该步骤S302，直至本层所有掩膜信号均被使用，即从/>到/>都做了EMD得到第一级IMF分量。
S203：将第i层的所有中间信号所得到的第一级IMF分量叠加后取均值记为d(i)，当前层的输入信号减去d(i)作为下一层的输入信号，本层循环结束，层数i加1，返回S201重复以上步骤进入下一层循环，直至最后一层(第q层)中心频率得到的d(i)即为所需IMF分量，即期望频率fd对应脑电频带的时域信号分量。
本发明实施例中，最终分别提取出α(8-12Hz)、低β(13-21Hz)、高β(22-35Hz)和γ(>35Hz)四个频带的IMF分量。
图2给出了某一被试的脑电信号经过4次MEMD分解得到的结果，左边四个子图为每个期望频率提取的IMF，即对应所需要的脑电频带的时域信号分量(截取1s)，右边四个子图表示提取的4个IMF计算的功率谱密度，中心频率落在α(8-12Hz)、低β(13-21Hz)、高β(22-35Hz)和γ(>35Hz)频带内，即分解结果与所需频带实现对应。
步骤3)将步骤1)中预处理后的脑电信号作为隐马尔可夫HMM模型的输入，推断出预定数量的状态时间序列STC和维特比路径I*。本发明实施例给出HMM的模型介绍和具体推断过程：
HMM是一种时间序列概率模型，模型假设一个时间序列可以通过一个包含有限数量状态的隐藏序列来描述，在每个时间点推断出各状态出现的概率，并认为概率最大的状态为当前时间点的激活状态。
定义观测序列为O＝{o1,o2,,...,oT}，T为时刻，对每一个时刻t，ot所有可能的取值为集合{v1,v2,,...,vM}，M为可能的观测值个数；定义状态序列为I＝{i1,i2,,...,iT}，对每一个时刻t，it所有可能的取值为集合{q1,q2,,...,qN}，N为可能的状态个数。一般地，HMM记作λ＝[π,A,B]，模型三要素分别为：
①初始状态概率π，即模型初始时刻处于各状态的概率，通常记为π＝(π1,π2,...,πN),N为状态数，πi表示初始时刻处于状态i的概率，i取值为1到N；
②状态转移概率A，即模型在各个状态之间转换的概率，通常记为矩阵A[aij],aij表示t时刻状态为qi，t+1时刻状态为qj的概率；
③输出观测概率B，即模型根据当前状态获取各个观测值的概率，通常记为矩阵B[bij]，bij表示任意时刻t，当前状态为qi，则观测值为vj的概率。
HMM模型包含三个基本问题：学习问题、评估问题和预测问题，这里基于学习问题求解模型参数并得到STC，实现步态下的脑状态解码。
学习问题，即给定观测序列O＝{o1,o2,,...,oT}，这里观测序列为步骤2中得到的预处理后的任意单通道脑电信号，求解模型参数λ＝(A,B,π)使得P(O|λ)最大，P(O|λ)表示在模型参数为λ条件下，观测序列出现的概率，求解过程如下：
已知观测序列为O＝{o1,o2,,...,oT}，其中观测序列O的隐藏状态序列为I＝{i1,i2,,...,iT}，完全数据为(O,I)＝(o1,o2,...,oT,i1,i2,,...,iT)，则P(O,I|λ)指在模型参数为λ条件下，观测序列为O，状态序列为I的概率，其对数似然函数为logP(O,I|λ)。
定义γt(i)为给定模型参数λ和观测序列O，在时刻t状态为qi的概率，即
γt(i)＝P(st＝qi|O,λ) (2)
定义ξt(i,j)为给定模型λ和观测序列O，在时刻t状态为qi，在时刻t+1状态为qj的概率，即
ξt(i,j)＝P(st＝qi,st+1＝qj|O,λ) (3)
根据EM算法的E步，构造函数(EM算法包含在HMM模型里，是使用HMM模型公认的求解学习问题的经典算法)，这里Q只是一个表示函数的符号，相当于取个名字，也可以用任意一个字母，比如/>表示这个函数是关于λ和/>这两个变量的函数，在EM算法里默认用Q，其中/>是模型参数当前的估计值，λ是需要最大化的模型参数，则
根据EM算法的M步，极大化函数迭代收敛得到最终的模型参数λ＝(A,B,π)。
初始化模型参数，对n＝0选取得到λ(0)＝(A(0),B(0),π(0))；
对n＝1,2,...递推得到参数计算公式如下：
公式(6)中I(ot＝vj)表示当时刻t观测值为vj时，I(ot＝vj)＝1，否则I(ot＝vj)＝0；
终止递推，得到模型参数λ＝(A,B,π)。
定义t时刻各状态的前向概率为αt(i)＝P(o1,o2,...,ot,it＝qi|λ)，αt(i)表示在模型参数λ条件下，时刻t时观测序列为o1,o2,...,ot，状态为qi的概率；定义t时刻各状态的后向概率为βt(i)＝P(ot+1,ot+2,...,oT|it＝qi,λ)，βt(i)表示在模型参数为λ时刻t时状态为qi条件下，从时刻t+1到最后时刻T的观测序列为ot+1,ot+2,...,oT的概率，i的取值为1到N。
求解得到λ后，可根据前向概率和后向概率计算任意时刻t各状态出现的概率为：
P(it＝qi,O|λ)＝αt(i)βt(i) (8)
依据时间点递推得到所有时刻各状态出现的概率，每一个状态(共N个状态)都得到一个状态时间序列STC。
步骤4)计算步骤2)得到的所有IMF的希尔伯特包络，并与步骤3)中得到的STC计算互相关，根据最大正相关系数进行STC与脑电频带的对应，如果不存在正相关系数则归于背景频带，最终得到五个频带。
S401：希尔伯特包络计算过程如下：
首先对输入信号x(t)进行希尔伯特变换，得到然后将原输入信号作为实部，经过希尔伯特变换后的信号作为虚部，构造解析信号/>j为虚数单位，令x(t)＝A(t)cos(ω0t+θ(t))(ω0指正弦信号的角频率，是三角函数在数学上的固定概念)，带入解析信号整理得到/>A(t)为瞬时振幅，即包络，θ(t)为瞬时相位；解析信号的绝对值即为所求的希尔伯特包络。这里的输入信号x(t)为分解得到的每个IMF分量。
S402：相关系数r计算公式如下：
其中，X和Y分别表示两个时间序列，这里X为S501中得到的各IMF的希尔伯特包络，Y为步骤4中HMM模型得到的STC。
S403)每一个STC与4个IMF分别计算相关系数，得到4个相关系数值，如果4个相关系数均为负值，则该STC属于背景活动频带；如果4个相关系数中存在正值，则最大的正相关系数所对应的IMF频带就是该STC的最终所属频带。比如：若STC1与IMF1-α，IMF2-低β，IMF3-高β，IMF4-γ计算得到的相关系数分别为0.1，0.2，0.15，0.3，则STC1属于γ频带；若STC2与IMF1-α，IMF2-低β，IMF3-高β，IMF4-γ计算得到的相关系数分别为-0.1，-0.2，-0.15，-0.3，则STC2属于背景频带。
步骤5)分别计算不同频带下分配的STC的时域特征，提取不同频带下状态占用率FO特征和状态生命时间LT特征。则通过对不同频带下的FO特征和LT特征进行分析可以判断输入的脑电信号是否存在步态障碍：
①状态占用率(Fractional Occupancy,FO)，定义为每个状态所有出现的时间在整个信号时间长度上占的比例；对于状态k的占用率计算公式为：
其中，T是所有状态序列的时间点数，t为指定状态所对应的时间点数。
②状态生命时间(Life Times,LT)，定义为某状态未转移至其他状态前该状态的平均持续时间；对于状态k的生命时间计算公式为：
上式中occ(k)表示状态k的出现次数，计算公式为：
其中，threshold为门限值，可以通过设置门限来摒弃过短的时间片段，不计入状态k的出现次数。
图3给出某一状态依据STC计算时间特征LT的方法示意，对检测到的爆发活动设置0.2s的门限。
本发明已在帕金森步态脑电数据集上进行验证。基于所提出算法对患者在三种声音提示下踏步时记录的脑电数据进行分析，实现了脑电状态的识别以及不同频带下状态时域特征的比较，有效识别出步态障碍的病理特征并对比了不同声音条件下的特征差异，即声音提示显著降低β频带下爆发活动的FO及LT，与之相反，α频带下爆发活动的FO及LT增加。
因此可以看出，基于掩膜经验模态分解与隐马尔可夫模型的步态下脑电解码方法可以作为临床步态分析的有力支撑，有助于对步态障碍事件的检测，在康复训练***的开发中具有广泛的应用前景。
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (6)
1.基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，其特征在于，包括如下步骤：
采集脑电信号，并进行预处理；
预处理后的脑电信号一方面进行掩膜经验模态分解、提取所需频率范围对应的本征模函数IMF分量，并计算所有IMF分量的希尔伯特包络；
预处理后的脑电信号另一方面输入隐马尔科夫模型，利用隐马尔科夫模型，通过设定状态的数量，得到每个状态的状态时间序列STC；
将所有IMF分量的希尔伯特包络，与状态时间序列STC计算互相关系数，根据最大正相关系数进行STC与脑电频带的对应，如果不存在正相关系数则归于背景频带，最终得到不同频带下分配的STC；
分别计算不同频带下分配的STC的时域特征，提取不同频带下状态占用率FO特征和状态生命时间LT特征。
2.如权利要求1所述的基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，其特征在于，所述采集受试者步进过程中的脑电信号，并进行预处理，其中预处理的具体流程包括：对连续脑电数据进行50Hz陷波，然后进行8Hz的高通滤波，然后进行48Hz低通滤波，最后将原始数据降采样至100Hz，降采样后的采样率为SR。
3.如权利要求1所述的基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，其特征在于，所述预处理后的脑电信号一方面进行掩膜经验模态分解、提取所需频率范围对应的本征模函数IMF分量，具体采用如下步骤：
MEMD使用正弦信号作为掩膜信号，振幅为输入信号的标准差，每一层掩膜信号的中心频率小于信号所包含的最高频率分量，每一层掩膜信号的相位均匀分布于整个周期；具体构建过程为：
预设所需频带的期望频率fd，期望频率的选择与脑电频带的划分相关，在每个脑电频带上选取一个期望频率；
对同一个输入信号以不同的期望频率完成MEMD算法，即S301-S303，每一个期望频率fd获得一个IMF分量；
S201：输入信号为预处理后的脑电信号，即s(t)，对于一个期望频率fd，输入信号进行q层分解，q计算为：q＝(SR/2)/(fd×p)，p是经验参数；SR降采样后的采样率；
每一层对应一个掩膜信号的中心频率，i的取值为1到q，q个中心频率构成一个向量fc，第i个中心频率
第i层分解的信号记为x(i)(t)，i初值取1时，第一层分解信号为x(t)(1)＝s(t)；
预先设置每层掩膜信号的相位数为l，相位在整个周期上均匀分布，则第i层第k个掩膜信号的相位k的取值为1到l，每一层共有l个掩膜信号；第i层第k个掩膜信号构造为/> 的振幅A为输入信号s(t)的标准差，每一层的所有掩膜信号频率为中心频率/>
S202：构建中间信号后，对第i层第k个中间信号/>使用EMD经验模态分解得到第一级IMF分量/>循环该步骤S302，直至本层所有掩膜信号均被使用，即从/>到/>都做了EMD经验模态分解得到第一级IMF分量；
S203：将第i层的所有中间信号所得到的第一级IMF分量叠加后取均值记为d(i)，当前层的输入信号减去d(i)作为下一层的输入信号，本层循环结束，层数i加1，返回S201重复以上步骤进入下一层循环，直至最后一层中心频率得到的d(i)即为所需IMF分量，即期望频率fd对应脑电频带的时域信号分量。
4.如权利要求1～3任一所述的基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，其特征在于，所述预设所需频带的期望频率fd，期望频率的选择与脑电频带的划分相关，在每个脑电频带上选取一个期望频率，具体为：
脑电频带划分为四个频带，分别为：α频带为[8Hz,12Hz]、低β频带为[13Hz,21Hz]、高β频带[22Hz,35Hz]和γ频带35Hz以上；
对应四个频带分别取对应期望频率fd分别为10Hz，18Hz，28Hz，40Hz。
5.如权利要求1～3任一所述的基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，其特征在于，所述计算所有IMF分量的希尔伯特包络，具体为：
取希尔伯特包络的输入信号x(t)为分解得到的每个IMF分量，t为时刻；
首先对输入信号x(t)进行希尔伯特变换，得到然后将原输入信号作为实部，经过希尔伯特变换后的信号作为虚部，构造解析信号/>j为虚数单位；
令x(t)＝A(t)cos(ω0t+θ(t))，ω0指正弦信号的角频率，A(t)为瞬时振幅，即包络，θ(t)为瞬时相位；
带入解析信号整理得到解析信号的绝对值即为所求的希尔伯特包络。
6.如权利要求1～3任一所述的基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法，其特征在于，所述状态占用率FO特征为每个状态所有出现的时间在整个信号时间长度上占的比例；所述状态生命时间LT特征为一个状态未转移至其他状态前该状态的平均持续时间。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310776366.7A CN116662742B (zh) | 2023-06-28 | 2023-06-28 | 基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310776366.7A CN116662742B (zh) | 2023-06-28 | 2023-06-28 | 基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116662742A true CN116662742A (zh) | 2023-08-29 |
| CN116662742B CN116662742B (zh) | 2024-07-12 |
Family
ID=87726173
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310776366.7A Active CN116662742B (zh) | 2023-06-28 | 2023-06-28 | 基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116662742B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117814756A (zh) * | 2024-01-10 | 2024-04-05 | 湖北大学 | 一种基于深度学习检测模型的特征处理方法及*** |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150199010A1 (en) * | 2012-09-14 | 2015-07-16 | Interaxon Inc. | Systems and methods for collecting, analyzing, and sharing bio-signal and non-bio-signal data |
| CN105068663A (zh) * | 2015-09-18 | 2015-11-18 | 中国石油大学(华东) | 一种基于脑电信号的对象选择方法及装置 |
| CN107256393A (zh) * | 2017-06-05 | 2017-10-17 | 四川大学 | 基于深度学习的一维生理信号的特征提取与状态识别 |
| CN108805089A (zh) * | 2018-06-14 | 2018-11-13 | 南京云思创智信息科技有限公司 | 基于多模态的情绪识别方法 |
| CN108899050A (zh) * | 2018-06-14 | 2018-11-27 | 南京云思创智信息科技有限公司 | 基于多模态情绪识别***的语音信号分析子*** |
| US20190113973A1 (en) * | 2012-09-14 | 2019-04-18 | Interaxon Inc | Systems and methods for collecting, analyzing, and sharing bio-signal and non-bio-signal data |
| CN110163128A (zh) * | 2019-05-08 | 2019-08-23 | 南京邮电大学 | 改进的emd算法结合小波包变换及csp算法的脑电信号分类方法 |
| CN110604578A (zh) * | 2019-09-04 | 2019-12-24 | 平顶山学院 | 基于semg的人手手内动作识别方法 |
| CN111712193A (zh) * | 2018-02-14 | 2020-09-25 | 塞勒尼恩有限公司 | 用于脑电图测量的设备和方法 |
| CN113180596A (zh) * | 2021-04-07 | 2021-07-30 | 中山大学 | 非接触式的睡眠分析方法、装置及存储介质 |
| CN114578963A (zh) * | 2022-02-23 | 2022-06-03 | 华东理工大学 | 一种基于特征可视化和多模态融合的脑电身份识别方法 |
-
2023
- 2023-06-28 CN CN202310776366.7A patent/CN116662742B/zh active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150199010A1 (en) * | 2012-09-14 | 2015-07-16 | Interaxon Inc. | Systems and methods for collecting, analyzing, and sharing bio-signal and non-bio-signal data |
| US20190113973A1 (en) * | 2012-09-14 | 2019-04-18 | Interaxon Inc | Systems and methods for collecting, analyzing, and sharing bio-signal and non-bio-signal data |
| CN105068663A (zh) * | 2015-09-18 | 2015-11-18 | 中国石油大学(华东) | 一种基于脑电信号的对象选择方法及装置 |
| CN107256393A (zh) * | 2017-06-05 | 2017-10-17 | 四川大学 | 基于深度学习的一维生理信号的特征提取与状态识别 |
| CN111712193A (zh) * | 2018-02-14 | 2020-09-25 | 塞勒尼恩有限公司 | 用于脑电图测量的设备和方法 |
| CN108805089A (zh) * | 2018-06-14 | 2018-11-13 | 南京云思创智信息科技有限公司 | 基于多模态的情绪识别方法 |
| CN108899050A (zh) * | 2018-06-14 | 2018-11-27 | 南京云思创智信息科技有限公司 | 基于多模态情绪识别***的语音信号分析子*** |
| CN110163128A (zh) * | 2019-05-08 | 2019-08-23 | 南京邮电大学 | 改进的emd算法结合小波包变换及csp算法的脑电信号分类方法 |
| CN110604578A (zh) * | 2019-09-04 | 2019-12-24 | 平顶山学院 | 基于semg的人手手内动作识别方法 |
| CN113180596A (zh) * | 2021-04-07 | 2021-07-30 | 中山大学 | 非接触式的睡眠分析方法、装置及存储介质 |
| CN114578963A (zh) * | 2022-02-23 | 2022-06-03 | 华东理工大学 | 一种基于特征可视化和多模态融合的脑电身份识别方法 |
Non-Patent Citations (4)
| Title |
|---|
| 于沐涵;陈峰;: "基于HCSP和模糊熵的脑电信号分类", 计算机工程与设计, no. 02, 16 February 2018 (2018-02-16) * |
| 周百新;梁忠诚;赵阳;王蔚;赵新红;: "EMD及其在脑电分析中的应用研究", 南京师范大学学报(工程技术版), no. 01, 20 March 2007 (2007-03-20) * |
| 张向刚;唐海;付常君;石宇亮;: "一种基于隐马尔科夫模型的步态识别算法", 计算机科学, no. 07, 15 July 2016 (2016-07-15) * |
| 王成龙;韦巍;李天永;: "基于IMF能量矩的脑电情绪特征提取研究", 现代电子技术, no. 20, 15 October 2018 (2018-10-15) * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117814756A (zh) * | 2024-01-10 | 2024-04-05 | 湖北大学 | 一种基于深度学习检测模型的特征处理方法及*** |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116662742B (zh) | 2024-07-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108714026B (zh) | 基于深度卷积神经网络和在线决策融合的细粒度心电信号分类方法 | |
| Stamkopoulos et al. | ECG analysis using nonlinear PCA neural networks for ischemia detection | |
| Yang et al. | Spatiotemporal differentiation of myocardial infarctions | |
| CN107837082A (zh) | 基于人工智能自学习的心电图自动分析方法和装置 | |
| CN108478209A (zh) | 心电信息动态监护方法和动态监护*** | |
| Miao et al. | A spatial-frequency-temporal optimized feature sparse representation-based classification method for motor imagery EEG pattern recognition | |
| CN108681396A (zh) | 基于脑-肌电双模态神经信号的人机交互***及其方法 | |
| CN116662742B (zh) | 基于隐马尔可夫模型与掩膜经验模态分解的脑电解码方法 | |
| Yang et al. | A novel deep learning scheme for motor imagery EEG decoding based on spatial representation fusion | |
| Liang et al. | Obstructive sleep apnea detection using combination of CNN and LSTM techniques | |
| CN110141215A (zh) | 降噪自编码器的训练方法、心电信号的降噪方法及相关装置、设备 | |
| CN105534517A (zh) | 一种去除三导联心电信号中车辆运动噪声的方法 | |
| Thomas et al. | Heart sound segmentation using fractal decomposition | |
| Azar | Neuro-fuzzy system for cardiac signals classification | |
| Lee et al. | Classification of both seizure and non-seizure based on EEG signals using hidden Markov model | |
| Narayana et al. | ECG based biometric authentication system using deep learning methods | |
| Zhang et al. | Movement recognition via channel-activation-wise sEMG attention | |
| Busia et al. | A Tiny Transformer for Low-Power Arrhythmia Classification on Microcontrollers | |
| CN113288050A (zh) | 基于图卷积网络的多维增强癫痫发作预测*** | |
| Murthy et al. | Design and implementation of hybrid techniques and DA-based reconfigurable FIR filter design for noise removal in EEG signals on FPGA | |
| Bhattacharya | Detection of weak chaos in infant respiration | |
| Siekierski et al. | Heart beats classification method using a multi-signal ECG spectrogram and convolutional neural network with residual blocks | |
| Yao et al. | Combining Rhythm Information between Heartbeats and BiLSTM‐Treg Algorithm for Intelligent Beat Classification of Arrhythmia | |
| Jyotishi et al. | An Attention Based Hierarchical LSTM Architecture for ECG Biometric System | |
| Jayanthi et al. | Hybrid CNN-RNN in Motor Imagery Identification of Brain-Computer Interface |
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 |