rddata.m

可以分别读取MIT-BIH心电图库的头文件

% 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: [e-mail]),
%
%      algorithm is based on a program written by
%      Klaus Rheinberger (University of Innsbruck)
%      (email: [e-mail])
%
% -------------------------------------------------------------------------
clc; clear all;

%------ SPECIFY DATA ------------------------------------------------------
PATH= 'D:\matlab701\work'; % path, where data are saved
HEADERFILE= '118e18.txt';      % header-file in text format
ATRFILE= '118e18.atr';         % attributes-file in binary format
DATAFILE='118e18.dat';         % data-file
SAMPLES2READ=50000;         % 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=[];
    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 A 是三列3000行的矩阵
fclose(fid2);
M2H= bitshift(A(:,2), -4);        %字节向右移四位，取每一行中的第二列的字节的高四位，即取字节的高四位，比如56(111000)，右移四位就是(00000011),为3
M1H= bitand(A(:,2), 15);          %取已经移动了的字节的低四位，就是原来字节的低四位，bit(a,b)a,b的二进制对应位求与运算，比如bitand(3(0011),10(1010)),之后为(0010)，即为2
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;%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;%time是3000个采样点的时间位置
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')';%A有2311行，2列
fclose(fid3);
ATRTIME=[];
ANNOT=[];%保存每个心跳的注释代码
sa=size(A);%sa为[2311，2]
saa=sa(1);%saa为2311
i=1;
while i<=saa  %i=1:2311
    annoth=bitshift(A(i,2),-2);%查看第一个储存单元的高六位，实际阅读的顺序因该是A(i,2)A(i,1)
    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+2列的第一个字节，第i+2列的第二个字节的
        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));%find函数返回的是下标，是注释文件中的时间长度要小于你所采集的时间长度，%
                                %比如，你采集了3000个点，时间到8.3秒，那么ATRTIME的时间就要小于8.3秒
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.7,num2str(ANNOTD(k)));%ATRTIME相当于x轴，0是y轴，num2str是将其转化成字符
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');