CN111568409B - 基于双谱分析和图傅里叶变换的心电信号特征提取方法 - Google Patents

Abstract

本发明公开了一种基于双谱分析和图傅里叶变换的心电信号特征提取方法，属于心电信号特征提取方法技术领域。本方法首先利用高阶谱算法将时域信号转换到高阶谱域，再利用图傅里叶变换，将双谱矩阵转换至特征值域，然后从特征谱矩阵中直接提取图谱特征在得到特征向量后可以根据AAMI标准，对图谱特征进行分类。本方法克服了现有双谱矩阵特征提取方法计算量大的问题，并且准确度高，能够对各种心电信号进行有效分析。

Description

基于双谱分析和图傅里叶变换的心电信号特征提取方法

技术领域

本发明涉及一种基于双谱分析和图傅里叶变换的心电信号特征提取方法，属于心电信号特征提取方法技术领域。

背景技术

心血管疾病已逐渐成为危害人类生命的最常见疾病之一，因此如何诊断和预防这类疾病已成为当今医学界面临的重要问题。心电信号不仅可以用来分析和鉴别心律失常、心肌梗死等多种疾病，还可以反映心肌细胞损伤程度、发育过程、心房和心室功能结构等，已成为心血管疾病诊断中一种简明有效的工具。

信号的高阶谱是一种非平稳信号分析工具，双谱作为高阶谱中阶次最低的，不仅包含了高阶谱的所有特性，并且计算最简单，因此得到广泛应用。在实际操作中，双谱矩阵较为复杂，目前大多数特征提取方法是在得到双谱矩阵的基础上，需要进一步地降维处理，例如围线积分、双谱切片、主成分分析、独立成分分析、核主成分分析等。

作为代数图论中一个非常重要的研究方向，图谱理论早在上个世纪五十年代就开始发展。图谱理论简单的说就是在图和矩阵之间建立起对应的关系，这是图谱理论的基本思想。这样把图的问题转移到存储矩阵之后，就可以通过研究矩阵谱的相关属性来研究图的问题。存储矩阵主要包括邻接矩阵和拉普拉斯矩阵。邻接矩阵的研究已经具有较长的历史，是一门相对成熟的学科。相较于邻接矩阵，近年来出现了更多关于图的拉普拉斯矩阵谱的研究。并且由图谱理论衍生的图信号处理也是近几年快速发展的一个应用。图傅里叶变换是图信号关于图拉普拉斯矩阵特征函数的展开，也是图信号处理的基础。

通过图傅里叶变换将建立在双谱基础上的双谱矩阵转化到特征值域，从特征值的统计量中提取特征，利用图谱特征构造特征向量的方法计算简单，可区分度高。

发明内容

为解决现有技术的不足，本发明提供出了一种基于双谱分析和图傅里叶变换的心电信号特征提取方法，不仅克服了现有双谱矩阵特征提取方法计算量大的问题，并且准确度高，能够对各种心电信号进行有效分析，实现对心律不齐等心血管疾病问题的有效诊断。

本发明为解决其技术问题采用如下技术方案：

一种基于双谱分析和图傅里叶变换的心电信号特征提取方法，包含如下步骤：

步骤A.心电信号预处理，使用滤波方法去除信号中的噪声，分离得到单个心拍数据；

步骤B.对心电信号进行高阶谱分析，计算得到信号的双谱矩阵S；

步骤C.利用二维图傅里叶变换将双谱矩阵映射到特征值域；

步骤D.对获取的特征值谱

进行图谱特征提取，并组成特征向量，包括图谱平坦度、图谱亮度和图谱滚降度。

所述步骤A的具体过程如下：

步骤(A-1).将原始心电信号减去平均值消除直流分量；

步骤(A-2).通过中值滤波器消除基线漂移；

步骤(A-3).通过低通滤波器消除工频干扰和肌电噪声；

步骤(A-4).通过高通滤波器消除低频噪声；

步骤(A-5).寻找信号QRS复合波中的峰值R波位置，以峰值点为原点，该点和之前以及之后各取K个数据，形成长度为2K点的单个心拍数据h(t)。

所述步骤B的具体过程如下：

步骤(B-1).计算出单个心拍数据h(t)的三阶积累量

R_3s(τ₁,τ₂)＝E[h(t)h(t+τ₁)h(t+τ₂)]，τ₁，τ₂∈(-∞，∞)

其中：h(t)为单个心拍数据，τ₁和τ₂为时间延迟；h(t+τ₁)和h(t+τ₂)为经过时间延迟τ₁，τ₂后的时移信号，最后得到三阶累计量R_3s(τ₁,τ₂)；

步骤(B-2).计算出单个心拍数据h(t)的双谱

式中，S为h(t)的双谱矩阵，ω₁，ω₂是两个独立的频率。

所述步骤C的具体过程如下：

步骤(C-1).计算对角矩阵D，点d(m,n)的计算公式为：

对角矩阵D由d(m,n)组合而成，w(i,j)为节点i和节点j之间的连接关系；

步骤(C-2).计算拉普拉斯矩阵L

L＝D-W

式中：W为加权邻接矩阵，由元素w(m,n)组合而成，w(m,n)显示了节点m和节点n之间的连接关系，大小与节点之间距离有关；L是实对称矩阵；

步骤(C-3).计算对应于L的拉普拉斯特征值λ_l和拉普拉斯特征矩阵X_l，λ_l和X_l满足：

LX_l＝λ_lX_l,l＝0,1,...,N-1

其中，x_l(m,n)为矩阵X_l中第m行、第n列的元素；

步骤(C-4).计算双谱矩阵S的二维图傅里叶变换

得到N阶方阵

为矩阵中第m行，第n列的元素，x_l(i,j)为矩阵X_l中第i行，第j列的元素，s(i,j)为双谱矩阵S中第i行，第j列的元素。

所述步骤D的具体过程如下：

步骤(D-1).计算图谱平坦度GSF

式中，N是特征谱矩阵的点数，

为矩阵G中第p行第q列的元素；

步骤(D-2).计算图谱亮度GSB

F为给定的边界特征值；

步骤(D-3).计算图谱滚降度GSR

这里β是系数。

最后，由特征值构成特征向量，根据AAMI标准对图谱特征进行分类。

本发明的有益效果如下：

1、本发明在预处理的基础上，采用双谱作为特征提取的第一步，利用双谱的优良性质，理论上消除了高斯噪声的影响。

2、使用二维图傅里叶变换作为第二步，将双谱矩阵转换为特征值谱，各类特征差距变大，可区分度提高。在特征值谱上直接提取特征，避免了特征的二次提取引入的误差和计算量。

附图说明

图1是心电信号特征提取步骤框图。

图2(a)是第100号记录的单个心拍信号预处理前的时域波形图，图2(b)是第100号记录的单个心拍信号预处理后的时域波形图。

图3是第100号记录的心电信号的单个心拍双谱图。

图4是经过二维图傅里叶变换之后，第100号记录的单个心拍的特征值图谱。

图5(a)是N类心电信号经过二维图傅里叶变换，转换到特征值域后提取的特征的箱型图，图5(b)是S类心电信号经过二维图傅里叶变换，转换到特征值域后提取的特征的箱型图，图5(c)是V类心电信号经过二维图傅里叶变换，转换到特征值域后提取的特征的箱型图。

具体实施方式

下面结合附图对本发明创造做进一步的说明。以MIT-BIH数据库中的数据为实施例，实施例包括以下步骤：

步骤A.心电信号预处理，使用多种滤波器方法去除信号中的噪声，分离得到单个心拍数据。

步骤(A-1).将原始心电信号减去平均值来消除直流分量。

步骤(A-2).通过中值滤波器来消除基线漂移。

步骤(A-3).通过低通滤波器来消除工频干扰和肌电噪声。

步骤(A-4).通过高通滤波器来消除低频噪声。

步骤(A-5).寻找信号QRS(QRS波代表两个心室兴奋传播过程的电位变化)复合波中的峰值R波位置。以峰值点为原点，该点和之前以及之后各取K＝100个数据，形成长度为200点的单个心拍数据h(t)。

如图1所示为信号特征分类过程的框图。共分为四个步骤，首先是预处理得到单个心拍的信号，之后经过双谱计算得到对应的双谱矩阵，再采用二维图傅里叶变换，最后提取图谱特征。

如图2所示为数据库中第100号记录的心电信号预处理前后对比，采样频率为360Hz。

步骤B.对心电信号进行高阶谱分析，计算得到信号的双谱矩阵S,矩阵中第m行，第n列元素为s(m,n)。

步骤(B-1).计算出单个心拍数据h(t)的三阶积累量。

R_3s(τ₁,τ₂)＝E[h(t)h(t+τ₁)h(t+τ₂)]，τ₁，τ₂∈(-∞，∞)

其中：h(t)为单个心拍数据，τ₁和τ₂为时间延迟；h(t+τ₁)和h(t+τ₂)为经过时间延迟τ₁，τ₂后的时移信号，最后得到三阶累计量R_3s(τ₁,τ₂)；

步骤(B-2).计算出单个心拍数据h(t)的双谱。

式中，S为h(t)的双谱矩阵，ω₁，ω₂是两个独立的频率。

得到单个心拍数据信号的双谱分析结果，如图3所示为数据库中第100号记录的心电信号的单个心拍双谱图。

步骤C.利用二维图傅里叶变换将双谱矩阵映射到特征值域。在二维图傅里叶变换中，无向图G(V,E,W),其中的V代表顶点的集合，E代表边集合。W为加权邻接矩阵，由元素w(m,n)组合而成，w(m,n)显示了节点m和节点n之间的连接关系，大小与节点之间距离有关。这里选择：

式中，m,n＝1,2,...,N。

步骤(C-1).计算对角矩阵D。点d(m,n)的计算公式为：

对角矩阵D由d(m,n)组合而成，w(i,j)为节点i和节点j之间的连接关系。

步骤(C-2).计算拉普拉斯矩阵L。

L＝D-W

式中L是实对称矩阵。

步骤(C-3)计算对应于L的拉普拉特征值λ_l和拉普拉斯特征矩阵X_l，λ_l和X_l满足：

LX_l＝λ_lX_l,l＝0,1,...,N-1

x_l(m,n)为矩阵X_l中第m行，第n列的元素。

步骤(C-4).计算双谱矩阵S的二维图傅里叶变换。

得到N阶特征谱矩阵

x_l(i,j)为矩阵X_l中第i行，第j列的元素，s(i,j)为双谱矩阵S中第i行，第j列的元素。

得到单个心拍数据双谱图的图傅里叶变换结果，如图4所示为数据库中第100号记录的单个心拍数据双谱图经过图傅里叶变换之后的特征值谱图。

步骤D.对获取的特征值谱

进行图谱特征提取，并组成特征向量。

包括：图谱平坦度(GSF)、图谱亮度(GSB)和图谱滚降度(GSR)。步骤(D-1).计算图谱平坦度

式中，N是特征谱矩阵的点数，这里N＝100，

为矩阵G中第p行第q列的元素。

步骤(D-2).计算图谱亮度

F为给定的边界特征值，这里取为25。

步骤(D-3).计算图谱滚降度

这里β是系数，选择为0.3。

得到单个心拍数据信号在特征值域的图谱特征，图5所示为MIT-BIH数据库中根据美国医疗器械促进协会(Association for the Advancement of MedicalInstrumentation,AAMI)标准划分得到四类心电信号的图谱平坦度，图谱亮度和图谱滚降度的箱型分布图，其中：N类(正常或者束支传导阻滞节拍)共67994个；S类(室上性异常节拍)共2577个；V类(心室异常节拍)共4249个；F类(融合节拍)共784个。从图中可以看出，四类信号的图谱平坦度和图谱滚降度特征有较明显的区分。最后，由特征值构成特征向量，根据AAMI标准对图谱特征进行分类。

Claims (4)

1.一种基于双谱分析和图傅里叶变换的心电信号特征提取方法，其特征在于，包含如下步骤：

步骤A.心电信号预处理，使用滤波方法去除信号中的噪声，分离得到单个心拍数据；

步骤B.对心电信号进行高阶谱分析，计算得到信号的双谱矩阵S；步骤C.利用二维图傅里叶变换将双谱矩阵映射到特征值域，具体过程如下：

步骤(C-1).计算对角矩阵D，点d(m,n)的计算公式为：

对角矩阵D由d(m,n)组合而成，w(i,j)为节点i和节点j之间的连接关系；

步骤(C-2).计算拉普拉斯矩阵L

L＝D-W

式中：W为加权邻接矩阵，由元素w(m,n)组合而成，w(m,n)显示了节点m和节点n之间的连接关系，大小与节点之间距离有关；L是实对称矩阵；

步骤(C-3).计算对应于L的拉普拉斯特征值λ_l和拉普拉斯特征矩阵X_l，λ_l和X_l满足：

LX_l＝λ_lX_l,l＝0,1,...,N-1

其中，x_l(m,n)为矩阵X_l中第m行、第n列的元素；

步骤(C-4).计算双谱矩阵S的二维图傅里叶变换

得到N阶方阵

为矩阵中第m行，第n列的元素，x_l(i,j)为矩阵X_l中第i行，第j列的元素，s(i,j)为双谱矩阵S中第i行，第j列的元素；

步骤D.对获取的特征值谱

进行图谱特征提取，并组成特征向量，包括图谱平坦度、图谱亮度和图谱滚降度。

2.根据权利要求1所述的基于双谱分析和图傅里叶变换的心电信号特征提取方法，其特征在于，所述步骤A的具体过程如下：

步骤(A-1).将原始心电信号减去平均值消除直流分量；

步骤(A-2).通过中值滤波器消除基线漂移；

步骤(A-3).通过低通滤波器消除工频干扰和肌电噪声；

步骤(A-4).通过高通滤波器消除低频噪声；

步骤(A-5).寻找信号QRS复合波中的峰值R波位置，以峰值点为原点，该点和之前以及之后各取K个数据，形成长度为2K点的单个心拍数据h(t)。

3.根据权利要求1所述的基于双谱分析和图傅里叶变换的心电信号特征提取方法，其特征在于，所述步骤B的具体过程如下：

步骤(B-1).计算出单个心拍数据h(t)的三阶积累量

R_3s(τ₁,τ₂)＝E[h(t)h(t+τ₁)h(t+τ₂)]，τ₁，τ₂∈(-∞，∞)

其中：h(t)为单个心拍数据，τ₁和τ₂为时间延迟；h(t+τ₁)和h(t+τ₂)为经过时间延迟τ₁，τ₂后的时移信号，最后得到三阶累计量R_3s(τ₁,τ₂)；

步骤(B-2).计算出单个心拍数据h(t)的双谱

式中，S为h(t)的双谱矩阵，ω₁，ω₂是两个独立的频率。

4.根据权利要求1所述的基于双谱分析和图傅里叶变换的心电信号特征提取方法，其特征在于，所述步骤D的具体过程如下：

步骤(D-1).计算图谱平坦度GSF

式中，N是特征谱矩阵的点数，

为矩阵G中第p行第q列的元素；

步骤(D-2).计算图谱亮度GSB

F为给定的边界特征值；

步骤(D-3).计算图谱滚降度GSR

这里β是系数；

最后，由特征值构成特征向量，根据AAMI标准对图谱特征进行分类。

Images