首页 > 代码库 > ECG信号读取，检测QRS，P，T 波（基于小波去噪与检测），基于BP神经网络的身份识别

ECG信号读取，检测QRS，P，T 波（基于小波去噪与检测），基于BP神经网络的身份识别

2024-07-09 01:24:57   215人阅读

   这学期选了神经网络的课程，最后作业是处理ECG信号，并利用神经网络进行识别。

1  ECG介绍与读取ECG信号

1）ECG介绍

   具体ECG背景应用就不介绍了，大家可以参考百度 谷歌。只是简单说下ECG的结构：


   一个完整周期的ECG信号有 QRS P T 波组成，不同的人对应不用的波形，同一个人在不同的阶段波形也不同。我们需要根据各个波形的特点，提取出相应的特征，对不同的人进行身份识别。

2）ECG信号读取

首先需要到MIT-BIH数据库中下载ECG信号，具体的下载地址与程序读取内容介绍可以参考一下地址（讲述的很详细）：http://blog.csdn.net/chenyusiyuan/article/details/2027887。
   读取代码（基于MATLAB）如下：
  
clc; clear all;
%------ SPECIFY DATA ------------------------------------------------------
%%选择文件名
stringname='111';
%选择你要处理的信号点数
points=10000; 
PATH= 'F:\ECG\MIT-BIH database directory'; % path, where data are saved
HEADERFILE= strcat(stringname,'.hea');      % header-file in text format
ATRFILE= strcat(stringname,'.atr');        % attributes-file in binary format
DATAFILE=strcat(stringname,'.dat');        % data-file
SAMPLES2READ=points;         % number of samples to be read
                      % in case of more than one signal:
                            % 2*SAMPLES2READ samples are read
   
%------ LOAD HEADER DATA --------------------------------------------------
fprintf(1,'\\n$> WORKING ON %s ...\n', HEADERFILE);
signalh= fullfile(PATH, HEADERFILE);
fid1=fopen(signalh,'r');
z= fgetl(fid1);
A= sscanf(z, '%*s %d %d %d',[1,3]);
nosig= A(1);  % number of signals
sfreq=A(2);   % sample rate of data
clear A;
for k=1:nosig
    z= fgetl(fid1);
    A= sscanf(z, '%*s %d %d %d %d %d',[1,5]);
    dformat(k)= A(1);           % format; here only 212 is allowed
    gain(k)= A(2);              % number of integers per mV
    bitres(k)= A(3);            % bitresolution
    zerovalue(k)= A(4);         % integer value of ECG zero point
    firstvalue(k)= A(5);        % first integer value of signal (to test for errors)
end;
fclose(fid1);
clear A;

%------ LOAD BINARY DATA --------------------------------------------------
if dformat~= [212,212], error('this script does not apply binary formats different to 212.'); end;
signald= fullfile(PATH, DATAFILE);            % data in format 212
fid2=fopen(signald,'r');
A= fread(fid2, [3, SAMPLES2READ], 'uint8')';  % matrix with 3 rows, each 8 bits long, = 2*12bit
fclose(fid2);
M2H= bitshift(A(:,2), -4);
M1H= bitand(A(:,2), 15);
PRL=bitshift(bitand(A(:,2),8),9);     % sign-bit
PRR=bitshift(bitand(A(:,2),128),5);   % sign-bit
M( : , 1)= bitshift(M1H,8)+ A(:,1)-PRL;
M( : , 2)= bitshift(M2H,8)+ A(:,3)-PRR;
if M(1,:) ~= firstvalue, error('inconsistency in the first bit values'); end;
switch nosig
case 2
    M( : , 1)= (M( : , 1)- zerovalue(1))/gain(1);
    M( : , 2)= (M( : , 2)- zerovalue(2))/gain(2);
    TIME=(0:(SAMPLES2READ-1))/sfreq;
case 1
    M( : , 1)= (M( : , 1)- zerovalue(1));
    M( : , 2)= (M( : , 2)- zerovalue(1));
    M=M';
    M(1)=[];
    sM=size(M);
    sM=sM(2)+1;
    M(sM)=0;
    M=M';
    M=M/gain(1);
    TIME=(0:2*(SAMPLES2READ)-1)/sfreq;
otherwise  % this case did not appear up to now!
    % here M has to be sorted!!!
    disp('Sorting algorithm for more than 2 signals not programmed yet!');
end;
clear A M1H M2H PRR PRL;
fprintf(1,'\\n$> LOADING DATA FINISHED \n');

%------ LOAD ATTRIBUTES DATA ----------------------------------------------
atrd= fullfile(PATH, ATRFILE);      % attribute file with annotation data
fid3=fopen(atrd,'r');
A= fread(fid3, [2, inf], 'uint8')';
fclose(fid3);
ATRTIME=[];
ANNOT=[];
sa=size(A);
saa=sa(1);
i=1;
while i<=saa
    annoth=bitshift(A(i,2),-2);
    if annoth==59
        ANNOT=[ANNOT;bitshift(A(i+3,2),-2)];
        ATRTIME=[ATRTIME;A(i+2,1)+bitshift(A(i+2,2),8)+...
                bitshift(A(i+1,1),16)+bitshift(A(i+1,2),24)];
        i=i+3;
    elseif annoth==60
        % nothing to do!
    elseif annoth==61
        % nothing to do!
    elseif annoth==62
        % nothing to do!
    elseif annoth==63
        hilfe=bitshift(bitand(A(i,2),3),8)+A(i,1);
        hilfe=hilfe+mod(hilfe,2);
        i=i+hilfe/2;
    else
        ATRTIME=[ATRTIME;bitshift(bitand(A(i,2),3),8)+A(i,1)];
        ANNOT=[ANNOT;bitshift(A(i,2),-2)];
   end;
   i=i+1;
end;
ANNOT(length(ANNOT))=[];       % last line = EOF (=0)
ATRTIME(length(ATRTIME))=[];   % last line = EOF
clear A;
ATRTIME= (cumsum(ATRTIME))/sfreq;
ind= find(ATRTIME <= TIME(end));
ATRTIMED= ATRTIME(ind);
ANNOT=round(ANNOT);
ANNOTD= ANNOT(ind);

%------ DISPLAY DATA ------------------------------------------------------
figure(1); clf, box on, hold on ;grid on ;
plot(TIME, M(:,1),'r');
if nosig==2
    plot(TIME, M(:,2),'b');
end;
for k=1:length(ATRTIMED)
    text(ATRTIMED(k),0,num2str(ANNOTD(k)));
end;
xlim([TIME(1), TIME(end)]);
xlabel('Time / s'); ylabel('Voltage / mV');
string=['ECG signal ',DATAFILE];
title(string);
fprintf(1,'\\n$> DISPLAYING DATA FINISHED \n');
% -------------------------------------------------------------------------
fprintf(1,'\\n$> ALL FINISHED \n');

以MIT-BIH数据库中111.dat 为例。

2 去除高频噪声与基线漂移

   ECG读取完后，原始ECG信号含有高频噪声和基线漂移，利用小波方法可以去除相应噪声。具体原理如下：将一维的ECG信号进行8层的小波分解后（MATLAB下wavedec函数，小波类型是bior2.6）得到相应的细节系数与近似系数。根据小波原理我们可以知道，1,2层的细节系数包含了大部分高频噪声，8层的近似系数包含了基线漂移。基于此，我们将1,2层的细节系数（即高频系数置0），8成的近似系数（低频系数）置0；在对应进行小波重构，重构后我们可以明显得到去噪信号，信号无基线漂移。下面通过图片与代码进一步讲解：
   小波去噪代码：（matlab） 
  
%%%%%%%%%%%%%%%%%%%去除噪声和基线漂移%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
level=8; wavename='bior2.6';
ecgdata=http://www.mamicode.com/ECGsignalM1;>
去噪前后对比图像如下：
去噪前：




去噪后：




3 QRS 检测
   QRS检测是处理ECG信号的基础，无论最后实现什么样的功能，QRS波的检测都是前提。所以准确的检测QRS波是特征提取的前提，我采用基于二进样条4层小波变换，在3层的细节系数中利用极大极小值方法可以很好的检测出R波。3层细节系数的选择是基于R波在3层系数下表现的与其他噪声差别最大；具体实现如下：
二进样条小波滤波器：  低通滤波器：[1/4 3/4 3/4 1/4]
                      高通滤波器：[-1/4 -3/4 3/4 1/4]
在第3层细节系数中首先找到极大极小值对：
1）找极大值方法：找出斜率大于0的值，并赋值为1，其余为0，极大值就在序列类似1， 0这样的点，即前面一个值比后面的大的值对应的位置点。
2）找极小值方法：类似极大值,找出斜率<0的值对应的位置，并赋值为1，其余的为0，极小值就在类似1,0的序列中对应的位置，即前面一个值比后面的大的值对应的位置点。
检测出的极大极小值如下：


3）设置阈值，提取出R波。我们可以看出，R波的值要明显大于其他位置的值，其在3层细节系数的特点也类似于此。这样我们就可以设置一个可靠的阈值（将所有点分为4部分，求出每部分最大值的平均值T，阈值为T/3）来提取一组相邻的最大最小值对。这样最大最小值间的过0点就是对应于原始信号的R波点。
R波对应的极大极小值对如下：




4）补偿R波点。由于在二进样条小波变换的过程中，3层细节系数与原始信号的对应的位置有10个点的漂移，在程序中需要补偿。（这个在程序中会给出）。
5）找Q S 波。基于R波的位置，在R波位置（在1层细节系数下）的前3个极点为Q波，在R波位置（1细节系数下）的后3个极点为S波。这样我们就将QRS波定位出来。
6）由于不同的情况，可能造成R波的漏检和错检（把T波检测为R波），我们根据相邻R波的距离进行检测漏检与错检。当相邻R波的距离<0.4 mean(RR)平均距离时，这是错检，这样去除值小的R波。当相邻R波的距离>1.6mean(RR)时，在两个RR波间找到一个最大的极值对，定位R波，这是防止漏检。


经过上述方法，一个鲁棒性很好的QRS检测方法就出来了，经过测试，QRS检测能达到98%。检测结果R波用红线标注，Q S 波用黑线标注。


4  T P 波检测
     P T 波的检测与R波检测有很大的相同性。只不过 P T 波在4层细节系数中可以表述出更好的特性。同样根据根据极大极小值原理，可以分别检测出T P波，以及他们的起始点与终止点，即TB,TE,PB PE。具体程序我会在稍后的程序中给出。

各波段检测结果如下：



具体QRS T P波检查代码如下：
<pre name="code" class="cpp">level=4;    sr=360; 
%读入ECG信号
%load ecgdata.mat;
%load ECGsignalM1.mat;
%load Rsignal.mat
mydata = http://www.mamicode.com/DenoisingSignal;>


4特征提取
将各波段的位置提取出来后，我们根据15个距离特征与6个幅值特征作为身份识别的特征。具体信息简下表：
距离特征：
R-QR-SR-P
P-PBR-PER-T
R-TBR-TEPB-PE
TB-TEQ-PS-T
P-TQ-PBS-TE
幅值特征：
Q-RS-R
PB-PP-Q
T-TBT-S


我们将MIT-BIH中的101.dat、103.dat、105.dat、106.dat、111.dat分别取出10个这样的特征，其中5个作为训练样本、5个作为测试样本，送入神经网络进行训练。

特征提取代码：
%%%%%%%%%%%%%%%%%%%%%%%%%提取特征向量，进行神经网络训练%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%特征向量根据你需要检测部位的不同，选取特征向量。
%%%%%%%%%%%%%%%本例进行身份识别，选取5组信号，即5个同的人，每组数据采取10例ECG信号，  
%%%%%%%%%%%%%%%提取每例的15个距离特征向量、6个幅值特征向量作为特征数据
%%%%%%%%%%%%%%%距离特征：R-Q R-S R-P R-PBegin R-PEnd R-T R-TBegin R-TEnd
%%%%%%%%%%%%%%% PBegin-PEnd TBegin-TEnd Q-P S-T P-T Q-PBegin S-TEnd
%%%%%%%%%%%%%%%幅值特征： Q-R S-R PBegin-P P-Q T-TBegin T-S
%%%%%%%%%%%%%%每组的10例信号中5个训练5个测试
%%%%%%%%%%%%%%10组信号取第 2 4 6 8 10 12 14 16 18 20个周期， 2 6 10 14 18训练，其余测试


%%%%首先找到R Q S P T峰值， 起点 终点 的位置
locatedR=find(countR);
locatedQ=find(countQ);
locatedS=find(countS);
locatedP=find(countP);
locatedPBegin=find(countPbegin);
locatedPEnd=find(countPend);
locatedTBegin=find(countTbegin);
locatedTEnd=find(countTend);
locatedT=find(countT);
%%%%%%初始化各种特征值
RQ=0;RS=0;RP=0;RPB=0;RPE=0;RT=0;RTB=0;RTE=0;
PBPE=0;TBTE=0;QP=0;ST=0;PT=0;QPB=0;STE=0;
ampQR=0;ampSR=0;ampPBP=0;ampPQ=0;ampTTB=0;ampTS=0;
testECG=zeros(5,21);
counttest=1;
trainECG=zeros(5,21);
counttrain=1;
%%%%%%%%%%%%%%%%%开始计算
for i=2:2:20
    %距离特征
    RQ=abs(locatedR(i)-locatedQ(i));
    RS=abs(locatedS(i)-locatedR(i));
    RP=abs(locatedR(i)-locatedP(i-1));
    RPB=abs(locatedR(i)-locatedPBegin(i-1));
    RPE=abs(locatedR(i)-locatedPEnd(i-1));
    RT=abs(locatedR(i)-locatedT(i));
    RTB=abs(locatedR(i)-locatedTBegin(i));
    RTE=abs(locatedR(i)-locatedTEnd(i));
    PBPE=abs(locatedPBegin(i-1)-locatedPEnd(i-1));
    TBTE=abs(locatedTBegin(i)-locatedTEnd(i));
    QP=abs(locatedQ(i)-locatedP(i-1));
    ST=abs(locatedS(i)-locatedT(i));
    PT=abs(locatedP(i-1)-locatedT(i));
    QPB=abs(locatedQ(i)-locatedPBegin(i-1));
    STE=abs(locatedS(i)-locatedTEnd(i));
    %幅值特征
    ampQR=ecgdata(locatedR(i))-ecgdata(locatedQ(i));
    ampSR=ecgdata(locatedR(i))-ecgdata(locatedS(i));
    ampPBP=ecgdata(locatedP(i-1))-ecgdata(locatedPBegin(i-1));
    ampPQ=ecgdata(locatedQ(i))-ecgdata(locatedP(i-1));
    ampTTB=ecgdata(locatedT(i))-ecgdata(locatedTBegin(i));
    ampTS=ecgdata(locatedT(i))-ecgdata(locatedS(i));
    %%%%组成向量，并归一化
    featureVector=[RQ,RS,RP,RPB,RPE,RT,RTB,RTE,PBPE,TBTE,QP,ST,PT,QPB,STE];
    maxFeature=max(featureVector);
    minFeature=min(featureVector);
    for j=1:length(featureVector)
        featureVector(j)=2*(featureVector(j)-minFeature)/(maxFeature-minFeature)-1;
    end
    amplitudeVector=[ampQR,ampSR,ampPBP,ampPQ,ampTTB,ampTS];
    maxAmplitude=max(amplitudeVector);
    minAmplitued=min(amplitudeVector);
    for j=1:length(amplitudeVector)
        amplitudeVector(j)=2*(amplitudeVector(j)-minAmplitued)/(maxAmplitude-minAmplitued)-1;
    end
    if rem(i,4)==0
        testECG(counttest,:)=[featureVector,amplitudeVector];
        counttest=counttest+1;
    else
        trainECG(counttrain,:)=[featureVector,amplitudeVector];
        counttrain=counttrain+1;
    end
    clear amplitudeVector  featureVector; 
end
save testsample111.mat  testECG
save trainsample111.mat trainECG


5 BP神经网络训练与检测
我相信很多人对神经网络比较熟悉了，这里我就不多讲了，在matlab中，主要有三个函数， newff 负责建立网络， train 负责训练网络， sim 负责进行仿真。调整好参数。就可以进行训练与测试啦。

具体代码如下：


clear all;
load testsample101.mat;
test101=testECG;
load testsample103.mat;
test103=testECG;
load testsample105.mat;
test105=testECG;
load testsample106.mat;
test106=testECG;
load testsample111.mat;
test111=testECG;

load trainsample101.mat;
train101=trainECG;
load trainsample103.mat;
train103=trainECG;
load trainsample105.mat;
train105=trainECG;
load trainsample106.mat;
train106=trainECG;
load trainsample111.mat;
train111=trainECG;
%训练样本
train_sample=[ train103' train101' train105' train106' train111']; %21*25
%测试样本
test_sample=[test103' test101' test105' test106' test111'];
%输出类别
t=[2 2 2 2 2 1 1 1 1 1 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5];
train_result=ind2vec(t);
test_result=ind2vec(t);
pr(1:21,1)=-1;
pr(1:21,2)=1;
net = newff(pr,[21,5],{'tansig' 'purelin'},'traingdx','learngdm');
net.trainParam.epochs=1000;
net.trainParam.goal=0.0002;
net.trainParam.lr=0.0003;
net = train(net,train_sample,train_result);
result_sim=sim(net,test_sample);
result_sim_ind=vec2ind(result_sim);
correct=0;
for i=1:length(t)
    if result_sim_ind(i)==t(i);
        correct=correct+1;
    end
end
disp('正确率：');correct/length(t)




运行结果：正确率为 0.96 左右，效果还不错。



6：本次ECG实现的所有代码与相关原理信息的下载地址（0积分）：http://download.csdn.net/detail/yuansanwan123/7530687



希望大家给予批评。有错误之处务必指正。最后感谢能够坚持看到最后的人们！



勉励自己一句话：勤学如春起之苗，不见其长，日有所赠；
                辍学如磨刀之石，不见其损，日有所亏。
                                       奋斗吧--碗。



































        

            
                
                  
                
                

                    联系

                    我们