%% find the QRS complexes
%  This method will try a different number of methods to find the QRS
%  complexes in case one fails
%
%  Inputs
%   -   ecg     The ECG signal
%   -   fs      The sampling frequency in Hz
%  Outputs
%   -   QRS     The QRS complex locations
%
%  Pietro Bonizzi
%  Joel Karel
%  Maastricht University
%  2017
function QRS = findQRScomplexes(ecg,fs)

% qrs detection
[QRS] = qrs_detect2(ecg',0.25,0.6,fs);
% plot(tm,ecg,tm(QRS),ecg(QRS),'or')
if isempty(QRS)
    [~, ~, ~, QRS, ~, ~]  = rpeakdetect(ecg,fs);
end
if length(QRS)<length(ecg)/fs*0.66
    if QRS(1)>fs
    [QRS1] = qrs_detect2(ecg(1:QRS(1)-1)',0.25,0.6,fs);
    else 
        QRS1 = [];
    end
    if length(ecg)-(QRS(end)+round(fs*0.2))>fs
    [QRS3] = qrs_detect2(ecg(QRS(end)+round(fs*0.2):end)',0.25,0.6,fs);
    else
        QRS3 = [];
    end
    QRS = [QRS1 QRS QRS3+QRS(end)+round(fs*0.2)-1];
end

QRSfinal = [];
for i = 1:length(QRS)-1
    clear QRStemp
    if QRS(i+1)-QRS(i) > fs*1.5
         [QRStemp] = qrs_detect2(ecg(QRS(i)+fs/10:QRS(i+1)-fs/10)',0.25,0.6,fs);
         QRSfinal = [QRSfinal QRS(i) QRStemp+QRS(i)+fs/10-1];
    else 
        QRStemp = [];
        QRSfinal = [QRSfinal QRS(i) QRStemp];
    end
end

QRSfinal = [QRSfinal QRS(end)];

clear QRS
QRS = QRSfinal;

RRcheck = diff(QRS);
% outlier detection
Q1 = prctile(RRcheck,5)-1;
Q3 = prctile(RRcheck,95)+1;
IQR = Q3-Q1;
% upperbound = Q3+10*IQR;
% lowerbound = Q1-10*IQR;
upperbound = Q3;
lowerbound = Q1;
% l(peakFWA>upperbound) = [];
QRS(find(RRcheck>upperbound)+1) = [];
RRcheck(find(RRcheck>upperbound)) = [];
% l(peakFWA<lowerbound) = [];
QRS(find(RRcheck<lowerbound)+1) = [];

QRS(find(diff(QRS)<fs*0.2)+1) = [];

% x = 1:length(ecg);
% plot(x,ecg,x(QRS),ecg(QRS),'or')
