MIT-BIH心率失常数据提取及部分MATLAB程序解释

首先，在PhysioBank ATM中下载心率失常数据（网址：http://write.blog.csdn.net/postedit），第一步，在database中选择

第二步，点击，在接下来的页面中就可看到48组数据。

本文，使用100号数据进行数据的读取和代码的分析。

点击下载到你想得文件夹里。其中100.hea文件，单击右键下载，左键点击你可以直接看见里面的数据。

下文是matlab代码：

%—— SPECIFY DATA ——————————————————
PATH= ‘E:\MIT_BIHdata’; %这里需要改为你的数据所存的位置
HEADERFILE= ‘100.hea’;      % header-file in text format
ATRFILE= ‘100.atr’;         % attributes-file in binary format
DATAFILE=’100.dat’;         % data-file
SAMPLES2READ=3000;         % 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);%输出\n> WORKING ON+HEADERFILE的值
signalh= fullfile(PATH, HEADERFILE);
fid1=fopen(signalh,’r’);
z= fgetl(fid1);
A= sscanf(z, ‘%*s %d %d %d’,[1,3]);%*s为不存入
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;%如果不是212格式就返回error
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);%获得A的行列
saa=sa(1);%获得sa的一个数
i=1;

%以下代码为读取标注的一种译码
while i<=saa
    annoth=bitshift(A(i,2),-2);%取高6位
    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)];%第一个16位变化+第二个16位变化总和为发送时间。
        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);%正shift为左移。
        hilfe=hilfe+mod(hilfe,2);%判断奇偶，奇数则跳过一个16位
        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;%cunsum各行相加
ind= find(ATRTIME <= TIME(end));%找到小于TIME的值
ATRTIMED= ATRTIME(ind);
ANNOT=round(ANNOT);
ANNOTD= ANNOT(ind);
%—— DISPLAY DATA ——————————————————
figure(1); clf, box on, hold 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)));%图中标记ATRTIMED里的值
end;
xlim([TIME(1), TIME(end)]);%设置x轴上下线
xlabel(‘Time / s’); ylabel(‘Voltage / mV’);
string=[‘ECG signal ‘,DATAFILE];
title(string);
fprintf(1,’\\n >DISPLAYINGDATAFINISHED\n′);fprintf(1,′n > ALL FINISHED \n’);

该代码核心只需知道，我们如何读取心率失常的数据，然后如何将它的标注显示出来，这为我们后续的工作做好准备。标注译码的过程可以参考该文章

http://wenku.baidu.com/link?url=ithj7n-yModPue-HfBixyXkNAO094Q7YPUuSVba_KJrKojn_VJL5uY7q8Mg3FA2s87cUpvPvbsRikHdzbTjUNwQHY7sivakq1qLAw49Pigy

以下为MATLAB运行结果：

其中，标注为1的是正常心率，而8和28为异常。具体数字对应何异常参考http://wenku.baidu.com/view/9a61f02fbd64783e09122bcf.html

最后解释以下，ECG48个 信号介绍：1）有23个记录是从一个 包含4000个24小时动态ECG记录的 数据集中随机选取的，这些动态ECG记录是波士顿Beth Israel 医院从住院病人（约60%）和门诊病人（约40%）这样一个混合人群里收集的。2）其余25个记录是从同一个数据集中选取的，用来包含那些不常见的但是临床上有重要意义的在一个小随机样本中不能很好的表示的心律失常。