function [M,annotsig,heasig] = load_physio(fname,samples,offset)
% funkce nacita signaly ulozene ve formatu 212 (viz www.physionet.org).
%
% [M,annotsig,heasig]=load_physio(fname); - nacteni celeho souboru
% [M,annotsig,heasig]=load_physio(fname,samples); - nacteni daneho poctu vzorku od zacatku souboru
% [M,annotsig,heasig]=load_physio(fname,samples,offset); - nacteni daneho poctu vzorku od stanoveneho vzorku
%
% VYSTUPY:
% M ............ 1. sloupec je cas, 2. a 3. jsou signaly nactene ze souboru *.dat
% annotsig ..... anotace nactene ze souboru *.atr
% heasig ....... informace o signalu nactene ze souboru *.hea
%
% VSTUPY:
% fname ... nazev souboru (vcetne cesty a pripony) cesta/*.dat 
%           je mozne zadata kterykoliv ze tri souboru signalu (atr,dat,hea)
% samples .. delka vystupniho vektoru (pocet vzorku)
%           / (Inf - read to the end of the file)
% offset .. cislo vzorku od ktereho se signal nacte
%
%
% PRIKLAD POUZITI:
%     [M,annotsig,heasig] = load_physio('x:/mypath/myfile.dat',inf,0)
%     [M,annotsig,heasig] = load_physio('x:/mypath/myfile.hea')
%

% 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 ANNOTD.
% 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.
%
%      2010-03-21, Jan Hrubes, DBME, FEEC, VUT Brno
%      program is script modified to function
%      
%      script was 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)
%
% -------------------------------------------------------------------------

%------ INPUTS ------------------------------------------------------
if ~exist('offset','var')
    offset = 0;
end

fname = fname(1:end-3); % cut off extension

%------ LOAD HEADER DATA --------------------------------------------------
% fprintf(1,'\\n$> WORKING ON %s ...\n', [fname '*']);
fid1=fopen([fname 'hea'],'r');
z= fgetl(fid1);
heasig.name = char(sscanf(z, '%s %*d %*d %*d',[1,inf]));
A= sscanf(z, '%*s %d %d %d',[1,3]);
heasig.nsig = A(1);  % number of signals
heasig.freq = A(2);   % sample rate of data
heasig.nsamp = A(3);   % number of samples

clear A;
for k=1:heasig.nsig
    z= fgetl(fid1);
    A= sscanf(z, '%*s %d %d %d %d %d %d %d %*s',[1,7]);
    heasig.fmt(k) = A(1);           % format; here only 212 is allowed
    heasig.gain(k) = A(2);              % number of integers per mV
    heasig.bitres(k) = A(3);            % bitresolution
    heasig.zerovalue(k) = A(4);         % integer value of ECG zero point
    heasig.firstvalue(k) = A(5);        % first integer value of signal (to test for errors)
    heasig.cksum(k) = A(6);
    heasig.bsize(k) = A(7);
    heasig.lead{k} = char(sscanf(z, '%*s %*d %*d %*d %*d %*d %*d %*d %s',[1,inf]));
    heasig.sigfile{k} = char(sscanf(z, '%s %*d %*d %*d %*d %*d %*d %*d %*s',[1,inf]));
end;
status = fclose(fid1);
% clear A z fid1



%------ LOAD BINARY DATA --------------------------------------------------
if ~exist('samples','var')
    samples = heasig.nsamp-offset; % default loads whole signal
end
if (isinf(samples)) || ((samples+offset) > heasig.nsamp)
    samples = heasig.nsamp-offset;
end

SAMPLES2READ=samples+offset;
if heasig.fmt~= [212,212], error('this script does not apply binary formats different to 212.'); end;
fid2=fopen([fname 'dat'],'r');
A= fread(fid2, [3, SAMPLES2READ], 'uint8')';  % matrix with 3 rows, each 8 bits long, = 2*12bit
status = 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,:) ~= heasig.firstvalue, error('inconsistency in the first bit values'); end;
switch heasig.nsig
case 2
    M( : , 1)= (M( : , 1)- heasig.zerovalue(1))/heasig.gain(1);
    M( : , 2)= (M( : , 2)- heasig.zerovalue(2))/heasig.gain(2);
    TIME=(0:(SAMPLES2READ-1))/heasig.freq;
case 1
    M( : , 1)= (M( : , 1)- heasig.zerovalue(1));
    M( : , 2)= (M( : , 2)- heasig.zerovalue(1));
    M=M';
    M(1)=[];
    sM=size(M);
    sM=sM(2)+1;
    M(sM)=0;
    M=M';
    M=M/heasig.gain(1);
    TIME=(0:2*(SAMPLES2READ)-1)/heasig.freq;
%     clear sM
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;
M = [TIME' M];
M = M(offset+1:end,:); % skipping offset
% clear fid2 A M1H M2H PRR PRL TIME
% fprintf(1,'\\n$> LOADING DATA FINISHED \n');

%------ LOAD ATTRIBUTES DATA ----------------------------------------------
fid3=fopen([fname 'atr'],'r');
A= fread(fid3, [2, inf], 'uint8')';
status = 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))/heasig.freq;
ind= find(ATRTIME >= M(1,1) & ATRTIME <= M(end,1));
annotsig = [ATRTIME(ind),round(ANNOT(ind))];
% clear annoth hilfe sa saa fid3 status i ind ANNOT ATRTIME ATRTIMED

%------ DISPLAY DATA ------------------------------------------------------
% figure(1); clf, box on, hold on
% plot(M(:,1), M(:,2),'r');
% if heasig.nsig==2
%     plot(M(:,1), M(:,3),'b');
% end;
% for k=1:size(annotsig,1)
%     text(annotsig(k,1),0,num2str(annotsig(k,2)));
% end;
% xlim([M(1,1), M(end,1)]);
% xlabel('Time / s'); ylabel('Voltage / mV');
% string=['ECG signal ',heasig.name];
% title(string);
% fprintf(1,'\\n$> DISPLAYING DATA FINISHED \n');
% clear k string

% -------------------------------------------------------------------------

% fprintf(1,'\\n$> ALL FINISHED \n');