matlab_采用差分阈值法对心电信号QRS波进行特征提取ecg0

PATH= 'C:\MATLAB7\work\ecg0'; %指定数据的储存路径 HEADERFILE= '100.hea'; %.hea 格式，头文件，可用记事本打开 DATAFILE='100.dat'; %.dat 格式，ECG 数据 SAMPLES2READ=3600; %指定需要读入的样本数,若.dat文件中存储有两个通道的信号:则读入 2*SAMPLES2READ 个数据 %读取头文件中的信息 %fprintf(1,'\\n$> WORKING ON %s ...\n', HEADERFILE); % 在Matlab命令行窗口提示当前工作状态 signalh=fullfile(PATH, HEADERFILE); % 通过函数 fullfile 获得头文件的完整路径 fid1=fopen(signalh,'r'); % 打开头文件，其标识符为 fid1 ，属性为'r'--“只读” z=fgetl(fid1); % 读取头文件的第一行数据，字符串格式 A=sscanf(z, '%*s %d %d %d',[1,3]); % 按照格式 '%*s %d %d %d' 转换数据并存入矩阵 A 中 nosig= A(1); % 信号通道数目 sfreq=A(2); % 数据采样频率 clear A; % 清空矩阵 A ，准备获取下一行数据 for k=1:nosig % 读取每个通道信号的数据信息 z= fgetl(fid1); A= sscanf(z, '%*s %d %d %d %d %d',[1,5]); dformat(k)= A(1); % 信号格式; 这里只允许为 212 格式 gain(k)= A(2); % 每 mV 包含的整数个数 bitres(k)= A(3); % 采样精度（位分辨率） zerovalue(k)= A(4); % ECG 信号零点相应的整数值 firstvalue(k)= A(5); % 信号的第一个整数值 (用于偏差测试) end; fclose(fid1); clear A; %读取data数据 if dformat~= [212,212], error('this script does not apply binary formats different to 212.'); end; signald= fullfile(PATH, DATAFILE); % 读入 212 格式的 ECG 信号数据 fid2=fopen(signald,'r'); A= fread(fid2, [3,SAMPLES2READ], 'uint8'); % matrix with 3 rows, each 8 bits long, = 2*12bit 矩阵A共有SAMPLES2READ行、3列，每列数据都是以uint8格式读入，注意这时数据通过uint8的读入方式已经成为十进制数了 fclose(fid2); M2H= bitshift(A(2,:), -4); % 字节向右移四位，即取字节的高四位，属于信号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;% 将M1H、M2H分别左移8位，即乘以2^8，再分别加上A(:,1)，A(:,2)， M( 2 , :)= bitshift(M2H,8)+ A( 2 , : )-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 转置为2行SAMPLES2READ列 %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'); %画出心电信号波形 X=M(1,:); figure(1); clf, box on, hold on plot(TIME,X,'r'); %if nosig==2 %plot(TIME, M(2,:),'b'); %end; %for k=1:length(ATRTIMED) %text(ATRTIMED(k),0,num2str(ANNOTD(k))); %end; plot(TIME,0,'g'); %标出x轴 xlim([TIME(1), TIME(end)]); ylim([-1.1, 1.1]); 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'); %EMG滤波第一步，LPF1??有必要吗？幅值减小 %for i=3:2048 % X(i)=(X0(i)+3*X0(i-1)+X0(i-2))/5; %end %X(1)=X0(1); %X(2)=X0(2); %figure(2); %plot(TIME,X); %hold on %plot(TIME,X0,'r'); %plot(TIME,0,'g'); %标出x轴 %xlim([TIME(1), TIME(end)]); %ylim([-1.1, 1.1]); %确定QRS波群 for n=2:2047 Y1(n)=X(n+1)-X(n-1); %一阶导数 end for n=3:2046 Y2(n)=X(n+2)-2*X(n)+X(n-2); %二阶导数 end max1=max(Y1); max2=max(Y2); for n=3:2046 Y3(n)=max1*Y1(n)+max2*Y2(n); end a=0.1; max3=max(Y3); QRS=[]; for i=3:2046 if Y3(i)>max3*a Y4(i)=1; QRS=[QRS,i]; else Y4(i)=0; end end figure(3); plot(TIME,X,'r'); hold on plot(TIME(3:2046),Y4(3:2046)); plot(TIME,0,'g'); xlim([TIME(1), TIME(end)]); ylim([-1.1, 1.1]); %对Y4的处理1，去毛刺 for i=3:2046 if (Y4(i)==1&&Y4(i-1)==0&&Y4(i+1)==0) Y4(i)=0; end end figure(4) plot(TIME,X,'r'); hold on plot(TIME(3:2046),Y4(3:2046)); plot(TIME,0,'g'); xlim([TIME(1), TIME(end)]); ylim([-1.1, 1.1]); %对Y4的处理2，使QRS波群范围内一阶导数全为1，无“零点”，最终确定QRS波群 L=480; %窗口大小，根据两个R峰之间的点数来确定 for i=3:1950 if (Y4(i)==0&&Y4(i-1)==1) for j=1:L/5 if Y4(i+j)==1 Y4(i)=1; break; end end end end figure(5) plot(TIME,X,'r'); hold on plot(TIME(3:2046),Y4(3:2046)); plot(TIME,0,'g'); xlim([TIME(1), TIME(end)]); ylim([-1.1, 1.1]); %确定R峰 B=[]; for i=3:2046 %找出QRS波群的起始点 if (Y4(i)==1&&Y4(i-1)==0) B=[B,i]; elseif (Y4(i)==1&&Y4(i+1)==0) B=[B,i]; end end s=size(B); W=[]; MAXR=[]; for j=1:2:s(2)-2 a=B(j); b=B(j+1); [maxR,w0]=max(X(1,a:b)); %找出QRS波群内幅值最大的一点，即为R峰 w=a+w0-1; %得出的最大值的位置是以a开头向后数（w0-1）位 W=[W,w]; MAXR=[MAXR,maxR]; end [maxRL,w0L]=max(X(1,B(s(2)-1):2048)); %最后一个QRS波群是不完整的，再对Y4进行处理2时，算不到最后的98个点 wL=B(s(2)-1)+w0L-1; W=[W,wL]; MAXR=[MAXR,maxRL]; %figure(6) %plot(TIME,X,'r'); %hold on %plot(TIME(W),X(W),'o'); %plot(TIME,0,'g'); %xlim([TIME(1), TIME(end)]); %ylim([-1.1, 1.1]); %求峰域值 dsum=[]; for n=2:2047 d=Y1(n).^2; dsum=[dsum,d]; end dert=sum(dsum)/2046; %方差 se=26*dert; figure(7) plot(TIME(2:2047),Y1(2:2047));%画出一阶导数的波形 hold on plot(TIME,X,'k'); plot(TIME,se,'r'); plot(TIME,-se,'r'); plot(TIME,0,'g'); xlim([TIME(1), TIME(end)]); %R峰起点 Rnum=size(W); Rstart=[]; for i=1:Rnum(2) s1=0; for j=(W(i)-1):-1:1 if abs(Y1(j))>se s1=j; %R峰左边斜率首先超过se线的点为起点s1 break; end end for m=s1:-1:1 if (abs(Y1(m))<se) %第一个斜率穿过se线的即为R峰起点 Rstart=[Rstart,m]; break; end end end figure(8) plot(TIME,X,'r'); hold on plot(TIME(Rstart),X(Rstart),'o') plot(TIME,0,'g'); xlim([TIME(1), TIME(end)]); ylim([-1.1, 1.1]); %Q峰位置 RSnum=size(Rstart); Q=[]; for i=1:RSnum(2) for j=Rstart(i):-1:1 if (Rstart(i)-j)<0.04*sfreq %由R峰起点向前0.04s内（因为Q波的宽度小于0.04s） if (abs(X(j))>=abs(X(j-1))&&abs(X(j))>=abs(X(j+1))) %第一个绝对值极大值即为Q峰 Q=[Q,j]; break; end end end end %figure(9) %plot(TIME,X,'r'); %hold on %plot(TIME(Q),X(Q),'o') %plot(TIME,0,'g'); %xlim([TIME(1), TIME(end)]); %ylim([-1.1, 1.1]); %R峰终点 Rend=[]; for i=1:Rnum(2)-1 s2=0; c1=0; c2=0; for j=(W(i)+1):2047 if abs(Y1(j))>se %R峰右边斜率首先超过se线的点为起点s2 s2=j; break; end end %for c=(s2+1):2047 % if(abs(Y1(c))>abs(Y1(c+1))&&abs(Y1(c))>abs(Y1(c-1))) % c1=c; % break; % end % end for k=(s2+1):2047 if (abs(Y1(k))>se&&abs(Y1(k+1))<se) %第一个自下而上穿越-se线的即为R峰终点 c2=k;