这学期的课程选择神经网络。最后的作业处理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=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.mat mydata = 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); end j=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); %斜率大于0 pdw=((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; %求出极值点的值,其它点置0 wpeak=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; end end %找出极值对 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; end i=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 end end end num1=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 end end %画出原图及标出检測结果 %%%%%%%%%%%%%%%%%%%%%%%%%%開始求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; end end %找出极值对 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; end end %%%%%%%%%%%%%%%%%%%%%%%%找到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; end end %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; end end %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 on plot(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'); end plot(countP,'r'); plot(countT,'r'); plot(countPbegin,'k'); plot(countPend,'k'); plot(countTbegin,'k'); plot(countTend,'k'); hold off4特征提取
将各波段的位置提取出来后，我们依据15个距离特征与6个幅值特征作为身份识别的特征。详细信息简下表：
距离特征：
| R-Q | R-S | R-P |
| P-PB | R-PE | R-T |
| R-TB | R-TE | PB-PE |
| TB-TE | Q-P | S-T |
| P-T | Q-PB | S-TE |
幅值特征：
| Q-R | S-R |
| PB-P | P-Q |
| T-TB | T-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
希望大家给予批评。有错误之处务必指正。最后感谢能够坚持看到最后的人们！
勉励自己一句话：勤学如春起之苗，不见其长。日有所赠；
辍学如磨刀之石，不见其损，日有所亏。
奋斗吧--碗。
版权声明：本文博客原创文章，博客，未经同意，不得转载。
时间： 2024-11-07 08:52:05