% function [p,Im,th,w,s,md,mx,Y1,Y2,Y3,g,alfa] = detqrs3(x,Fs,n,Ft,pth)
function [QRS] = detqrs3(x,Fs,n,Ft,pth)

% Funzione che effettua la detezione dei punti fiduciari dei complessi QRS
% presenti in un segnale ECG in ingresso
% x segnale ecg di cui effettuare la detezione dei complessi qrs
% Fs freq. di campionamento del segnale
% n ordine del filtro passa - basso
% Ft freq. di taglio del filtro passa - basso
% pth percentuale dell'ampiezza del qrs con ampiezza max, su cui viene
% definita la soglia

if nargin < 3
    n = 4;
    Ft = 60;
    pth = 50;
end

wt = Ft/(Fs/2);   % blocco filtraggio passa - basso
B = fir1(n,wt);
Y1 = filter(B,1,x);
alfa = Y1;

for i=1:4
    Y1(i) = mean(Y1);
end


for i = 2:1:length(x)
    Y2(i) = (Y1(i)-Y1(i-1))*100;
end
for i = 1:1:length(Y2)
    if Y2(i) > 0
        Y3(i) = Y2(i)^2;
    else Y3(i) = 0;
    end
end

th = max(Y3)*pth/100;
mx = max(Y3);
md = mean(Y3);
p = zeros(1,length(Y1));
w = floor(0.16*Fs);
s = floor(0.3*Fs);
for i = w+1:1:length(Y1)-w
    if i > w+1
        if Y1(i) > Y1(i-w:i-1)
            if Y1(i) > Y1(i+1:i+w)
                if i <= s  
                    if p(i-w:i) == 0
                        if max(Y3(i-w:i+w)) > th
                            p(i) = 1;
                        end
                    end
                elseif i > s
                    if p(i-s:i) == 0
                        if max(Y3(i-w:i+w)) > th
                            p(i) = 1;
                        end
                    end
                end
            end
        end
    elseif i == w+1
        if Y1(i) > Y1(i-w:i-1)
            if Y1(i) > Y1(i+1:i+w)
                if p(i-w:i) == 0
                    if Y3(i-w:i+w) > th
                        p(i) = 1;
                    end
                end
            end
        end
    end
end

if p(1:s) == 0;
    for i = 2:s
        if Y1(i) > Y1(1:i-1)
            if Y1(i) > Y1(i+1:s)   
                %if Y3(i) > th
                    p(i) = 1;
                    break
                    %end
            end
        end
    end
end

g = zeros(1,length(p));
for i = 1:length(p)
    if p(i) == 1
        g(i) = 3/2*max(Y1);
    end
    if p(i) == 0
        g(i) = min(Y1);
    end
end

QRS = find(p);

% t=1:length(Y1);
% plot(t,x,'b',t,g,'r:')
% pause
% close
                    
        