%**************************************************************************
% James Lamberg 28Apr05
% EE 5821 Biomedical Modeling
% Imports a .dat file containing two signals in PhysioNet format 212
% Outputs a .MAT matrix
% The matrix can then be used in the Wavelets toolbox.
% This program is based on rddata.m by Robert Tratnig.
% The algorithm was based on a program written by Klaus Rheinberger.
% 10May05 - Updated so output file name was how I wanted it
% 17 May 2005 - Added check that input contains exactly two signals -- GBM
%**************************************************************************
clc; clear all; close all;
%**************************************************************************
% Get Data & User Inputs
%**************************************************************************
PATH = input('Path Where Data Is Stored, Else Leave Blank = ','s');
FILE = input('ECG File Name = ','s');       % Input String Filename
HEADERFILE = strcat(FILE,'.hea');           % Header In TXT Format
ATRFILE = strcat(FILE,'.atr');              % Attributes In Binary Format
DATAFILE = strcat(FILE,'.dat');             % ECG Data File
SAMPLESTART = input('ECG Start Time = ','s');
SAMPLEEND = input('ECG End Time = ','s');   
SAMPLESTART_1 = str2num(SAMPLESTART);       % The Start Time In Seconds
SAMPLEEND_1 = str2num(SAMPLEEND);           % The End Time In Seconds
%**************************************************************************
% 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
if nosig ~= 2
    error('Error: Input must have exactly 2 signals');
end;
sfreq = A(2);                               % Sample Rate Of Data
SAMPLESTART_2 = sfreq * SAMPLESTART_1;      % Samples/Second * Seconds
SAMPLEEND_2 = sfreq * SAMPLEEND_1;          % Samples/Second * Seconds
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 - 212 Only
    gain(k) = A(2);                         % Integers Per mV
    bitres(k) = A(3);                       % Bit Resolution
    zerovalue(k) = A(4);                    % Value Of ECG Zero Point
    firstvalue(k) = A(5);                   % First Value Of Signal
end;
fclose(fid1);
clear A;
%**************************************************************************
% Fix Minor Data Problems
%**************************************************************************
if SAMPLESTART_2 == 0                       % If User Starts At Zero
    SAMPLESTART_2 = 1;                      % We Start At The First Entry
end;
if SAMPLEEND_1 == 0                         % If User Wants To End At Zero
    error('Error: Results In No ECG Data');
end;
if SAMPLESTART_1 < 0 | SAMPLEEND_1 < 0      % If User Enters Negative Time
    error('Error: ECG Start & End Time Must Be Positive');
end;
if SAMPLESTART_1 > SAMPLEEND_1              % If User Sets End Before Start
    error('Error: ECG Start Time Should Be Before ECG End Time');
end;
%**************************************************************************
% Load Binary Data
%**************************************************************************
if dformat ~= [212,212], error('Error: File Not 212 Binary Format'); end;
signald = fullfile(PATH, DATAFILE);         % Data In 212 Format
fid2 = fopen(signald,'r');
A = fread(fid2, [3, SAMPLEEND_2], 'uint8')';
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('Error: 1st Bit Values Inconsistent'); end;
switch nosig
case 2
    M( : , 1) = (M( : , 1) - zerovalue(1))/gain(1);
    M( : , 2) = (M( : , 2) - zerovalue(2))/gain(2);
    TIME = (0:(SAMPLEEND_2 - 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*(SAMPLEEND_2)-1)/sfreq;
otherwise 
    disp('Error: Sorting Algorithm For > 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);
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
    elseif annoth == 61
    elseif annoth == 62
    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);
%**************************************************************************
% Manipulate Data So We Only Look At What The User Wants
%**************************************************************************
ECG_1_Temp = M(:,1);
ECG_1 = ECG_1_Temp(SAMPLESTART_2 : SAMPLEEND_2);
if nosig == 2
    ECG_2_Temp = M(:,2);
    ECG_2 = ECG_2_Temp(SAMPLESTART_2 : SAMPLEEND_2);
end;
Time_Adjusted = TIME(SAMPLESTART_2 : SAMPLEEND_2);
%**************************************************************************
% Display Data
%**************************************************************************
figure(1); clf, box on, hold on
plot(Time_Adjusted, ECG_1,'r');
if nosig == 2
    plot(Time_Adjusted, ECG_2,'b');
end;
for k = 1:length(ATRTIMED)
    text(ATRTIMED(k),0,num2str(ANNOTD(k)));
end;
xlim([Time_Adjusted(1), Time_Adjusted(end)]);
xlabel('Time (Seconds)'); ylabel('Voltage (mV)');
string = ['ECG Signal ',DATAFILE];
title(string);
fprintf(1,'\\n$> DISPLAYING DATA FINISHED \n');
%**************************************************************************
% Output Data File Into Current Working Directory
%**************************************************************************
save(strcat(FILE,'_ECG_',SAMPLESTART,'_',SAMPLEEND) ...
    , 'ECG_1' , 'ECG_2' , 'Time_Adjusted');
fprintf(1,'\\n$> ALL FINISHED \n');
%**************************************************************************
% End Of Code
%**************************************************************************