function [qrs,ect_qrs] = matchedQRS(signal,templates,qrsref,btype)

% Matched filtering for QRS detection

%

%

%

% --

% ECG classification from single-lead segments using Deep Convolutional Neural 

% Networks and Feature-Based Approaches - December 2017

% 

% Released under the GNU General Public License

%

% Copyright (C) 2017  Fernando Andreotti, Oliver Carr

% University of Oxford, Insitute of Biomedical Engineering, CIBIM Lab - Oxford 2017

% fernando.andreotti@eng.ox.ac.uk

%

% 

% For more information visit: https://github.com/fernandoandreotti/cinc-challenge2017

% 

% Referencing this work

%

% Andreotti, F., Carr, O., Pimentel, M.A.F., Mahdi, A., & De Vos, M. (2017). 

% Comparing Feature Based Classifiers and Convolutional Neural Networks to Detect 

% Arrhythmia from Short Segments of ECG. In Computing in Cardiology. Rennes (France).

%

% Last updated : December 2017

% 

% This program is free software: you can redistribute it and/or modify

% it under the terms of the GNU General Public License as published by

% the Free Software Foundation, either version 3 of the License, or

% (at your option) any later version.

% 

% This program is distributed in the hope that it will be useful,

% but WITHOUT ANY WARRANTY; without even the implied warranty of

% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

% GNU General Public License for more details.

% 

% You should have received a copy of the GNU General Public License

% along with this program.  If not, see <http://www.gnu.org/licenses/>.

THRES = 0.7;  % correlation threshold, relative to priori knowledge on QRS

THRESgof = 0.4; % NMSE threshold

WIN = round(2*length(templates{1}));

qrs = cell(size(templates));

xcors = cell(size(templates));

ect_qrs = [];

%% Find new peaks by cross correlation

for i = 1:length(templates)

    padt = [templates{i}; zeros(length(signal)-length(templates{i}),1)];

    xcor = xcorr(padt,signal)/(norm(padt)*norm(signal));

    xcor = xcor(length(signal):-1:1); xcor = circshift(xcor,[round(length(templates{i})/2) 1]);

    xcor(xcor<0) = 0;

    thresh=nanmedian(THRES*xcor(qrsref));

    [~,qrstmp] = findpeaks(xcor,'MinPeakDistance',round(0.8*WIN),'MinPeakProminence',thresh);

%     qrstmp = sort([qrstmp; qrsref(btype == i)]); % Rescue previously detected abnomallies

%     qrstmp(diff(qrstmp)<WIN) = [];

    qrstmp(qrstmp<length(templates{i})/2|qrstmp>(length(signal)-length(templates{i})/2)) = []; % remove extremities

    segs = arrayfun(@(x) signal(x-floor(length(templates{i})/2):x+floor(length(templates{i})/2)),qrstmp,'UniformOutput',0);

    gof = cellfun(@(x)  1 - sum((x -templates{i}).^2)/sum((templates{i}-mean(templates{i})).^2),segs); % calculates goodness of fit for each beat agains template

    qrstmp = qrstmp(gof>THRESgof);

    qrs{i} = qrstmp;

    xcors{i} = xcor(qrstmp).*gof(gof>THRESgof);

    clear xcortmp qrstmp

end

%% Separate beats by type (avoid double peaks)

%== Removing peaks of low amplitude (likely noise)

% if length(templates) > 1

%     if 0.3*abs(median(signal(qrs{1}))) > abs(median(signal(qrs{2})))

%         qrs(2) = [];

%         templates(2) = [];

%     elseif 0.3*abs(median(signal(qrs{2}))) > abs(median(signal(qrs{1})))

%         qrs{1} = qrs{2};        

%         templates{1} = templates{2};

%         qrs(2) = [];

%         templates(2) = [];

%     end            

% end

%== Highest correlation/gof is kept.

if length(templates) > 1

        % matching closest indices

        if length(qrs{1}) > length(qrs{2})

            qrs1 = qrs{1};

            qrs2 = qrs{2};

        else

            qrs2 = qrs{1};

            qrs1 = qrs{2};

            tmp = xcors{1};

            xcors{1} =xcors{2};

            xcors{2} = tmp;

        end

        [A,dist] = dsearchn(qrs1,qrs2);         % closest ref for each point in test qrs

        keepA = (dist<WIN);

        if any(keepA)

            tmpmat = [A(keepA) find(keepA)];

            [~,idx]=max([xcors{1}(tmpmat(:,1)), xcors{2}(tmpmat(:,2))],[],2);

            qrs1(tmpmat(idx ~= 1,1)) = [];

            qrs2(tmpmat(idx ~= 2,2)) = [];

        end

        if numel(qrs1)>numel(qrs2)

            qrs = qrs1;          

            ect_qrs = qrs2;

        else

            qrs = qrs2;          

            ect_qrs = qrs1;

        end

%     end

else

    qrs = qrs{1};

end