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

   这学期选了神经网络的课程，最后作业是处理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 savedHEADERFILE= strcat(stringname,'.hea');      % header-file in text formatATRFILE= strcat(stringname,'.atr');        % attributes-file in binary formatDATAFILE=strcat(stringname,'.dat');        % data-fileSAMPLES2READ=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 signalssfreq=A(2);   % sample rate of dataclear 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 212fid2=fopen(signald,'r');A= fread(fid2, [3, SAMPLES2READ], 'uint8')';  % matrix with 3 rows, each 8 bits long, = 2*12bitfclose(fid2);M2H= bitshift(A(:,2), -4);M1H= bitand(A(:,2), 15);PRL=bitshift(bitand(A(:,2),8),9);     % sign-bitPRR=bitshift(bitand(A(:,2),128),5);   % sign-bitM( : , 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 nosigcase 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 datafid3=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 = EOFclear 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=ECGsignalM1;figure(2);plot(ecgdata(1:points));grid on ;axis tight;axis([1,points,-2,5]);title('原始ECG信号');%%%%%%%%%%进行小波变换8层[C,L]=wavedec(ecgdata,level,wavename);%%%%%%%提取尺度系数，A1=appcoef(C,L,wavename,1);A2=appcoef(C,L,wavename,2);A3=appcoef(C,L,wavename,3);A4=appcoef(C,L,wavename,4);A5=appcoef(C,L,wavename,5);A6=appcoef(C,L,wavename,6);A7=appcoef(C,L,wavename,7);A8=appcoef(C,L,wavename,8);%%%%%%%提取细节系数D1=detcoef(C,L,1);D2=detcoef(C,L,2);D3=detcoef(C,L,3);D4=detcoef(C,L,4);D5=detcoef(C,L,5);D6=detcoef(C,L,6);D7=detcoef(C,L,7);D8=detcoef(C,L,8);%%%%%%%%%%%%重构A8=zeros(length(A8),1); %去除基线漂移,8层低频信息RA7=idwt(A8,D8,wavename);RA6=idwt(RA7(1:length(D7)),D7,wavename);RA5=idwt(RA6(1:length(D6)),D6,wavename);RA4=idwt(RA5(1:length(D5)),D5,wavename);RA3=idwt(RA4(1:length(D4)),D4,wavename);RA2=idwt(RA3(1:length(D3)),D3,wavename);D2=zeros(length(D2),1); %去除高频噪声，2层高频噪声RA1=idwt(RA2(1:length(D2)),D2,wavename);D1=zeros(length(D1),1);%去除高频噪声，1层高频噪声DenoisingSignal=idwt(RA1,D1,wavename);figure(3);plot(DenoisingSignal);title('去除噪声的ECG信号'); grid on; axis tight;axis([1,points,-2,5]);clear ecgdata;

去噪前后对比图像如下：

去噪前：

去噪后：

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.matmydata = DenoisingSignal;ecgdata=mydata';swa=zeros(4,points);%存储概貌信息swd=zeros(4,points);%存储细节信息signal=ecgdata(0*points+1:1*points); %取点信号%算小波系数和尺度系数%低通滤波器 1/4 3/4 3/4 1/4%高通滤波器 -1/4 -3/4 3/4 1/4%二进样条小波for i=1:points-3   swa(1,i+3)=1/4*signal(i+3-2^0*0)+3/4*signal(i+3-2^0*1)+3/4*signal(i+3-2^0*2)+1/4*signal(i+3-2^0*3);   swd(1,i+3)=-1/4*signal(i+3-2^0*0)-3/4*signal(i+3-2^0*1)+3/4*signal(i+3-2^0*2)+1/4*signal(i+3-2^0*3);endj=2;while j<=level   for i=1:points-24     swa(j,i+24)=1/4*swa(j-1,i+24-2^(j-1)*0)+3/4*swa(j-1,i+24-2^(j-1)*1)+3/4*swa(j-1,i+24-2^(j-1)*2)+1/4*swa(j-1,i+24-2^(j-1)*3);     swd(j,i+24)=-1/4*swa(j-1,i+24-2^(j-1)*0)-3/4*swa(j-1,i+24-2^(j-1)*1)+3/4*swa(j-1,i+24-2^(j-1)*2)+1/4*swa(j-1,i+24-2^(j-1)*3);   end   j=j+1;end%画出原信号和尺度系数，小波系数%figure(10);%subplot(level+1,1,1);plot(ecgdata(1:points));grid on ;axis tight;%title('ECG信号在j=1,2,3,4尺度下的尺度系数对比');%for i=1:level%    subplot(level+1,1,i+1);%    plot(swa(i,:));axis tight;grid on; xlabel('time');ylabel(strcat('a  ',num2str(i)));%end%figure(11);%subplot(level,1,1); plot(ecgdata(1:points)); grid on;axis tight;%title('ECG信号及其在j=1,2,3,4尺度下的尺度系数及小波系数');%for i=1:level%    subplot(level+1,2,2*(i)+1);%    plot(swa(i,:)); axis tight;grid on;xlabel('time');%    ylabel(strcat('a   ',num2str(i)));%    subplot(level+1,2,2*(i)+2);%    plot(swd(i,:)); axis tight;grid on;%    ylabel(strcat('d   ',num2str(i)));%end%画出原图及小波系数%figure(12);%subplot(level,1,1); plot(real(ecgdata(1:points)),'b'); grid on;axis tight;%title('ECG信号及其在j=1,2,3,4尺度下的小波系数');%for i=1:level%    subplot(level+1,1,i+1);%    plot(swd(i,:),'b'); axis tight;grid on;%    ylabel(strcat('d   ',num2str(i)));%end%**************************************求正负极大值对**********************%ddw=zeros(size(swd));pddw=ddw;nddw=ddw;%小波系数的大于0的点posw=swd.*(swd>0);%斜率大于0pdw=((posw(:,1:points-1)-posw(:,2:points))<0);%正极大值点pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);%小波系数小于0的点negw=swd.*(swd<0);ndw=((negw(:,1:points-1)-negw(:,2:points))>0);%负极大值点nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);%或运算ddw=pddw|nddw;ddw(:,1)=1;ddw(:,points)=1;%求出极值点的值,其他点置0wpeak=ddw.*swd;wpeak(:,1)=wpeak(:,1)+1e-10;wpeak(:,points)=wpeak(:,points)+1e-10;%画出各尺度下极值点%figure(13);%for i=1:level%    subplot(level,1,i);%    plot(wpeak(i,:)); axis tight;grid on;%ylabel(strcat('j=   ',num2str(i)));%end%subplot(4,1,1);%title('ECG信号在j=1,2,3,4尺度下的小波系数的模极大值点');interva2=zeros(1,points);intervaqs=zeros(1,points);Mj1=wpeak(1,:);Mj3=wpeak(3,:);Mj4=wpeak(4,:);%画出尺度3极值点figure(14);plot (Mj3);%title('尺度3下小波系数的模极大值点');posi=Mj3.*(Mj3>0);%求正极大值的平均thposi=(max(posi(1:round(points/4)))+max(posi(round(points/4):2*round(points/4)))+max(posi(2*round(points/4):3*round(points/4)))+max(posi(3*round(points/4):4*round(points/4))))/4;posi=(posi>thposi/3);nega=Mj3.*(Mj3<0);%求负极大值的平均thnega=(min(nega(1:round(points/4)))+min(nega(round(points/4):2*round(points/4)))+min(nega(2*round(points/4):3*round(points/4)))+min(nega(3*round(points/4):4*round(points/4))))/4;nega=-1*(nega<thnega/4);%找出非0点interva=posi+nega;loca=find(interva);for i=1:length(loca)-1    if abs(loca(i)-loca(i+1))<80       diff(i)=interva(loca(i))-interva(loca(i+1));    else       diff(i)=0;    endend%找出极值对loca2=find(diff==-2);%负极大值点interva2(loca(loca2(1:length(loca2))))=interva(loca(loca2(1:length(loca2))));%正极大值点interva2(loca(loca2(1:length(loca2))+1))=interva(loca(loca2(1:length(loca2))+1));intervaqs(1:points-10)=interva2(11:points);countR=zeros(1,1);countQ=zeros(1,1);countS=zeros(1,1);mark1=0;mark2=0;mark3=0;i=1;j=1;Rnum=0;%*************************求正负极值对过零，即R波峰值，并检测出QRS波起点及终点*******************%while i<points    if interva2(i)==-1       mark1=i;       i=i+1;       while(i<points&interva2(i)==0)          i=i+1;       end       mark2=i;%求极大值对的过零点       mark3= round((abs(Mj3(mark2))*mark1+mark2*abs(Mj3(mark1)))/(abs(Mj3(mark2))+abs(Mj3(mark1))));%R波极大值点       R_result(j)=mark3-10;%为何-10？经验值吧       countR(mark3-10)=1;%求出QRS波起点       kqs=mark3-10;       markq=0;     while (kqs>1)&&( markq< 3)         if Mj1(kqs)~=0            markq=markq+1;         end         kqs= kqs -1;     end  countQ(kqs)=-1;  %求出QRS波终点    kqs=mark3-10;  marks=0;  while (kqs<points)&&( marks<3)      if Mj1(kqs)~=0         marks=marks+1;      end      kqs= kqs+1;  end  countS(kqs)=-1;  i=i+60;  j=j+1;  Rnum=Rnum+1; endi=i+1;end%************************删除多检点，补偿漏检点**************************%num2=1;while(num2~=0)   num2=0;%j=3,过零点   R=find(countR);%过零点间隔   R_R=R(2:length(R))-R(1:length(R)-1);   RRmean=mean(R_R);%当两个R波间隔小于0.4RRmean时,去掉值小的R波for i=2:length(R)    if (R(i)-R(i-1))<=0.4*RRmean        num2=num2+1;        if signal(R(i))>signal(R(i-1))            countR(R(i-1))=0;        else            countR(R(i))=0;        end    endendendnum1=2;while(num1>0)   num1=num1-1;   R=find(countR);   R_R=R(2:length(R))-R(1:length(R)-1);   RRmean=mean(R_R);%当发现R波间隔大于1.6RRmean时,减小阈值,在这一段检测R波for i=2:length(R)    if (R(i)-R(i-1))>1.6*RRmean        Mjadjust=wpeak(4,R(i-1)+80:R(i)-80);        points2=(R(i)-80)-(R(i-1)+80)+1;%求正极大值点        adjustposi=Mjadjust.*(Mjadjust>0);        adjustposi=(adjustposi>thposi/4);%求负极大值点        adjustnega=Mjadjust.*(Mjadjust<0);        adjustnega=-1*(adjustnega<thnega/5);%或运算        interva4=adjustposi+adjustnega;%找出非0点        loca3=find(interva4);        diff2=interva4(loca3(1:length(loca3)-1))-interva4(loca3(2:length(loca3)));%如果有极大值对,找出极大值对        loca4=find(diff2==-2);        interva3=zeros(points2,1)';        for j=1:length(loca4)           interva3(loca3(loca4(j)))=interva4(loca3(loca4(j)));           interva3(loca3(loca4(j)+1))=interva4(loca3(loca4(j)+1));        end        mark4=0;        mark5=0;        mark6=0;    while j<points2         if interva3(j)==-1;            mark4=j;            j=j+1;            while(j<points2&interva3(j)==0)                 j=j+1;            end            mark5=j;%求过零点            mark6= round((abs(Mjadjust(mark5))*mark4+mark5*abs(Mjadjust(mark4)))/(abs(Mjadjust(mark5))+abs(Mjadjust(mark4))));            countR(R(i-1)+80+mark6-10)=1;            j=j+60;         end         j=j+1;     end    end endend%画出原图及标出检测结果%%%%%%%%%%%%%%%%%%%%%%%%%%开始求PT波段%对R波点前的波用加窗法，窗口大小为100，然后计算窗口内极大极小的距离%figure(20);%plot(Mj4);%title('j=4 细节系数'); hold on%%%%%%%还是直接求j=4时的R过零点吧Mj4posi=Mj4.*(Mj4>0);%求正极大值的平均Mj4thposi=(max(Mj4posi(1:round(points/4)))+max(Mj4posi(round(points/4):2*round(points/4)))+max(Mj4posi(2*round(points/4):3*round(points/4)))+max(Mj4posi(3*round(points/4):4*round(points/4))))/4;Mj4posi=(Mj4posi>Mj4thposi/3);Mj4nega=Mj4.*(Mj4<0);%求负极大值的平均Mj4thnega=(min(Mj4nega(1:round(points/4)))+min(Mj4nega(round(points/4):2*round(points/4)))+min(Mj4nega(2*round(points/4):3*round(points/4)))+min(Mj4nega(3*round(points/4):4*round(points/4))))/4;Mj4nega=-1*(Mj4nega<Mj4thnega/4);Mj4interval=Mj4posi+Mj4nega;Mj4local=find(Mj4interval);Mj4interva2=zeros(1,points);for i=1:length(Mj4local)-1    if abs(Mj4local(i)-Mj4local(i+1))<80       Mj4diff(i)=Mj4interval(Mj4local(i))-Mj4interval(Mj4local(i+1));    else       Mj4diff(i)=0;    endend%找出极值对Mj4local2=find(Mj4diff==-2);%负极大值点Mj4interva2(Mj4local(Mj4local2(1:length(Mj4local2))))=Mj4interval(Mj4local(Mj4local2(1:length(Mj4local2))));%正极大值点Mj4interva2(Mj4local(Mj4local2(1:length(Mj4local2))+1))=Mj4interval(Mj4local(Mj4local2(1:length(Mj4local2))+1));mark1=0;mark2=0;mark3=0;Mj4countR=zeros(1,1);Mj4countQ=zeros(1,1);Mj4countS=zeros(1,1);flag=0;while i<points    if Mj4interva2(i)==-1       mark1=i;       i=i+1;       while(i<points&Mj4interva2(i)==0)          i=i+1;       end       mark2=i;%求极大值对的过零点,在R4中极值之间过零点就是R点。       mark3= round((abs(Mj4(mark2))*mark1+mark2*abs(Mj4(mark1)))/(abs(Mj4(mark2))+abs(Mj4(mark1))));       Mj4countR(mark3)=1;       Mj4countQ(mark1)=-1;       Mj4countS(mark2)=-1;       flag=1;    end    if flag==1        i=i+200;        flag=0;    else        i=i+1;    endend%%%%%%%%%%%%%%%%%%%%%%%%找到MJ4的QRS点后，这里缺少对R点的漏点检测和冗余检测，先不去细究了。%%%%%%%%%%对尺度4下R点检测不够好，需要改进的地方%%%%%%%figure(200);%plot(Mj4);%title('j=4');%hold on;%plot(Mj4countR,'r');%plot(Mj4countQ,'g');%plot(Mj4countS,'g');%%%%%%%%%%%%%%%%%%%%%%%%%%Mj4过零点找到%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Rlocated=find(Mj4countR);Qlocated=find(Mj4countQ);Slocated=find(Mj4countS);countPMj4=zeros(1,1);countTMj4=zeros(1,1);countP=zeros(1,1);countPbegin=zeros(1,1);countPend=zeros(1,1);countT=zeros(1,1);countTbegin = zeros(1,1);countTend = zeros(1,1);windowSize=100;%%%%%%%%%%%%%%%%%%%%%%%P波检测%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Rlocated Qlocated 是在尺度4下的坐标for i=2:length(Rlocated)    flag=0;    mark4=0;    RRinteral=Rlocated(i)-Rlocated(i-1);    for j=1:5:(RRinteral*2/3)       % windowEnd=Rlocated(i)-30-j;       windowEnd=Qlocated(i)-j;        windowBegin=windowEnd-windowSize;        if windowBegin<Rlocated(i-1)+RRinteral/3            break;        end        %求窗内的极大极小值        %windowposi=Mj4.*(Mj4>0);        %windowthposi=(max(Mj4(windowBegin:windowBegin+windowSize/2))+max(Mj4(windowBegin+windowSize/2+1:windowEnd)))/2;        [windowMax,maxindex]=max(Mj4(windowBegin:windowEnd));        [windowMin,minindex]=min(Mj4(windowBegin:windowEnd));        if minindex < maxindex &&((maxindex-minindex)<windowSize*2/3)&&windowMax>0.01&&windowMin<-0.1            flag=1;            mark4=round((maxindex+minindex)/2+windowBegin);            countPMj4(mark4)=1;            countP(mark4-20)=1;            countPbegin(windowBegin+minindex-20)=-1;            countPend(windowBegin+maxindex-20)=-1;        end        if flag==1            break;         end    end    if mark4==0&&flag==0 %如果没有P波，在R波左间隔1/3处赋值-1        mark4=round(Rlocated(i)-RRinteral/3);        countP(mark4-20)=-1;    endend %plot(countPMj4,'g');       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%T波检测%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear windowBegin windowEnd maxindex minindex windowMax windowMin mark4 RRinteral;windowSizeQ=100;for i=1:length(Rlocated)-1;    flag=0;    mark5=0;    RRinteral=Rlocated(i+1)-Rlocated(i);    for j=1:5:(RRinteral*2/3)       % windowBegin=Rlocated(i)+30+j;       windowBegin=Slocated(i)+j;        windowEnd  =windowBegin+windowSizeQ;        if windowEnd >Rlocated(i+1)-RRinteral/4            break;        end        %%%%%求窗口内的极大极小值        [windowMax,maxindex]=max(Mj4(windowBegin:windowEnd));        [windowMin,minindex]=min(Mj4(windowBegin:windowEnd));        if minindex < maxindex &&((maxindex-minindex)<windowSizeQ)&&windowMax>0.1&&windowMin<-0.1            flag=1;            mark5=round((maxindex+minindex)/2+windowBegin);            countTMj4(mark5)=1;            countT(mark5-20)=1;%找到T波峰值点            %%%%%确定T波起始点和终点            countTbegin(windowBegin+minindex-20)=-1;            countTend(windowBegin+maxindex-20)=-1;        end        if flag==1            break;        end    end    if mark5==0 %如果没有T波，在R波右 间隔1/3处赋值-2        mark5=round(Rlocated(i)+ RRinteral/3);        countT(mark5)=-2;    endend%plot(countTMj4,'g');%hold off;        figure(4);plot(ecgdata(0*points+1:1*points)),grid on,axis tight,axis([1,points,-2,5]);title('ECG信号的各波波段检测');hold onplot(countR,'r');plot(countQ,'k');plot(countS,'k');for i=1:Rnum    if R_result(i)==0;        break    end    plot(R_result(i),ecgdata(R_result(i)),'bo','MarkerSize',10,'MarkerEdgeColor','g');endplot(countP,'r');plot(countT,'r');plot(countPbegin,'k');plot(countPend,'k');plot(countTbegin,'k');plot(countTend,'k');hold off



4特征提取
将各波段的位置提取出来后，我们根据15个距离特征与6个幅值特征作为身份识别的特征。具体信息简下表：
距离特征：
R-QR-SR-PP-PBR-PER-TR-TBR-TEPB-PETB-TEQ-PS-TP-TQ-PBS-TE幅值特征：
Q-RS-RPB-PP-QT-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; endsave testsample111.mat  testECGsave 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;    endenddisp('正确率：');correct/length(t)
运行结果：正确率为 0.96 左右，效果还不错。



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



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






勉励自己一句话：勤学如春起之苗，不见其长，日有所赠；

                          辍学如磨刀之石，不见其损，日有所亏。

                                       奋斗吧--碗。