本发明涉及医学信号处理技术领域,更确切地说一种心电图特征向量提取方法。
背景技术：
与传统的基于密码或id卡的认证方式相比，生物识别技术(biometrics)具有更高的安全度与可靠度，已逐渐成为国际研究的热点，并逐渐在不同的领域发挥着举足轻重的作用。
近年来借助于人体内蕴的心电信号ecg(electrocardiogram)进行身份识别的方法广受关注。ecg是从人体体表采集的反映心脏心动的电位信号，人体的生理条件差异使得ecg具有许多个体特征。相较指纹、语音以及手掌，ecg作为一种活体生物信号，具备易检测、难复制的特点。利用心电信号进行身份识别的定义为：给定一条心电信号，判定该信号所属人的身份。利用心电信号进行自动身份识别的系统或装置通常基于心电信号的模式识别技术进行实现，需要实现心电信号的特征提取和识别两个必要环节。
在已公开的利用心电信号进行身份识别的装置和方法的发明中，在提取心电图特征向量时，特征向量往往通过一段较长时间的心电信号提取，存在两方面问题，一方面用户需要较长时间的等待进行心电信号采集，另一方面较大的数据量需要耗费更多的计算成本。因此关于心跳信号的识别方法开始兴起，此类方法的主要优势在于从心电信号中分割一些独立的心跳信号来提取特征向量，提高了心电信号身份识别的实际应用效果。
技术实现要素：
本发明的目的是为解决用心电图进行身份识别时间长、成本高的问题，而提供一种心电图特征向量提取方法。
一种心电图特征向量提取方法，它包括：
1）导入数据
从外部采集n个人的单导联心电信号，每人采集两条，n大于50，存储为数组变量s,共计2n条。对2n条心电信号的2n个标签存储为变量label_l；
2）去除心电信号的噪声
对数组变量s内每一个数组执行小波分解重构算法，结果存储为数组变量s1。
3）提取r波
在数组变量s1内每一个数组中找到m个r波峰值点位置，共计2n*m个位置，在每个位置前取99点，在每个位置后取100点，包含当前位置点共计取200点作为一条r波。2n*m个位置提取到2n*m条r波，依次放入2n*m行、200列的矩阵中，记该变量为y，作为构建特征方法的输入信号；对应变量y内2n*m条r波的2n*m个标签存储为变量label。
4）构建特征
a．将y中的r波划分训练集和测试集,把从n个人的第一条信号中提取的n*m条r波作为训练集，第二条信号中提取的n*m条r波，作为测试集，分别存放在n*m行、200列的矩阵中，记为train和test；从变量label中划分出训练集和测试集对应的标签，分别存放在长度为n*m的数组中，训练集标签记为train_label，测试集标签记为test_label；将训练集train和测试集test分别分割为n*m*20段，每一段长度为10,分别存放在n*m*20行、10列的矩阵中，记为train_f和test_f；
（2）对train_f执行k-means聚类，设置聚类簇数为50，聚类后得到大小为50*10的矩阵，记为c1，将c1转置得到10*50的矩阵，记为c，对train_f和c计算欧式距离，计算结果存放在n*m*20行、50列的矩阵中，记为train_sample_f；对test_f和c计算欧式距离，计算结果存放在n*m*20行、50列的矩阵中，记为test_sample_f；
（3）将训练集样本train_sample_f和测试集样本test_sample_f分别按照列优先的原则变维为n*m行、1000列的矩阵并归一化到[0,1]之间，记为train_sample和test_sample；
train_sample即为构建好的测试集特征，test_sample即为构建好的测试集特征；
5）对样本进行自动识别
将训练集特征train_sample进行稀疏，把稀疏后的训练集特征以及对应的标签train_label全部输入到lr分类器中进行训练，lr分类器训练结束后，保留参数，再将测试集特征test_sample进行稀疏，把稀疏后的测试集特征以及对应的标签test_label输入到训练好的lr分类器进行测试；
识别结果为判断r波所属个体的判断准确率a；
6）输出分类结果a；
所述小波分解重构算法为：
a．小波分解：选择harr函数作为小波基。数组变量s内每一个数组分解得到近似系数向量ca1,细节系数向量cd1；对ca1分解得到近似系数向量ca2,细节系数向量cd2；对ca2分解得到近似系数向量ca3,细节系数向量cd3；对ca3分解得到近似系数向量ca4,细节系数向量cd4；对ca4分解得到近似系数向量ca5,细节系数向量cd5；对ca5分解得到近似系数向量ca6,细节系数向量cd6；对ca6分解得到近似系数向量ca7,细节系数向量cd7；对ca7分解得到近似系数向量ca8,细节系数向量cd8；
b．去除噪声：对需要去噪的信号cd2，采用软阈值方式去噪。各层系数降噪所需要的阈值不同，采用stdc作为噪声强度估计；
stdc公式：stdc2=median(|cd2|)/0.6745，其中median(m)为计算m的平均值；
计算出每一层的阈值：
thr2=(2*log(length(cd2)))^(1/2)*stdc2，其中length(cd2)为计算细节向量cd2的长度；
对需要去噪的信号cd2，将cd2中模值小于3*thr2的系数置零，将模值大于3*thr2的系数做处理：大于3*thr2的系数统一减去3*thr2，小于-3*thr2的系数统一加3*thr2
计算后得到去噪后的信号cd2；
对cd3、cd4、cd5、cd6、cd7、cd8也采用上述软阈值方式去除噪声，得到cd3、cd4、cd5、cd6、cd7、cd8；
c．小波重构：对分解得到的系数向量ca8置零并赋给ca8，对cd8完成小波重构得到信号ca7；对cd7完成小波重构得到信号ca6；对cd6完成小波重构得到信号ca5；对cd5完成小波重构得到信号ca4；对ca4完成小波重构得到信号ca3；对cd3完成小波重构得到信号ca2；对cd2完成小波重构得到信号ca1；ca1即为s内一个数组执行小波分解重构算法的结果。
与现有技术相比，本发明的优势是：本发明截取r波顶点前后1.2秒信号进行识别，有效减小了身份识别需要的信号采集时间。本方法中对训练集采用k-means聚类，原理比较简单，实现很容易，收敛速度快。
具体实施方式
实施例1一种心电图特征向量提取方法
下面结合具体的实施方式对本发明作进一步说明。
具体实例中使用国际通行心电图数据库ecg-id,该数据库的数据及使用说明公开于行业内周知的physionet.org网站。ecg-id数据集包含310录音,从90个主题获得(44个男性和46位女性)。每一个记录含有一个导联的信号,记录在30秒的持续时间,并数字化在500赫兹与12位分辨率超过标称±10mv范围。在本实例中，共使用55个人的单导联心电信号，每人采集两条通过工作在计算机上的软件系统和行业内所周知的matlab仿真环境进行实现。
本实施例的详细步骤如下：
一．导入数据
从外部采集55个人的单导联心电信号，每人采集两条，存储为数组变量s,所述的数组变量s存储共计110条心电信号。对应数组变量s内110条心电信号的110个标签存储为变量label_l。
二．去除心电信号的噪声
对数组变量s内每一个数组执行小波分解重构算法，结果存储为数组变量s1，s1内每个数组对应s内一个数组执行小波分解重构算法的结果。所述小波分解重构算法为：
（1）小波分解
选择harr函数作为小波基。对数组变量s内每一个数组分解得到近似系数向量ca1,细节系数向量cd1；对ca1分解得到近似系数向量ca2,细节系数向量cd2；对ca2分解得到近似系数向量ca3,细节系数向量cd3；对ca3分解得到近似系数向量ca4,细节系数向量cd4；对ca4分解得到近似系数向量ca5,细节系数向量cd5；对ca5分解得到近似系数向量ca6,细节系数向量cd6；对ca6分解得到近似系数向量ca7,细节系数向量cd7；对ca7分解得到近似系数向量ca8,细节系数向量cd8。
（2）去除噪声
对需要去噪的信号cd2，采用软阈值方式去噪。各层系数降噪所需要的阈值不同，采用stdc作为噪声强度估计
stdc公式：stdc2=median(|cd2|)/0.6745，其中median(m)为计算m的平均值
计算出每一层的阈值：
thr2=(2*log(length(cd2)))^(1/2)*stdc2，其中length(cd2)为计算细节向量cd2的长度
对需要去噪的信号cd2，将cd2中模小于3*thr2的系数全部置零，而将模大于3*thr2的做一个比较特殊的处理，大于3*thr2的系数统一减去3*thr2，小于-3*thr2的系数统一加3*thr2
计算后得到去噪后的信号cd2。
对cd3、cd4、cd5、cd6、cd7、cd8也采用上述软阈值方式去除噪声。得到cd3、cd4、cd5、cd6、cd7、cd8。
（3）小波重构
对分解得到的系数向量ca8直接置零并赋给ca8，对ca8,cd8完成小波重构得到信号ca7；对ca7,cd7完成小波重构得到信号ca6；对ca6,cd6完成小波重构得到信号ca5；对ca5,cd5完成小波重构得到信号ca4；对ca4,cd4完成小波重构得到信号ca3；对ca3,cd3完成小波重构得到信号ca2；对ca2,cd2完成小波重构得到信号ca1。ca1即为s内一个数组执行小波分解重构算法的结果。
三．提取r波
（1）在数组变量s1内每一个数组中找到13个r波峰值点位置（心电图数据库ecg-id有个别单导联心电图信号只有13个r波峰值点），共计1430个位置。在每个位置前取99点，在每个位置后取100点，包含当前位置点共计取200点作为一条r波，1430个位置提取到1430条r波，将这些r波依次放入1430行、200列的矩阵中，命名该变量为y。
（2）y将会作为构建特征方法的输入信号。对应变量y内1430条r波的1430个标签存储为变量label。
四．构建特征
所述构建特征方法的输入信号为y，y内存储着从110条心电信号中提取的1430条r波。
（1）将y中的r波划分训练集和测试集,将从每个人第一条信号中提取的715条r波作为训练集，存放在715行、200列的矩阵中并命名为train，将从每个人第二条信号中提取的715条r波，存放在715行、200列的矩阵中并命名为test。从变量label中划分出训练集和测试集对应的标签，得到训练集标签存放在长度为715的数组中，命名为train_label。得到测试集标签存放在长度为715的数组中，命名为test_label。将训练集train分割为14300段，每一段长度为10,存放在14300行、10列的矩阵中，命名为train_f；将测试集test分割为14300段，每一段长度为10,存放在14300行、10列的矩阵中，命命名为test_f。
（2）对train_f进行k-means聚类，聚类数为50，聚类后得到大小为50*10的矩阵，命名为c1。将c1转置得到10*50的矩阵，命名为c，对train_f和c计算欧式距离.计算结果存放在14300行、50列的矩阵中，命名为train_sample_f；对test_f和c计算欧式距离.计算结果存放在14300行、50列的矩阵中，命名为test_sample_f。
在matlab中kmeans的用法为：
[c1,~]=kmeans(train_f,50)
在matlab中dist函数计算欧式距离：
train_sample_f=dist(train_f,c);
（3）将训练集样本train_sample_f按照列优先的原则变维为715行、1000列的矩阵并归一化（在matlab中用mapminmax函数实现，下同）到[0,1]之间，命名为train_sample；将测试集样本train_sample_f按照列优先的原则变维为715行、1000列的矩阵并归一化到[0,1]之间，命名为test_sample。
train_sample即为构建好的测试集特征，test_sample即为构建好的测试集特征。
五．对样本进行自动识别
将训练集特征train_sample进行稀疏，把稀疏后的训练集特征以及对应的标签train_label全部输入到lr分类器中进行训练，lr分类器训练结束后，保留参数，再将测试集特征test_sample进行稀疏，把稀疏后的测试集特征以及对应的标签test_label输入到训练好的lr分类器进行测试。
测试结束后网络会自动给出识别结果。识别结果为判断r波所属个体的判断准确率a。
六、输出分类结果a。
结果显示在表1中