function [Q_on T_off QT] = get_QT(filename,channel,start,stop)
% new_X = remove_rep(Dif,X)
%
% Input Variables
%   filename - The record name
%   channel - The respective channel
%   start - The inital time for the reading
%   stop - The time for stopping
% Output Variables
%   Q_on- The Q onset
%   T_off- T offset
%   QT- The duration of the QT interval
%
% Description
%   Determines the QT interval
% Created
%   8/23/2006, Mohamed A Mneimneh, Marquette University
%
% This software is released under the terms of the GNU General
% Public License (http://www.gnu.org/copyleft/gpl.html).
% Modified


index = 1;
Q_on(1) = 0;
T_off(1) = 0;

cmd = ['rdsamp -r ' filename ];
cmd = [cmd sprintf(' -s %u -f %12.3f -t %12.3f', channel-1 , start, stop)];
[status, rrString] = system(cmd);
if status == 0
    signal(:,index) = sscanf(rrString, '%*d %d', inf);
else
    error('Couldn''t extract signal.');
end %if

t = start:1e-3:stop-1e-3;




close all

w = 0.03;
clear in_y y
[b,a] = butter(3,w);
y = filter(b,a,signal(:,index));

m = flipud(y);

y = flipud(filter(b,a,m));
y_filt = y;


[I status] = find_peak(y);

[s R_peaks Q_peaks S_peaks] = Calc_slopes(y,I,status);
flag_r =0;
wind_R = 600;
while length(R_peaks) < 3 && wind_R + 5000 < length(y)
    flag_r = 1;
    wind_R = wind_R + 400;
    [I status] = find_peak(y(wind_R:wind_R + 5000));
    [s R_peaks Q_peaks S_peaks] = Calc_slopes(y(wind_R:wind_R + 5000),I,status);
    [Q_peaks R_peaks S_peaks] = remove_rep (Q_peaks, R_peaks, S_peaks);
    I = wind_R + I;
    s = wind_R + s;
    R_peaks = R_peaks + wind_R;
    Q_peaks = Q_peaks + wind_R;
    S_peaks = S_peaks + wind_R;
end

if (flag_r)
    y = y(wind_R:wind_R + 5000);
end


R_sum = 0;
for i = 1:length(R_peaks)-1
    R_sum = R_sum + R_peaks(i+1) -R_peaks(i);
end

if length(R_peaks) >1
    R_avg(channel,index) = R_sum/(length(R_peaks)-1);
    wind = floor(0.7*R_sum/(length(R_peaks)-1));
else
    wind = 400;
end

if wind ==0
    Q_on(index) =0;
    T_off(index)=0;
    printf('unable to determine the QT interval');
end

in_y = movavg(y, wind,wind, 1);
in_y = flipud(movavg(flipud(in_y), wind,wind, 1));

t_org = t;
t = t_org(wind:end-wind);
y = -in_y(wind:end-wind)+ y(wind:end-wind);

if isempty(y)
    printf('unable to determine the QT interval');
end
[I status] = find_peak(y);
[s Rn_peaks Qn_peaks Sn_peaks] = Calc_slopes(y,I,status);

[Qn_peaks Rn_peaks Sn_peaks] = remove_rep (Qn_peaks, Rn_peaks, Sn_peaks);





[T_peaks P_peaks] = Find_peaks_TPn(y,Qn_peaks,Rn_peaks,Sn_peaks,I,status);



Q_on(channel,index) =0;
T_off(channel,index)=0;
if (length(T_peaks )>1) && ( (length(T_peaks ) == length(P_peaks)) || ((length(T_peaks ) == length(P_peaks) -1 ) ) )
    for i =1:length(T_peaks)-1
        if ( (T_peaks(i+1) - P_peaks(i)) < (Rn_peaks(i+1) - Rn_peaks(i)) )

            for j = T_peaks(i+1):P_peaks(i+1)
                S_T(j) =  integrate(t(T_peaks(i+1):j), y(T_peaks(i+1):j));
            end
            temp = intersect(S_T(T_peaks(i+1):P_peaks(i+1)),y(T_peaks(i+1):P_peaks(i+1)));
            T_off(channel,index) = wind + T_peaks(i+1) + temp(1);
            Q_lim = Qn_peaks(find (Qn_peaks > P_peaks(i)));
            for j = Q_lim(1):-1: P_peaks(i)
                S_Q(j) = integrate(t(Qn_peaks(i+1):-1:j), y(Qn_peaks(i+1):-1:j));
            end
            temp = intersect(S_Q(P_peaks(i) : Qn_peaks(i+1))/max(abs(S_Q)),y(P_peaks(i) : Qn_peaks(i+1))/max(abs(y(P_peaks(i) : P_peaks(i+1)))));
            Q_on(channel,index) = wind + P_peaks(i) + temp(end);

            break
        end
    end
end

QT(channel,index) = T_off(channel,index) - Q_on(channel,index);