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