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

%     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 amplitude **** micro Volts
%
%     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 twa50; 

[p q] = size(data); 

samplingrate = 500; 

if q<5
    ecg = data(:,2)'; 
else
    ecg = data(:,10)';
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(data(:,1)',ecg), axis([0 10 min(ecg)-0.1 max(ecg)+0.1])
xlabel('Time (seconds)')
ylabel('ECG amplitude (milli Volts)')
title('Plot of given ECG signal')

subplot(212),plot(data(:,1)',filtered),axis([0 10 min(ecg)-.1 max(ecg)+.1])
xlabel('Time (seconds)')
ylabel('ECG amplitude (milli Volts)')
title('Plot of filtered ECG signal')

%% ************************************************************************
% 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)
    s = t ;
    t = r ;
    r = s ;
end

figure(2)
subplot(211), stem(r)
xlabel('sequence no.')
ylabel('R-peak amplitude (milli Volts)')
title('Plot of detected R-peaks')

subplot(212), stem(t)
xlabel('sequence no.')
ylabel('T-wave peak amplitude (milli Volts)')
title('Plot of detected T-peaks')

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

twave = todd-teven;

figure(3)
stem(twave)
xlabel('sequence no.')
ylabel('T-wave bit ti bit differences (milli Volts)')
title('Plot of TWA measurement')

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

if z<0.4*length(twave) & HR>72 & HR<180 & var(twave)<0.05
    ta = 1000*max(abs(twave(3:end-2)))
else
    ta=0
end

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