%% ************************************************************************

%     twa.m - Detecting and Quantifying T-Wave Alternans (TWA)
%
%     Copyright (C) 2008  Jubair Sieed <jubir@ymail.com>
%     Department of Electrical and Electronic Engineering,
%     Bangladesh University of Engineering and Technology.
%
%     This software is released under the terms of the GNU General
%     Public License (http://www.gnu.org/copyleft/gpl.html).
%
%     Interpretation: ta = 0 ==> No TWA detected
%                  ta = **** ==> TWA detected of maximum amplitude **** mV
%
%     Any error message implies that the data set is not appropriate for
%     this algorithm or mostly absence of TWA in it.
   
%% ************************************************************************     

clc, clear all, close all;

%% ************************************************************************
% Load matfile named twa**.mat in the working directory which was imported
% from corresponding ecg data file in text format. Original database is
% available at: http://physionet.org/pn3/twadb/
% *************************************************************************

load twa98; 

[p q] = size(data); 

samplingrate = 500; 

if q<5
    ecg = data(:,2)'; 
    threshold = 0.55*max(ecg);
else
    ecg = data(:,10)';
    threshold = 0.85*max(ecg);
end

l = length(ecg);

%% ************************************************************************
% Filtering (moderate filtering is used as it may vary the alternan values)
%**************************************************************************

fresult = fft(ecg);
fresult(1 : round(length(fresult)*1/samplingrate)) = 0;
fresult(end - round(length(fresult)*1/samplingrate) : end) = 0;
filtered = real(ifft(fresult));

figure(1)
subplot(211), 
plot(ecg), axis([1 5e3 min(ecg)-0.1 max(ecg)+0.1])
subplot(212), plot(filtered), axis([1 5e3 min(ecg)-0.1 max(ecg)+0.1])

%% ************************************************************************
% Find Local peaks (R/T) and position
%**************************************************************************

ecg = filtered;
peak = zeros(1,l);

threshold = 0.8*max(ecg);

for i = 2:l-1
    if ecg(i)>threshold & ecg(i)>ecg(i+1) & ecg(i)>ecg(i-1)
       peak(i) = ecg(i);
       i = i+10;
    end
end

pos = find(peak);
HR = length(pos);
r = ecg(pos);

%% ************************************************************************
% Find T-peak and Amplitude
%  ************************************************************************

m = 2*floor(HR/2); 
t = zeros(1,m);

for j=1:m-2
    if pos(j+1)-pos(j)>100
        t(j) = max(ecg(pos(j)+20:pos(j+1)-20));
    end
end

if q>4 & r(10)<1.75*t(10)
    t = r ;
end

figure(2)
subplot(211), stem(r)
subplot(212), stem(t)

todd = t(1:2:end-1); 
teven = t(2:2:end); 

twave = todd-teven;

figure(3)
stem(twave)

z=0;
for k=1:length(twave)-1
    if twave(k)/(twave(k+1)+eps)<0
        z=z+1;
    end
end

if z<0.35*length(twave) & HR>80
    ta = max(abs(twave(2:end-1)))
else
    ta=0
end

%% ************************************************************************