function twamag
% twamag.m - Detecting and Quantifying T-wave alternans 
% Copyright (C) 2008  Mahdi Zarrini  <mahdi.zrn63@gmail.com>
% Physionet Challenge 2008 

clc
clear all
close all
warning off all
PATH = input('Enter the record path:','s') ;
twann = input('Enter the number of record(0 to 99):','s') ;

if length(twann) == 1
    twann = strcat('0',twann) ;
end

FILE = strcat('twa',twann) ;
HEADERFILE = strcat(FILE,'.hea');           % Header In TXT Format
QRSFILE = strcat(FILE,'.qrs');              % QRS In Binary Format
DATAFILE = strcat(FILE,'.dat');             % ECG Data File

%**************************************************************************
% Load Header Data
%**************************************************************************
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]);
    gain(k) = A(2);                         % Integers Per mV
end;
fclose(fid1);
clear A;

%**************************************************************************
% Load Binary Data
%**************************************************************************
sigdata = fullfile(PATH,DATAFILE) ;
fid = fopen(sigdata);
A = fread(fid,'uint8');
fclose(fid);
n = nosig ;
size_a = size(A,1);

B1 = zeros(size_a/(2 * n) , 1);

for i=1:size_a/(2 * n)
    B1(i) = A(i*(2 * n) - 1) + A(i*(2 * n) ) * 256 ;
    
    if (B1(i) >= 32768)
        B1(i) = B1(i) - 65536;
    end
end

B1 = B1/gain(1) ;

%**************************************************************************
SAMPLEEND = length(B1) ;
TIME = (0:(SAMPLEEND)-1)/sfreq;
% figure; plot(TIME,B1) ;         %Plot ECG in mV/Sec
%**************************************************************************
% Load Binary Data
%**************************************************************************
% Eliminate Baseline Drift
s1 = B1 ;
% clear B1 
s2 = smooth(s1,150) ;
ecgsmooth = s1 - s2 ;
% Apply WT
[C,L] = wavedec(ecgsmooth,8,'db4');
[d1,d2,d3,d4,d5,d6,d7,d8] = detcoef(C,L,[1,2,3,4,5,6,7,8]) ;

% Denoise
[thr,sorh,keepapp] = ddencmp('den','wv',ecgsmooth) ;
ecg = wdencmp('gbl',C,L,'db4',8,thr,sorh,keepapp) ;
%**************************************************************************
% Load Attributes Data
%**************************************************************************
atrd = fullfile(PATH, QRSFILE);
fid3 = fopen(atrd,'r');
A = fread(fid3, [2, inf], 'uint8')';
fclose(fid3);
ATRTIME = [];
ANNOT = [];
sa = size(A);
saa = sa(1);
i = 1;
success = 1 ;
while i <= saa
    annoth = bitshift(A(i,2),-2);
    if annoth == 59
        success = 0 ;
    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;
if success == 1
    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);
    Rpeak = floor(ATRTIMED*sfreq) ;
    Beats = length(Rpeak) ;
end
%**************************************************************************
% R-detection if reading from Attributes failed
%**************************************************************************
if success == 0
    pos = ecg>0 ;
    pos = ecg.*pos ;
    R_detect_squared=pos.^2 ; 

    %  Thresholding 
    temp = R_detect_squared(1:500) ;
    max_value = max(temp) ;
    mean_value = mean(R_detect_squared) ;
    thr = (max_value - mean_value)/2 ;

    Beats = 0;
    b = 1 ;
    while b<length(R_detect_squared)
    if (R_detect_squared(b)<thr)&&(R_detect_squared(b+1)>thr)
        Beats = Beats + 1 ;
        c = b+1;
        while R_detect_squared(c)>thr
            c = c + 1 ;
        end
        t = ecg(b:c);
        [m Rpeak(Beats)] = max(t) ;
        Rpeak(Beats) = Rpeak(Beats) + b - 1 ;
        b = b + 100 ;
    else
        b = b + 1 ;
    end
    end
end

%**************************************************************************
% Separate The Beats And Detect T-peak
%**************************************************************************
Tp = zeros(Beats,1) ;

for i=1:Beats
    if i~=1
        RR = Rpeak(i) - Rpeak(i-1);
        LB = Rpeak(i) - 0.4*RR ; 
    else
        LB = 1 ;
    end
    if i~=Beats
        RR = Rpeak(i+1) - Rpeak(i);
        UB = 0.6*RR + Rpeak(i) ;
    else
        UB = length(ecg) ;
    end
    OneBeat = ecg((LB+20):UB); 

        Tp(i)=tdetect(OneBeat) ;
        Tp(i) = floor(Tp(i)+ LB+20) ;
        
        % Cancel Floor error 
        a = ecg(Tp(i)) ;
        b = ecg(Tp(i)+1) ;
        if Tp(i)~= 1
            c = ecg(Tp(i)-1) ;
            if b > a
                Tp(i) = Tp(i)+1 ;
            elseif c > a
                Tp(i) = Tp(i)-1 ;
            end
        end
end

figure; plot(B1) ;

for i=1:Beats
        text(Tp(i),B1(Tp(i)),'*','Color','r') ;
end
p = max(ecg(Tp)) - min(ecg(Tp))


%% Twave Detection
function Tp = tdetect(onebeat)

Crit = .008*(max(onebeat)-min(onebeat)) ; %Condition
[peak R]=max(onebeat) ;               %Rpeak detection
Le1 = length(onebeat) ;
Tp   = 0 ;

% Linearity & Slope Condition
for i=R:Le1
    f=0 ;
    a=onebeat(i) ;
    for j=0:19
        if (i+j)<=length(onebeat)
            if (a-onebeat(i+j))>Crit
                f=1;
            end
        end
    end
    if f==0
        break;
    end   
end

Sis = i + 19;

% Finding Tpeak
if Sis<(Le1-16)
    for i=Sis:Le1-16
        w1 = onebeat(i-16)-onebeat(i);
        w2 = onebeat(i)-onebeat(i+16) ;
        wings(i) = w1*w2 ;
    end
    [minwings Tp] = min(wings);
    Tp = Tp;
else
    return 
end