心电信号预处理_滤波去噪_QRs波检测 P波_T波检测_matlab

% This programm reads ECG data which are saved in format 212. % (e.g., 100.dat from MIT-BIH-DB, cu01.dat from CU-DB,...) % The data are displayed in a figure together with the annotations. % The annotations are saved in the vector ANNOT, the corresponding % times (in seconds) are saved in the vector ATRTIME. % The annotations are saved as numbers, the meaning of the numbers can % be found in the codetable "ecgcodes.h" available at www.physionet.org. % % ANNOT only contains the most important information, which is displayed % with the program rdann (available on www.physionet.org) in the 3rd row. % The 4th to 6th row are not saved in ANNOT. % % % created on Feb. 27, 2003 by % Robert Tratnig (Vorarlberg University of Applied Sciences) % (email: rtratnig@gmx.at), % % algorithm is based on a program written by % Klaus Rheinberger (University of Innsbruck) % (email: klaus.rheinberger@uibk.ac.at) % % ------------------------------------------------------------------------- clc; clear all; %M2=load ('D:\Cjprg\MIT-BIH-ECG-Database\102_1.txt'); %------ SPECIFY DATA ------------------------------------------------------ %PATH= 'M:\Profile\Desktop2k\MIT_CD\mitdb'; % path, where data are saved PATH= 'E:\人工神经网络\心电波检测\data'; % path, where data are savedD:\Cjprg\MIT-BIH-ECG-Database HEADERFILE= '118.hea'; % header-file in text format ATRFILE= '118.atr'; % attributes-file in binary format DATAFILE='118.dat'; % data-file SAMPLES2READ=648000; % 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); % ------------------------------------------------------------------------- fprintf(1,'\\n$> ALL FINISHED \n'); %BP网络 %输入样本 N1=369; % N s1=M(N1-140:N1+140,1); [cA1,cD1]=dwt(s1,'db1'); %采用db1基本小波来分解信号 N2=21868; % N s2=M(N2-140:N2+140,1); [cA2,cD2]=dwt(s2,'db1'); N3=79980; % N s3=M(N3-140:N3+140,1); [cA3,cD3]=dwt(s3,'db1'); N4=38488; % PVC s4=M(N4-140:N4+140,1); [cA4,cD4]=dwt(s4,'db1'); N5=165956; % PVC s5=M(N5-140:N5+140,1); [cA5,cD5]=dwt(s5,'db1'); N6=489189; % PVC s6=M(N6-140:N6+140,1); [cA6,cD6]=dwt(s6,'db1'); N7=226001; % N s7=M(N7-140:N7+140,1); [cA7,cD7]=dwt(s7,'db1'); N8=344565; % N s8=M(N8-140:N8+140,1); [cA8,cD8]=dwt(s8,'db1'); N9=444078; % N s9=M(N9-140:N9+140,1); [cA9,cD9]=dwt(s9,'db1'); N10=542658; % N s10=M(N10-140:N10+140,1); [cA10,cD10]=dwt(s10,'db1'); %figure(1); %subplot(2,1,1);plot(cA1);title('Approximation A1'); %显示分解结果 %subplot(2,1,2);plot(cA4);title('Approximation A4'); %显示分解结果 %训练网络 P = [cA1 cA2 cA4 cA3 cA7 cA5 cA8 cA9 cA6 cA10]; T = [0 0 1 0 0 1 0 0 1 0]; pr(1:141,1)=-4; pr(1:141,2)=4; net = newff(pr,[55 1],{'tansig' 'purelin'},'traingdx','learngdm'); net.trainParam.epochs = 1000; net.trainParam.goal=0.0002; net.trainParam.lr=0.0003; net = train(net,P,T); %测试一段心电波形 T=0.2; Rnum=0; counter=0; flag=0; i=140; while i<648000 if M(i,1)>T T=M(i,1); flag=1; i=i+1; elseif flag==1 %检测到一个R波峰 max=T; T=0.2; Rnum=Rnum+1; s=M(i-141:i+139,1); [cA,cD]=dwt(s,'db1'); %小波变换 P = cA; Y = sim(net,P); %送入网络仿真 y=round(Y); if y==1 counter=counter+1; %PVC的个数 fprintf(1,'发生PVC的样本点：%d\n',i-1) end; i=i+140; flag=0; else i=i+1; end; end; fprintf(1,'心电波个数：%d\n',Rnum) fprintf(1,'PVC个数：%d\n',counter)