這學期的課程選擇神經網絡。最后的作業處理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 off
4特征提取
將各波段的位置提取出來后，我們依據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 左右。效果還不錯。
希望大家給予批評。有錯誤之處務必指正。最后感謝能夠堅持看到最后的人們！
勉勵自己一句話： 勤學如春起之苗，不見其長。日有所贈；
輟學如磨刀之石，不見其損，日有所虧。
奮斗吧--碗。
版權聲明：本文博客原創文章，博客，未經同意，不得轉載。