function [R, QRSon, QRSend, T, Ton, Tend, P, Pon, Pend] = detectW(signal)
% 
%   [R, QRSon, QRSend, T, Ton, Tend, P, Pon, Pend] = detectW(signal)
% 
%   This function detects the apperance of R, P and T waves in ECG
%   recordings sampled at 250Hz. It also delineates the ECG, marking the
%   beginning and ending of each individual wave. To do so, it uses  
%   wavelet analysis.
% 
%   Input:
% 
%   'signal' is a column vector of doubles containing the signal that will be
%   processed. The signal must have a sampling frequency of 250Hz
% 
%   Outputs:
% 
%   'R' is a vector containing the R wave detection results. The vector
%   contains the indexes of the samples where detections were made.
%
%   'QRSon' is a vector containing the indexes of the samples where the QRS
%   complex beginnings where detected
%
%   'QRSend' is a vector containing the indexes of the samples where the QRS
%   complex endings where detected
%
%   'T' is a vector containing the T wave detection results. The vector
%   contains the indexes of the samples where detections were made.
%
%   'Ton' is a vector containing the indexes of the samples where the T
%   waves beginnings where detected
%
%   'Tend' is a vector containing the indexes of the samples where the T
%   waves endings where detected
%
%   'P' is a vector containing the P wave detection results. The vector
%   contains the indexes of the samples where detections were made.
%
%   'Pon' is a vector containing the indexes of the samples where the P
%   waves beginnings where detected
%
%   'Pend' is a vector containing the indexes of the samples where the T
%   waves endings where detected
%
%   This routine was developed by the Simon Bolivar University's Applied
%   Biophysics and Bioengineering Group (GBBA-USB) for free and public use. 
%   It is protected under the terms of the Creative Commons Attribution-Non
%   Comercial 4.0 International License, no profit may come from any 
%   applicaton that uses this function and proper credit must be given to 
%   the GBBA-USB. For further information consult 
%   http://creativecommons.org/licenses/by-nc/4.0/
%
%   Author: Carlos Alberto Ledezma

fs = 250;

w = 0 : 2*pi/length(signal) : 2*pi - 1/length(signal);

%Construction of highpass filters
G1(:,1) = 4 * 1j * exp(1j .* w ./ 2) .* sin(w./2);
G2(:,1) = 4 * 1j * exp(1j * 2.* w ./ 2) .* sin(2*w./2);
G3(:,1) = 4 * 1j * exp(1j *4.* w ./ 2) .* sin(4*w./2);
G4(:,1) = 4 * 1j * exp(1j *8.* w ./ 2) .* sin(8*w./2);

%Construction of lowpass filters
H1(:,1) = exp(1j .* w ./ 2).*(cos(w./2)).^3;
H2(:,1) = exp(1j *2 .* w ./ 2).*(cos(2*w./2)).^3;
H3(:,1) = exp(1j *4 .* w ./ 2).*(cos(4*w./2)).^3;


%Signal decomposition

S1 = fft(signal);
scale1 = ifft(S1 .* G1);

S2 = S1 .* H1;
scale2 = ifft(S2 .* G2);

S3 = S2 .* H2;
scale3 = ifft(S3 .* G3);

S4 = S3 .* H3;
scale4 = ifft(S4 .* G4);


%Definition of detection variables

%QRS Delineation constants
eps = [rms(scale1), rms(scale2), rms(scale3), 0.5*rms(scale4)];
epsT = 0.25*rms(scale4);
epsP = 0.02*rms(scale4);
annT = 1; %Counter for T annotations
annP = 1; %Counter for P annotations
%Annotation positions
n1 = 0;
n1f = 0;
n2 = 0;
n2f = 0;
n3 = 0;
n3f = 0;
n4 = 0;
n4f = 0;
%Scale 4 refractory period
win4 = (200/1000)*fs;


%Detection results
R = 0; %R peaks
QRSon = 0;%QRS beginning
QRSend = 0;%QRS end
T = 0;%T peak
Ton = 0;%T beginning
Tend = 0;%T end
P = 0;%P peak
Pon = 0;%P beginning
Pend = 0;%P end

%QRS detection
%Neighborhood
neigh = floor((50/1000)*fs);
%Detection of modulus maxima in scale 4
k = 1;
for i = 1 : length(scale4)
    if scale4(i) > eps(4)
        if length(n4) == 1
            [~,ind] = max(scale4(i:i+neigh));
            n4(k) = i + ind;
            k = k+1;
        else
            if i - n4(k-1) > win4;
                if i + neigh <= length(scale4)
                    [~,ind] = max(scale4(i:i+neigh));
                    n4(k) = i + ind;
                    k = k+1;
                else[~,ind] = max(scale4(i:length(scale4)));
                    n4(k) = i + ind;
                    k = k+1;
                end
                    
            end
        end
    end
end
%Searching scale 3 for modulus maxima
for i = 1 : length(n4)
    if n4(i)-neigh > 0
        if n4(i) + neigh <= length(scale3)
            if ~isempty(find(scale3(n4(i)-neigh:n4(i)+neigh) > eps(3),1))
                [~,ind] = max(abs(scale3(n4(i)-neigh:n4(i)+neigh)));
                n3(i) = n4(i) - neigh + ind;
            else
                n3(i) = 0;
            end
        else
            if ~isempty(find(scale3(n4(i)-neigh:length(scale3)) > eps(3),1))
                [~,ind] = max(abs(scale3(n4(i)-neigh:length(scale3))));
                n3(i) = n4(i) - neigh + ind;
            else
                n3(i) = 0;
            end
        end
            
    else
        n3(i) = 0;
    end
end

%Searching scale 2 for modulus maxima
for i = 1 : length(n4)
    if n4(i)-neigh > 0
        if n4(i) + neigh <= length(scale2)
            if ~isempty(find(scale2(n4(i)-neigh:n4(i)+neigh) > eps(2),1))
                [~,ind] = max(abs(scale2(n4(i)-neigh:n4(i)+neigh)));
                n2(i) = n4(i) - neigh + ind;
            else
                n2(i) = 0;
            end
        else
            if ~isempty(find(scale2(n4(i)-neigh:length(scale2)) > eps(2),1))
                [~,ind] = max(abs(scale2(n4(i)-neigh:length(scale2))));
                n2(i) = n4(i) - neigh + ind;
            else
                n2(i) = 0;
            end
        end
            
    else
        n2(i) = 0;
    end
end

%Searching scale 1 for modulus maxima
for i = 1 : length(n4)
    if n4(i)-neigh > 0
        if n4(i) + neigh <= length(scale1)
            if ~isempty(find(scale1(n4(i)-neigh:n4(i)+neigh) > eps(1),1))
                [~,ind] = max(abs(scale1(n4(i)-neigh:n4(i)+neigh)));
                n1(i) = n4(i) - neigh + ind;
            else
                n1(i) = 0;
            end
        else
            if ~isempty(find(scale1(n4(i)-neigh:length(scale1)) > eps(1),1))
                [~,ind] = max(abs(scale1(n4(i)-neigh:length(scale1))));
                n1(i) = n4(i) - neigh + ind;
            else
                n1(i) = 0;
            end
        end
    else
        n1(i) = 0;
    end
end

%Error correction
k = 1;
for i = 1 : length(n4)
    if n3(i) ~= 0 && n2(i) ~= 0 && n1(i) ~= 0
        n4f(k) = n4(i);
        n3f(k) = n3(i);
        n2f(k) = n2(i);
        n1f(k) = n1(i);
        k = k+1;
    end
end

%R-peak detecton in scale 1
%Neighborhood
neigh = floor((100/1000)*fs);
k=1;
for i = 1 : length(n1f)
    if scale1(n1f(i)) > 0
        if n1f(i) - neigh > 0
            if n1f(i) + neigh <= length(scale1)
                temp = find(scale1(n1f(i)-neigh:n1f(i)+neigh) < 0);
                temp = n1f(i) - neigh + temp;
            else
                temp = find(scale1(n1f(i)-neigh:length(scale1)) < 0);
                temp = n1f(i) - neigh + temp;
            end
        else
            temp = find(scale1(1:n1f(i)+neigh) < 0);
        end
                
        [~,ind] = max(abs(scale1(temp)));
        ind = temp(ind);
        if ind > n1f(i)
            [~,ann] = min(abs(scale1(n1f(i):ind)));
            R(k) = n1f(i) + ann;
            k=k+1;
        else
            [~,ann] = min(abs(scale1(ind:n1f(i))));
            R(k) = ind + ann;
            k=k+1;
        end
    else
        if n1f(i) - neigh > 0
            if n1f(i) + neigh <= length(scale1)
                temp = find(scale1(n1f(i)-neigh:n1f(i)+neigh) > 0);
                temp = n1f(i) - neigh + temp;
            else
                temp = find(scale1(n1f(i)-neigh:length(scale1)) > 0);
                temp = n1f(i) - neigh + temp;
            end
        else
            temp = find(scale1(1:n1f(i)+neigh) > 0);
        end

            
        [~,ind] = max(abs(scale1(temp)));
        ind = temp(ind);
        if ind > n1f(i)
            [~,ann] = min(abs(scale1(n1f(i):ind)));
            R(k) = n1f(i) + ann;
            k=k+1;
        else
            [~,ann] = min(abs(scale1(ind:n1f(i))));
            R(k) = ind + ann;
            k=k+1;
        end
    end
end

%Beat Delineation


for i = 1 : length(R)
    
    neigh = floor((70/1000)*fs);
    closeneigh = floor((5/1000)*fs);
    %QRS beginning
    if R(i) - neigh > 0
        if R(i) + neigh <= length(scale2)
            QRSpre = 0.06*max(abs(scale2(R(i)-neigh:R(i)+neigh)));
            nfirst = R(i) - neigh - 1 + find(abs(scale2(R(i) - neigh:R(i))) > QRSpre, 1, 'first');
        else
            QRSpre = 0.06*max(abs(scale2(R(i)-neigh:length(scale2))));
            nfirst = R(i) - neigh - 1 + find(abs(scale2(R(i) - neigh:R(i))) > QRSpre, 1, 'first');
        end
    else
        QRSpre = 0.06*max(abs(1:R(i)+neigh));
        nfirst = find(abs(1:R(i)) > QRSpre, 1, 'first');
    end
        
    if isempty(nfirst)
        nfirst = R(i) - neigh;
    end
    [~,ind] = max(abs(scale2(nfirst:nfirst+closeneigh)));
    nfirst = nfirst - 1 + ind;
    
    if scale2(nfirst) < 0
        xQRSon = 0.07*scale2(nfirst);
    else
        xQRSon = 0.05*scale2(nfirst);
    end
    
    if nfirst - neigh > 0
        if ~isempty(find(abs(scale2(nfirst - neigh : nfirst)) < xQRSon, 1, 'last'))
            QRSon(i) = nfirst - neigh - 1 + find(abs(scale2(nfirst - neigh : nfirst)) < xQRSon, 1, 'last');
        else
            QRSon(i) = nfirst;
        end
    else
        if ~isempty(find(abs(scale2(1 : nfirst)) < xQRSon, 1, 'last'))
            QRSon(i) = find(abs(scale2(1 : nfirst)) < xQRSon, 1, 'last');
        else
            QRSon(i) = nfirst;
        end
    end
        
    %QRS ending
    if R(i) - neigh > 0
        if R(i) + neigh <= length(scale2)
            QRSpost = 0.09*max(abs(scale2(R(i)-neigh:R(i)+neigh)));
            nlast = R(i) - 1 + find(abs(scale2(R(i):R(i)+neigh)) > QRSpost, 1, 'last');
        else
            QRSpost = 0.09*max(abs(scale2(R(i)-neigh:length(scale2))));
            nlast = R(i) - 1 + find(abs(scale2(R(i):length(scale2))) > QRSpost, 1, 'last');
        end
    else
        QRSpost = 0.09*max(abs(1:R(i)+neigh));
        nlast = R(i) - 1 + find(abs(scale2(R(i):R(i)+neigh)) > QRSpost, 1, 'last');
    end
    if isempty(nlast)
        nlast = R(i) + neigh;
    end
    [~,ind] = max(abs(scale2(nlast-closeneigh:nlast)));
    nlast = nlast - closeneigh - 1 + ind;
    
    if scale2(nlast) < 0
        xQRSend = 0.71*scale2(nlast);
    else
        xQRSend = 0.125*scale2(nlast);
    end
    
    if nlast + neigh <= length(scale2);
        if ~isempty(find(abs(scale2(nlast : nlast + neigh)) < xQRSend, 1, 'first'))
            QRSend(i) = nlast - 1 + find(abs(scale2(nlast : nlast + neigh)) < xQRSend, 1, 'first');
        else
            QRSend(i) = nlast;
        end
    else
        if ~isempty(find(abs(scale2(nlast : length(scale2))) < xQRSend, 1, 'first'))
            QRSend(i) = nlast - 1 + find(abs(scale2(nlast : length(scale2))) < xQRSend, 1, 'first');
        else
            QRSend(i) = nlast;
        end
    end
    
    %T detection and delineation
    
    neigh = floor((200/1000)*fs);
    win = floor((15/1000)*fs);
    if QRSend(i) + neigh < length(scale4)
        %Finding if T wave is present
        pos = QRSend(i)+win - 1 + find(scale4(QRSend(i)+win:QRSend(i)+win+neigh) > epsT);
        neg = QRSend(i)+win - 1 + find(scale4(QRSend(i)+win:QRSend(i)+win+neigh) < -1*epsT);
        if ~isempty(pos) && ~isempty(neg)
            gT = 0.125*max(abs(scale4(QRSend(i):QRSend(i)+neigh)));
            if (max(scale4(pos)) > gT) && (min(scale4(neg))<-1*gT)
                [~,indPos] = max(scale4(pos));
                indPos = pos(1) - 1 + indPos;
                [~,indNeg] = min(scale4(neg));
                indNeg = neg(1) - 1 + indNeg;
                if indPos > indNeg
                    [~,Tpos] = min(abs(scale3(indNeg:indPos)));
                    T(annT) = indNeg - 1 + Tpos;
                    nfirst = indNeg;
                    nlast = indPos;
                else
                    [~,Tpos] = min(abs(scale3(indPos:indNeg)));
                    T(annT) = indPos - 1 + Tpos;
                    nfirst = indPos;
                    nlast = indNeg;
                end
                %T beginning
                xTon = 0.5*abs(scale4(nfirst));
                neigh = floor((150/1000)*fs);
                if ~isempty(find(abs(scale4(nfirst - neigh : nfirst)) <= xTon, 1, 'last'))
                    Ton(annT) = nfirst - neigh + find(abs(scale4(nfirst - neigh : nfirst)) <= xTon, 1, 'last');
                else
                    Ton(annT) = nfirst;
                end
                
                %T ending
                xTend = 0.2*abs(scale4(nlast));
                if ~isempty(find(abs(scale4(nlast : nlast + neigh)) <= xTend, 1, 'first'))
                    Tend(annT) = nlast - 1 + find(abs(scale4(nlast : nlast + neigh)) <= xTend, 1, 'first');
                else
                    Tend(annT) = nlast;
                end
               annT = annT + 1;
            end
        end
    end
    %P wave detection and delineation
    neigh = floor((200/1000)*fs);
    win = floor((70/1000)*fs);
    closeneigh = floor((30/1000)*fs);
    if QRSon(i) - neigh > 0
        %Finding if P wave is present
        pos = QRSon(i) - neigh - 1 + find(scale4(QRSon(i)-neigh:QRSon(i)-win) > epsP);
        neg = QRSon(i) - neigh - 1 + find(scale4(QRSon(i)-neigh:QRSon(i)-win) < -1*epsP);
        if ~isempty(pos) && ~isempty(neg)
            gP = 0.05*max(abs(scale4(QRSon(i)-neigh:QRSon(i)-win)));
            if (max(scale4(pos)) > gP) && (min(scale4(neg))<-1*gP)
                [~,indPos] = max(scale4(pos));
                indPos = pos(1) - 1 + indPos;
                [~,indNeg] = min(scale4(neg));
                indNeg = neg(1) - 1 + indNeg;
                if indPos > indNeg
                    [~,Ppos] = min(abs(scale4(indNeg:indPos)));
                    Ppos = indNeg - 1 + Ppos;
                    [~,Pfinal] = min(abs(scale3(Ppos:Ppos+closeneigh)));
                    P(annP) = Ppos -1 + Pfinal;
                    
                else
                    [~,Ppos] = min(abs(scale4(indPos:indNeg)));
                    Ppos = indPos - 1 + Ppos;
                    [~,Pfinal] = min(abs(scale3(Ppos:Ppos+closeneigh)));
                    P(annP) = Ppos - 1 + Pfinal;
                   
                end
                %P beginning
                xPon = 0.5*abs(scale4(nfirst));
                neigh = floor((70/1000)*fs);
                [~,nfirst] = max(abs(scale3(P(annP)-neigh:P(annP))));
                nfirst = P(annP) - neigh - 1 + nfirst;
                if ~isempty(find(abs(scale3(nfirst - neigh : nfirst)) <= xPon, 1, 'last'))
                    Pon(annP) = nfirst - neigh + find(abs(scale3(nfirst - neigh : nfirst)) <= xPon, 1, 'last');
                else
                    Pon(annP) = nfirst;
                end
                
                %P ending
                xPend = 0.3*abs(scale4(nlast));
                [~,nlast] = max(scale3(P(annP):P(annP)+neigh));
                nlast = P(annP) - 1 + nlast;
                if ~isempty(find(abs(scale3(nlast : nlast + neigh)) <= xPend, 1, 'first'))
                    Pend(annP) = nlast - 1 + find(abs(scale3(nlast : nlast + neigh)) <= xPend, 1, 'first');
                else
                    Pend(annP) = nlast;
                end
               annP = annP + 1;
            end
        end
    end
    
end


%Plotting of results 

t = 0 : 1/fs : (length(signal)-1)/fs;

figure, 
plot(t,signal,t(R),signal(R),'xr',t(QRSon),signal(QRSon),'or',t(QRSend),signal(QRSend),'vr',t(T),signal(T),'xk',t(Ton),signal(Ton),'ok',t(Tend),signal(Tend),'vk',t(P),signal(P),'xg', t(Pon),signal(Pon),'og',t(Pend),signal(Pend),'vg');
legend('Signal','Rpeak','QRSbegin','QRSend', 'Tpeak','Tbegin','Tend', 'Ppeak','Pbegin','Pend'),
title('Detection Sample'),
xlabel('time (s)'),
ylabel('Amplitude'),
xlim([450 453]);


end
