% QT_gradient_detection.m - Find QT intervals based on the gradient rules
% Copyright (C) 2006 Vaclav Chudacek
% This software is released under the terms of the GNU General
% Public License (http://www.gnu.org/copyleft/gpl.html).
function QT_detection_using_Gradient_Rules_CinC_2006(path,plot_enabled)
if (nargin~=2),
disp(sprintf('This function requires two arguments. Setting to default'));
path = 'C:\!Vasek\Work\!CINC 2006\!Data\ptbdb';
plot_enabled = 0;
end
recordFile = 'RECORDS'; % Record names
load 'CINC_2006_Reference_results.txt'; % Official results for comparison as of 07.09.2006
CINC_2006_results = zeros(549,2); % Our results
CINC_2006_QT_lengths = zeros(1,549);
CINC_2006_differences = zeros(1,549); % Difference of our results to the official
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
recordFileFull = fullfile(path, recordFile);
fid_record = fopen(recordFileFull,'r');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% begin_at = 119;
% for i = 1:begin_at,
% fileName = fgetl(fid_record);
% end
for counter = 1:549,
fileName = fgetl(fid_record);
headerFile = fullfile(path,[fileName '.hea']);
% vcgFile = fullfile(path,[fileName '.xyz']);
ecgFile = fullfile(path,[fileName '.dat']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------ LOAD HEADER -------------------------------------------
fprintf(strcat('Loading new record #',num2str(counter),': ',fileName));
fprintf('\n');
fid_header = fopen(headerFile,'r');
z = fgetl(fid_header);
fclose(fid_header);
%------ LOAD DATA --------------------------------------------
fid_ecg = fopen(ecgFile,'r');
ECG = fread(fid_ecg, [12, inf], 'int16');
fclose(fid_ecg);
% fid_vcg = fopen(vcgFile,'r');
% VCG = fread(fid_vcg, [3, inf], 'int16');
% fclose(fid_vcg);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Very, very basic filtration - just to get rid of the worst...
for i = 1:12
herz = 1000/length(ECG(i,:));
A = fft(ECG(i,:));
A(1:round(0.66/herz)) = 0;
A(end-round(0.66/herz):end)=0;
A((round(50/herz)-5):(round(50/herz)+5)) = 0;
A(end-(round(50/herz)+5):end-(round(50/herz)-5)) = 0;
ECGf(i,:) = (real(ifft(A)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% R-peak detection based on PT_Hamilton algorithm %%%%%%%
% Source code distributed under GNU license can be found on
% http://www.eplimited.com/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% START %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nbr_heart_beats, r_peak_sample_list] = PanTomkinsHamilton(ECGf(5,:),1000);
[nbr_heart_beatsI, r_peak_sample_listI] = PanTomkinsHamilton(ECGf(1,:),1000);
[nbr_heart_beatsV4, r_peak_sample_listV4] = PanTomkinsHamilton(ECGf(10,:),1000);
r_peak_sample_list = r_peak_sample_list*5;
posR = r_peak_sample_list;
posRI = r_peak_sample_listI;
posRV4 = r_peak_sample_listV4;
RR_dist = zeros(1,length(min(100,length(posR)-1)));
selected_beat = 5; next_beat = 6;
for i = 1:min(100,length(posR)-1),
RR_dist(i) = (posR(i+1) - posR(i));
end
%%% Correction of detected peaks - R-T mess, wrong detection
if ((RR_dist(selected_beat))>(1.2*RR_dist(selected_beat-1))||(RR_dist(selected_beat))<(0.8*RR_dist(selected_beat-1)))||...
((RR_dist(selected_beat))>(1.2*median(RR_dist))||(RR_dist(selected_beat))), %%% added because of file #543
if ((RR_dist(selected_beat))>=1500), %%% added because of files #543, #315, #415
selected_beat = selected_beat - 1; next_beat = next_beat +1;
elseif (RR_dist(selected_beat+2))<=(1.2*RR_dist(selected_beat))&&...
(RR_dist(selected_beat+2))>=(0.8*RR_dist(selected_beat))&&...
(RR_dist(selected_beat))<425,
if RR_dist(selected_beat)<320,
if RR_dist(selected_beat)<500
next_beat = next_beat+1;
else
selected_beat = selected_beat+1; next_beat = next_beat+2;
end
else
selected_beat = selected_beat-1;
end
elseif (RR_dist(selected_beat+1))<=(1.2*RR_dist(selected_beat))&&(RR_dist(selected_beat+1))>=(0.8*RR_dist(selected_beat))
selected_beat = selected_beat+1; next_beat = next_beat +1;
end
end
posRpeak_out = posR(selected_beat);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% R-peak detection by PanTomkins&Hamilton algortihm END %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Q-on computation START %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q_window_span_left = 150; Q_window_span_right = 50;
R_window_span_left = 70; R_window_span_right = 10;
T_window_span_right_1 = 200; T_window_span_right_2 = 200;
Q_comps = 0; mistakeQ = 0; T_comps = 0; mistakeT = 0;
Q_on_abs = zeros(1,12);
for i = 1:12,
%%% Correction of R-peak for each lead
[temp correctionR_temp] = max((ECGf(i,posR(selected_beat)-R_window_span_left:posR(selected_beat)+R_window_span_right)));
correctionR = correctionR_temp - R_window_span_left;
posR_adjusted = posR(selected_beat) + correctionR;
%%% 1st step of Q-on computation - see function: FIND_Q_ON
try
Q_on_temp = find_Q_on(ECGf(i,posR_adjusted - Q_window_span_left:posR_adjusted + Q_window_span_right));
catch
disp(strcat(['!!!!! MISTAKE of: Q-on computation on lead:', (num2str(i)),'!!!!!']));
mistakeQ = 1;
Q_on_temp = 0;
end
if mistakeQ == 0,
Q_comps = Q_comps+1;
Q_on_rel(Q_comps) = Q_on_temp + correctionR + 1;
Q_on_abs(i) = posR(selected_beat) - Q_on_rel(Q_comps);
else
mistakeQ = 0;
end
end
% All possible Qons are sorted and the smallest are deleted
if min(Q_on_rel)<10,
Q_on_rel_tmp = sort(Q_on_rel,'ascend');
clear('Q_on_rel');
Q_on_rel = Q_on_rel_tmp(max(find(Q_on_rel_tmp<10))+1:end);
else
Q_on_rel = sort(Q_on_rel,'ascend');
end
% 2nd step of Q-on computations - see function PRESCISE_Q
for i = 2:2
Q_on_final_temp(i,:) = precise_Q_on(ECGf(i,posR(selected_beat) - 150 : posR(selected_beat) + Q_on_rel(end)+25 -150),Q_on_rel);
end
% Maximal score is computed to select the best fitting Q-on
maxScoreQ = max(sum(Q_on_final_temp(:,:)));
if maxScoreQ == 0,
Q_on_rel_final = median(Q_on_rel);
else
maxScoreSelected = find(sum(Q_on_final_temp(:,:))>=maxScoreQ);
Q_on_rel_final = Q_on_rel(maxScoreSelected(1));
end
if (median(Q_on_rel) - Q_on_rel_final) < -25,
Q_on_rel_final = round(median(Q_on_rel));
end
Q_on_rel_final = Q_on_rel_final - 150;
Q_on_abs_final = Q_on_rel_final + posR(selected_beat);
% 3rd step of Q-on computations - see function: CORRECT_Q
correctionQ = [];
leads = [1 2 8 9];
for i = 1:length(leads);
correctionQ = [correctionQ correct_Q_on(ECGf(leads(i),:),Q_on_abs_final,Q_on_abs_final-20,Q_on_abs_final+30)];
end
correctionQ = median(correctionQ);
Q_on_rel_final = Q_on_rel_final + correctionQ;
Q_on_abs_final = Q_on_abs_final + correctionQ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Q-on computation END %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R_pos_in_summedECG = 75; nextR_pos_summedECG_distance = 200;
summedECG = zeros(1,751);
for i = 1:12,
% Computation of summed and filtered ECG
if (i == 8)||(i == 9),
if(i == 8),
summedECG = slidingavg(slidingavg(slidingavg(ECGf(12,(posR(selected_beat)-R_pos_in_summedECG):(posR(next_beat)-nextR_pos_summedECG_distance)),50),10),3);
else
summedECG = summedECG + slidingavg(slidingavg(slidingavg(ECGf(12,(posR(selected_beat)-R_pos_in_summedECG):(posR(next_beat)-nextR_pos_summedECG_distance)),50),10),3);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Inverting signal when negative
if abs(min(summedECG))>max(summedECG),
summedECG = -1*summedECG;
end
%Correction of R-peak on summedECG in +-20 sample radius
[unused correctionR] = max(summedECG(R_pos_in_summedECG-20:R_pos_in_summedECG+20));
R_pos_in_summedECG = R_pos_in_summedECG - (20-correctionR+1);
Q_on_for_summedECG = R_pos_in_summedECG - Q_on_rel_final;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% T-off computation START %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_error = 0;
%Searching for first and second maxima/minima of the T-wave
[unused maxT] = max(abs(summedECG(R_pos_in_summedECG+min(150,round(RR_dist(selected_beat)/2)):(min(length(summedECG(R_pos_in_summedECG:end)),400+R_pos_in_summedECG))))-summedECG(Q_on_for_summedECG));
if ~exist('maxT','var'),
maxT = 0;
end
if summedECG(R_pos_in_summedECG+200+maxT)<0,
zaporne_T = 1;
maxT = maxT + R_pos_in_summedECG + 200;
[velikostT2 maxT2] = max(summedECG(maxT:end - min(length(summedECG)-maxT,100))-summedECG(maxT));
maxT2 = maxT2 -1;
if ((velikostT2<(summedECG(end - min(length(summedECG)-maxT,100))-summedECG(maxT))+100)&&(velikostT2>(summedECG(end - min(length(summedECG)-maxT,100))-summedECG(maxT))-100)),
maxT2=0;
end
else
maxT2 = 0;
maxT = maxT + R_pos_in_summedECG + 200;
zaporne_T =0;
end
if maxT>length(summedECG)-50,
maxT = length(summedECG)-200;
end
maxT2_rel = maxT2 - maxT;
%Signal used for T-off detection is shifted by value of Q-on
Tecg = summedECG(maxT:end-50)-summedECG(Q_on_for_summedECG);
%%% 1st step of T-off computation - see function: FIND_T_OFF
[Toff_temp, T_error] = find_T_off(Tecg,maxT, maxT2, maxT2_rel,zaporne_T);
% posTmax_out = maxT - R_pos_in_summedECG + posR(selected_beat);
T_off_rel_final = Toff_temp + maxT;
T_off_abs_final = Toff_temp + maxT - R_pos_in_summedECG + posR(selected_beat);
% 2nd step of T-off computations - see function: CORRECT_T_OFF
correctionT = [];
leads = [1 2 8 9 10];
for i = 1:length(leads);
correctionT(i) = correct_T_off(ECGf(leads(i),T_off_abs_final-200:T_off_abs_final));
end
correctionTfinal = min(correctionT);
T_off_abs_final = T_off_abs_final - correctionTfinal;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% T-off computation START %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting of the results START
if plot_enabled,
%%% Draws detected R(selected_beat):R(next_beat)
% figure(counter*100);plot(ECGf(1,posR(selected_beat):posR(next_beat)));
% Draws preliminary and corrected results
figure(counter);
plot(ECGf(2,Q_on_abs_final -100:T_off_abs_final+100),'r');
hold on;
plot(ECGf(1,Q_on_abs_final -100:T_off_abs_final+100));
plot(ECGf(8,Q_on_abs_final -100:T_off_abs_final+100));
plot(ECGf(9,Q_on_abs_final -100:T_off_abs_final+100));
plot(100,-500:500,'k');
plot(100 + correctionQ, -500:500,'g');
%%%%%%%%%%%
plot(length(ECGf(i,Q_on_abs_final -100:T_off_abs_final+100))-100+correctionTfinal,-500:500,'k');
plot(length(ECGf(i,Q_on_abs_final -100:T_off_abs_final+100))-100 - correctionTfinal,-500:500,'g');
end
% Saving final results
CINC_2006_results(counter,:) = [Q_on_abs_final,T_off_abs_final];
CINC_2006_QT_lengths(counter) = (T_off_abs_final-Q_on_abs_final);
CINC_2006_differences(counter) = abs(CINC_2006_Reference_results(counter)-CINC_2006_QT_lengths(counter));
clear ECG*; clear Q*; clear R*; clear T*; clear cor*; clear m*; clear n*; clear pos*; clear r*; clear z*; clear s*; clear t*; clear u*;
end
fclose(fid_record);
% Final score computation based on final reference QT measurements (not included)
omitted_results = find((CINC_2006_QT_lengths > 500) | (CINC_2006_QT_lengths < 250));
if length(omitted_results)>27
disp('Sorry you are omitting to many results. More than 27.')
end
square = (CINC_2006_differences.^2);
square(omitted_results) = -1;
RESULT = sqrt(mean(square(find(square ~= -1))));
RESULT = RESULT*((549-length(omitted_results))/549)
save('QTresults','CINC_2006_results')