multiscalemm.m

数学形态学用于技术分析ECG心律的探测

clc; clear all;

%------ SPECIFY DATA ------------------------------------------------------
PATH= 'E:\ecg\mit-cd'; % path, where data are saved
HEADERFILE= '100.hea';      % header-file in text format
ATRFILE= '100.atr';         % attributes-file in binary format
DATAFILE='100.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);

%------ DISPLAY DATA ------------------------------------------------------
figure(1); clf, box on, hold on
plot(TIME, M(:,1),'r');
figure;
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');




figure;
% clc;
% close all;
% clear;
% load 101.mat;
D101=M;
D101(:,2)=M(:,1);
x=1:1:648000;
y=M(:,1);
subplot(5,1,1);
plot(y);
axis([0 21600 -1 2]);
% above original ECG signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y3=y;
for v=1:1:3
D101(:,2)=y3;
n=0;
m=0;
E101=D101;
F101=D101;
ES=[1,1,1];
IS=[5,5,5];
for n=1:1:647997
    for m=1:1:3
    IS(m)=F101(n+m,2)-ES(m);
    end
    E101(n,2)=min(IS);
end
for n=647998:1:647999
  for m=1:1:3-(n-647997)
    IS(m)=F101(n+m,2)-ES(m);
    end
    E101(n,2)=min(IS);
end  
y0=E101(:,2);
subplot(5,1,2);
plot(x,y0);
axis([0 647999 -1 3]);

n=0;
m=0;
G101=E101;
H101=E101;
ES=[1,1,1];
IS=[-5,-5,-5];
for n=2:1:3
   for m=1:1:n-1
    IS(m)=H101(n-m,2)+ES(m);
    end
    G101(n,2)=max(IS);
end 
for n=4:1:648000
    for m=1:1:3
    IS(m)=H101(n-m,2)+ES(m);
    end
    G101(n,2)=max(IS);
end
y1=G101(:,2);
subplot(5,1,3);
plot(x,y1);
axis([0 647999 -1 0.8]);

n=0;
m=0;
E101=G101;
F101=G101;
ES=[1,1,1];
IS=[-5,-5,-5];
for n=2:1:3
   for m=1:1:n-1
    IS(m)=F101(n-m,2)+ES(m);
    end
    E101(n,2)=max(IS);
end 
for n=4:1:648000
    for m=1:1:3
    IS(m)=F101(n-m,2)+ES(m);
    end
    E101(n,2)=max(IS);
end


y2=E101(:,2);
subplot(5,1,4);
plot(x,y2);
axis([0 647999 -1 3]);

n=0;
m=0;
G101=E101;
H101=E101;
ES=[1,1,1];
IS=[5,5,5];
for n=1:1:647997
    for m=1:1:3
    IS(m)=H101(n+m,2)-ES(m);
    end
    G101(n,2)=min(IS);
end
ES=[1,1,1];
IS=[5,5,5];
for n=647998:1:647999
  for m=1:1:3-(n-647997)
    IS(m)=H101(n+m,2)-ES(m);
    end
    G101(n,2)=min(IS);
end  
y3=G101(:,2);

end
Subplot(5,1,5);
plot(x,y3);
axis([0 647999 -1 0.8]);
fprintf(1,'\\n$> ALL FINISHED \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y7=y;
for u=1:1:3
D101(:,2)=y7;
figure;
subplot(5,1,1);
plot(x,y);
axis([0 647999 -1 2]);

n=0;
m=0;
E101=D101;
F101=D101;
ES=[1,1,1];
IS=[-5,-5,-5];
for n=2:1:3
   for m=1:1:n-1
    IS(m)=F101(n-m,2)+ES(m);
    end
    E101(n,2)=max(IS);
end 
for n=4:1:648000
    for m=1:1:3
    IS(m)=F101(n-m,2)+ES(m);
    end
    E101(n,2)=max(IS);
end
y4=E101(:,2);
subplot(5,1,2);
plot(x,y4);
axis([0 647999 -1 3]);

n=0;
m=0;
G101=E101;
H101=E101;
ES=[1,1,1];
IS=[5,5,5];
for n=1:1:647997
    for m=1:1:3
    IS(m)=H101(n+m,2)-ES(m);
    end
    G101(n,2)=min(IS);
end
ES=[1,1,1];
IS=[5,5,5];
for n=647998:1:647999
  for m=1:1:3-(n-647997)
    IS(m)=H101(n+m,2)-ES(m);
    end
    G101(n,2)=min(IS);
end  

end
y5=G101(:,2);
Subplot(5,1,3);
plot(x,y5);
axis([0 647999 -1 0.8]);

n=0;
m=0;
E101=G101;
F101=G101;
ES=[1,1,1];
IS=[5,5,5];
for n=1:1:647997
    for m=1:1:3
    IS(m)=F101(n+m,2)-ES(m);
    end
    E101(n,2)=min(IS);
end
for n=647998:1:647999
  for m=1:1:3-(n-21597)
    IS(m)=F101(n+m,2)-ES(m);
    end
    E101(n,2)=min(IS);
end  
y6=E101(:,2);
subplot(5,1,4);
plot(x,y6);
axis([0 647999 -1 3]);

n=0;
m=0;
G101=E101;
H101=E101;
ES=[1,1,1];
IS=[-5,-5,-5];
for n=2:1:3
   for m=1:1:n-1
    IS(m)=H101(n-m,2)+ES(m);
    end
    G101(n,2)=max(IS);
end 
for n=4:1:648000
    for m=1:1:3
    IS(m)=H101(n-m,2)+ES(m);
    end
    G101(n,2)=max(IS);
end
y7=G101(:,2);
Subplot(5,1,5);
plot(x,y7);
axis([0 647999 -1 0.8]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\\n$> ALL FINISHED \n');
figure;
y8=(y3+y7)/2;
Subplot(5,1,1);
plot(x,y8);
axis([0 647999 -1 0.8]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S101=D101;
S101(:,2)=y8;
n=0;
m=0;
E101=S101;
F101=S101;
ES=[0,0.1,0.2,0.4,0.6,0.8,1,0.8,0.6,0.4,0.2,0.1,0];
IS=[5,5,5,5,5,5,5,5,5,5,5,5,5];
for n=1:1:647987
    for m=1:1:13
    IS(m)=F101(n+m,2)-ES(m);
    end
    E101(n,2)=min(IS);
end
IS=[5,5,5,5,5,5,5,5,5,5,5,5,5];
for n=647988:1:647999
    for m=1:1:13-(n-647987)
    IS(m)=F101(n+m,2)-ES(m);
    end
    E101(n,2)=min(IS);
end
y9=E101(:,2);
subplot(5,1,2);
plot(x,y9);
axis([0 647999 -1 3]);

n=0;
m=0;
G101=E101;
H101=E101;
ES=[0,0.1,0.2,0.4,0.6,0.8,1,0.8,0.6,0.4,0.2,0.1,0];
IS=[-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5];
for n=2:1:13
    for m=1:1:n-1
    IS(m)=H101(n-m,2)+ES(m);
    end
    G101(n,2)=max(IS);
end
IS=[-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5];
for n=14:1:648000
    for m=1:1:13
    IS(m)=H101(n-m,2)+ES(m);
    end
    G101(n,2)=max(IS);
end
y10=G101(:,2);
Subplot(5,1,3);
plot(x,y10);
axis([0 647999 -1 0.8]);

n=0;
m=0;
E101=G101;
F101=G101;
ES=[0,0.1,0.2,0.4,0.6,0.8,1,0.8,0.6,0.4,0.2,0.1,0];
IS=[-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5];
for n=2:1:13
    for m=1:1:n-1
    IS(m)=F101(n-m,2)+ES(m);
    end
    E101(n,2)=max(IS);
end
IS=[-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5];
for n=14:1:648000
    for m=1:1:13
    IS(m)=F101(n-m,2)+ES(m);
    end
    E101(n,2)=max(IS);
end
y11=E101(:,2);
subplot(5,1,4);
plot(x,y11);
axis([0 4 -1 3]);

n=0;
m=0;
G101=E101;
H101=E101;
ES=[0,0.1,0.2,0.4,0.6,0.8,1,0.8,0.6,0.4,0.2,0.1,0];
IS=[5,5,5,5,5,5,5,5,5,5,5,5,5];
for n=1:1:647987
    for m=1:1:13
    IS(m)=H101(n+m,2)-ES(m);
    end
    G101(n,2)=min(IS);
end
IS=[5,5,5,5,5,5,5,5,5,5,5,5,5];
for n=647988:1:647999
    for m=1:1:13-(n-647987)
    IS(m)=H101(n+m,2)-ES(m);
    end
    G101(n,2)=min(IS);
end
y12=G101(:,2);
Subplot(5,1,5);
% figure;
plot(x,y12);
axis([0 647999 -1 1]);
fprintf(1,'\\n$> ALL FINISHED \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
yy=y8-y12;
plot(x,yy);
axis([390000 400000 -1 1.5]);
n=0;
number1=0;
H101(:,2)=yy;
y13=yy;
for n=1:1:20
    y13(n)=0;
end
for n=647951:1:648000
    y13(n)=0;
end
% y13=0;
% for n=20:1:647950;
%     y13(n)=H101(n,2);
% end
number1=0;
T1=0.05*max(y13);
T2=0.1*max(y13);
k=-100;
for n=20:1:648000%7201%14401%
if H101(n,2)>=T1&(n-k)>150
number1=number1+1;
k=n;
end%if%
end %for%

HR=number1/30
number1 %#ok<NOPRT>
fprintf(1,'\\n$> ALL FINISHED \n');

%Deciosion rules   
    
% figure;
% subplot(4,1,1)
% plot(x,y);
% axis([0 20 -1 2]);
% subplot(4,1,2)
% plot(x,y8);
% axis([0 20 -1 2]);
% subplot(4,1,3)
% plot(x,y12);
% axis([0 20 -1 2]);
% subplot(4,1,4)
% plot(x,yy);
% axis([0 60 -1 2]);