基于分数阶傅里叶变换和张量分解的心电T波特征提取方法与流程

基于分数阶傅里叶变换和张量分解的心电t波特征提取方法
技术领域
1.本发明提出一种心电信号特征提取方法，适用于结合适合的分类器，建立分类模型，检测心电异常尤其是t波异常，属于心电信号处理领域。


背景技术：

2.体表心电图(electrocardiogram，ecg)是最常见的心脏状态无创检测方法，心电异常与恶性心律失常的发生有密切联系。基于心电图t波交替(t
‑
wave alternant,twa)现象因其基本发生机制与恶性心律失常在临床上具有相对应关系，目前被认为是一种对恶性心律失常具有重要预测作用的无创电生理检测指标。
3.t波交替(t
‑
wave alternant.twa)是指在规整心率时，体现在心电信号中t波或者st段的幅度或极性随心拍交替变化呈ababab
…
的一种心电变异现象。由于许多t波交替现象为微伏级，所以肉眼可以分辨的t波交替现象很少。所以寻找一种高精度的自动检测twa的方法就显得尤为重要。目前存在的t波交替检测方法包括频谱分析法，非线性法，符号变换法等。通常心电信号都存在噪声，在twa检测过程中要去除噪声。但是由于t波十分微弱，少量的噪声可能已经完全掩盖了t波，所以在去除噪声的过程中可能将部分t波去除。
4.张量分析是多因素分析方法之一，通过张量空间对高阶数据的分析处理，在数据挖掘、图像分析、机器学习、化学分析等领域都有普遍的运用。张量分析能够充分体现信号的结构特征及内在相关性和协同性，从而可作为特征提取和识别分类的依据和手段。因此本文提出一种基于分数阶傅里叶和张量分解的t波特征提取方法，旨在从单导联心电信号中检测t波异常，尤其是存在的微伏级的twa现象。


技术实现要素：
：
5.鉴于心电信号的时变特征与非线性特征。本发明结合frft(fractional fourier transform，分数阶傅里叶变换)时频分析的能力与张量分解的主成分分析能力，及信息熵对信号混乱程度的度量能力，提供了一种基于分数阶傅里叶变换和张量分解的心电信号特征提取方法，包括如下步骤。
6.步骤s1：连续采集包含p个心动周期的心电信号并进行预处理，获得信号d；该段心电所包含的心动周期个数p应大于等于8；
7.所述预处理包括去除心电信号中的工频干扰、肌电干扰和基线漂移。预处理还包括将信号升采样到统一的采样率。
8.步骤s2：对步骤s1预处理后的信号d进行心电的qrs波群定位，获得每个qrs波群的起点和止点及该qrs波群的r波峰值的位置信息；
9.步骤s3：构造t波矩阵
10.s3
‑
1：在每两个相邻的qrs波群之间，即在其中第一个qrs波群的止点和第二个qrs波群的起点之间，定位t波的起点和终点，并计算该段t波起点和终点之间的采样点的个数x；然后在该t波的起点和终点之间搜索极大值，以该极大值作为此t波峰值；
11.令所检测到的t波的个数为n，n≤p；令所检测到的各t波相应的x的最小值，为最小采样点数x
min
；
12.对检测到的每一个t波波峰，以t波波峰所在位置为中心采样点，分别向该t波峰值左侧和右侧各再取n个采样点，加上该t波波峰共2n+1个采样点，依时间顺序排列，构成该t波的2n+1点长的t波向量；
13.所述点数2n+1应小于等于t波起点和终点之间的最小采样点数x
min
；
14.s3
‑
2：对每个t波的2n+1点长的t波向量，分别进行m个不同阶次的分数阶傅里叶变换，第j个t波向量经过i阶的分数阶傅里叶变换后得到2n+1点长的信号d
ji
；其中j≤n，所用到的分数阶傅里叶变换的阶次0≤i≤1，并且对所有t波，所做的m个阶次分数傅里叶变换的组合是相同的；
15.s3
‑
3：计算各信号d
ji
的分数阶傅里叶幅度谱信号fv
ji
；
16.进一步地，作为优选：计算各信号d
ji
的分数阶傅里叶幅度谱信号fv
ji
之后，对每个t波，对每个分数阶傅里叶幅度谱信号fv
ji
以该阶次i下最大值最小值为参照，进行归一化处理，使每个fv
ji
的各点取值在0到1区间；得到归一化后的2n+1点长的分数阶傅里叶幅度谱信号。
17.s3
‑
4：对每个t波，将m个不同阶次的分数阶傅里叶幅度谱信号作为行向量，按照所用的分数傅里叶变换的阶次从小到大依次从上到下排列，构成m行
×
(2n+1)列的t波矩阵；步骤s3
‑
1共检测到n个t波因此有n个这样的矩阵；作为优选，m≤20；
18.步骤s4：构造t波张量空间
19.将n个t波矩阵按t波出现的顺序依次排列，得到m
×
(2n+1)
×
n的三阶张量；
20.步骤s5：张量分解
21.对上述三阶张量采用tucker分解的方法进行张量分解，获得核心张量g和三个投影方向上的投影矩阵，对应分数阶傅里叶变换阶次方向的投影矩阵u1、t波向量方向对应阶的投影矩阵u2，和心拍方向对应的投影矩阵u3；
22.步骤s6：
23.对心拍方向对应的投影矩阵u3，求熵值，以此作为该段心电的t波特征。
24.作为优选，步骤s6所述熵值是一个熵向量，计算方法为：
25.对投影矩阵u3每一行计算熵值，将获得的熵值依所在行号依次排列，获得该投影矩阵的熵向量。
26.或者，步骤s6所述熵值计算方法为：
27.对投影矩阵u3每一列求均值，获得一维的行向量，对该行向量求熵值，获得该投影矩阵的熵值。
28.投影矩阵u3每一行的熵值计算方法如下：
29.将投影矩阵u3各点值做归一化处理，即，使得各点值的取值范围在0到1区间；然后将0到1之间等分为m个分区；计算该行中各点的数值落在每个分区的概率，令第x行中各点落在第q个分区的概率为p
xq
，其值为该行落在第q个分区的点数除以该行总点数，1≤q≤m，根据如下公式计算第x行的熵e
x
30.31.式中对数取2为底，或者取其它对数底，它们间可用换底公式换算。
32.进一步地，当所提取的该段心电的t波特征用于心电信号分类时，作为特征向量，训练分类器，然后进行分类。
33.尤其是，当所提取的特征用于检测该段心电信号中是否出现t波交替twa现象时：首先建立出现了twa的心电和未出现twa现象的心电这两类样本的数据集，每段心电所包含的t波个数大于等于8，然后对数据集中的心电提取所述t波特征，用来训练分类器；
34.待检测心电信号所包含的t波个数应与训练样本所包含的t波个数一致；
35.当分类器完成训练后，对待检测心电信号提取所述t波特征，输入到该分类器，由分类器判断是否出现了t波交替twa现象。
36.对比现有技术，本发明有益效果在于：将多个心拍的心电连续起来处理，采用分数阶傅里叶变换对每一心拍的t波数据进行扩维获得t波矩阵，将连续心拍的t波矩阵构建三阶张量，再通过tucker分解获得第三方向上的投影分量，取其熵值，作为该段心电t波的特征，从而携带了t波的波动信息，可结合机器学习用于进行心电信号分类，尤其是检测twa现象。
附图说明
37.图1为本发明所述方法的流程图；
38.图2为将n个t波进行张量化的示意图。
具体实施方式
39.下面将对本发明加以详细说明，同时也叙述了本发明技术方案解决的技术问题及有益效果，需要指出的是，所描述的实例仅旨在便于对本发明的理解，而对其不起任何限定作用。
40.下面以进行twa现象检测为例，结合附图说明本发明的具体实施方式。算法流程图见图1。使用的数据库包括physionet：t
‑
wave alternans challenge database数据库和mit
‑
bih normal sinus rhythm database数据库。
41.从physionet：t
‑
wave alternans challenge database数据库中取出编号为01,06,09,13,15,17,21,25,28,29,30,33,34,35,50,51,64,67,69,70,72,73,76,79,82,88,91,97,98共29个的含twa的心电信号数据，心电含有12个通道约2分钟的记录，每条数据的采样频率为500hz，约60000点，将每个数据中12个通道中的信号单独取出，则共有29*12＝348条含有twa现象的心电信号数据。
42.从mit
‑
bih normal sinus rhythm database数据库中取出编号为16265,16272,16273,16420,16483,16539,16773,16786,16795,17052,17453,18177,18184,19088,19090,19093,19140,19830共18个数据，每个数据含有2个通道，每条数据的采样频率为128hz，将每个数据中的2个通道信号单独取出，将每条数据进行截取得到10条含60000点的数据。于是共获得18*2*10＝360条正常心电信号数据。由于含twa现象的心电信号的采样频率为500hz,为了使正常心电信号采样频率和含twa现象的心电信号采样频率相同，对正常心电信号进行升采样到500hz。实施例中对心电样本的使用依需要还进行了进一步切分。
43.步骤s1：连续采集包含p个心动周期的心电信号并进行预处理，获得信号d；该段心
电所包含的心动周期个数p应大于等于8；
44.所述预处理包括去除心电信号中的工频干扰、肌电干扰和基线漂移。预处理还包括将心电信号升采样到统一的采样率。
45.作为优选，在这里，采用fir带通滤波器去除ecg信号中50hz工频干扰、肌电干扰及基线漂移，滤波器截止频率设为5hz和15hz。
46.步骤s2：对步骤s1预处理后的信号d进行心电的qrs波群定位，获得每个qrs波群的起点和止点及该qrs波群的r波峰值的位置信息。qrs波群定位的算法很多，例如可采用经典的pam
‑
tompkins算法实现检测和定位qrs波群及r波峰值位置。
47.步骤s3：构造t波矩阵
48.s3
‑
1：在每两个相邻的qrs波群之间，即在其中第一个qrs波群的止点和第二个qrs波群的起点之间，定位t波的起点和终点，并计算该段t波起点和终点之间的采样点的个数x；然后在该t波的起点和终点之间搜索极大值，以该极大值作为此t波峰值；
49.定位t波的起点和终点有多种方法，如cn201810187726.9、cn201910437430.2、佘黎煌.心电图中t波电交替检测关键技术研究[d].2015.等。
[0050]
令所检测到的t波的个数为n，n≤p；令所检测到的各t波相应的x的最小值，为最小采样点数x
min
；
[0051]
对检测到的每一个t波波峰，以t波波峰所在位置为中心采样点，分别向该t波峰值左侧和右侧各再取n个采样点，加上该t波波峰共2n+1个采样点，依时间顺序排列，构成该t波的2n+1点长的t波向量；实施例中，n取30，即t波向量包括61个采样点。
[0052]
所述点数2n+1应小于等于t波起点和终点之间的最小采样点数x
min
；
[0053]
实施例对348条含有twa现象的心电信号数据和360条正常心电信号数据进行预处理和t波定位后，分别取31个t波采样点和61个t波采样点。
[0054]
s3
‑
2：对每个t波的2n+1点长的t波向量，分别进行m个不同阶次的分数阶傅里叶变换，第j个t波向量经过i阶的分数阶傅里叶变换后得到2n+1点长的信号d
ji
；其中j≤n，所用到的分数阶傅里叶变换的阶次0≤i≤1，并且对所有t波，所做的m个阶次分数傅里叶变换的组合是相同的；
[0055]
分数阶傅里叶变换的公式如下：
[0056][0057]
其中
[0058][0059]
其中α＝pπ/2，p为分数阶傅里叶变换的阶次(范围0～1)，f
p
表示分数阶傅里叶变换算子。
[0060]
s3
‑
3：计算各信号d
ji
的分数阶傅里叶幅度谱信号fv
ji
；进一步地，作为优选：计算各信号d
ji
的分数阶傅里叶幅度谱信号fv
ji
之后，对每个t波，对每个分数阶傅里叶幅度谱信号fv
ji
以该阶次i下最大值最小值为参照，进行归一化处理，使每个fv
ji
的各点取值在0到1区间；得到归一化后的2n+1点长的分数阶傅里叶幅度谱信号。
[0061]
s3
‑
4：对每个t波，将m个不同阶次的分数阶傅里叶幅度谱信号作为行向量，按照所用的分数傅里叶变换的阶次从小到大依次从上到下排列，构成m行
×
(2n+1)列的t波矩阵；步骤s3
‑
1共检测到n个t波因此有n个这样的矩阵；见附图2。
[0062]
实施例中用20个阶次下的分数阶傅里叶变换构造t波矩阵。
[0063]
步骤s4：构造t波张量空间
[0064]
将n个t波矩阵按t波出现的顺序依次排列，得到m
×
(2n+1)
×
n的三阶张量；见附图2。
[0065]
步骤s5：张量分解
[0066]
对上述三阶张量采用tucker分解的方法进行张量分解，获得核心张量g和三个投影方向上的投影矩阵，对应分数阶傅里叶变换阶次方向的投影矩阵u1、t波向量方向对应阶的投影矩阵u2，和心拍方向对应的投影矩阵u3；
[0067]
步骤s6：
[0068]
对心拍方向对应的投影矩阵u3，求熵值，以此作为该段心电的t波特征。对投影矩阵u3每一列求均值，获得一维的行向量，对该行向量求熵值，获得该投影矩阵的熵值。所述熵可以是香农熵、renyi熵、近似熵、样本熵等中的某一种熵。本实施例采用香农熵。
[0069]
或者，步骤s6所述熵值是一个熵向量，计算方法为：
[0070]
对投影矩阵u3每一行计算熵值，将获得的熵值依所在行号依次排列，获得该投影矩阵的熵向量。投影矩阵u3每一行的熵值计算方法如下：
[0071]
将投影矩阵u3各点值做归一化处理，即，使得各点值的取值范围在0到1区间；然后将0到1之间等分为m个分区；计算该行中各点的数值落在每个分区的概率，令第x行中各点落在第q个分区的概率为p
xq
，其值为该行落在第q个分区的点数除以该行总点数，1≤q≤m，根据如下公式计算第x行的熵e
x
[0072][0073]
式中对数取2为底，或者取其它对数底，它们间可用换底公式换算。也可按照上述思路对投影矩阵u3每一列计算熵值，将获得的熵值依所在的列号一次排列，获得该投影矩阵的熵向量。
[0074]
当所提取的特征用于检测该段心电信号中是否出现t波交替twa现象时：首先建立出现了twa的心电和未出现twa现象的心电这两类样本的数据集，每段心电所包含的t波个数大于等于8，然后对数据集中的心电提取所述t波特征，用来训练分类器；待检测心电信号所包含的t波个数应与训练样本所包含的t波个数一致；当分类器完成训练后，对待检测心电信号提取所述t波特征，输入到该分类器，由分类器判断是否出现了t波交替twa现象。
[0075]
由于在这个实例中，希望检测一段心电中是否存在twa现象，需要解决的是一个二分类问题，所以实施例采用了简单通用并且效率高的svm算法(当数据集样本量增大，采用其他机器学习算法建模可进一步提升分类效果)，结果如下：
[0076]
当所用到的20个分数阶傅里叶变换的阶次取不同范围时，每个样本心电长度为2分钟，twa的检测效果如下：
[0077][0078]
当构建三阶张量所用的t波个数分别为8、30、60时，20个分数分数阶傅里叶变换阶次范围取0.5
‑
1时的twa检测效果如下：
[0079][0080]
本方法结合机器学习算法在优化参数组合后可以进一步提升性能，也可以用于其他类型心电的识别和分类。
[0081]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换和替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。