function [qrs,feats]=multi_qrsdetect(signal,fs,fname)

% This function detects QRS peaks on ECG signals by taking the consensus of multiple algorithms, namely:

%    - gqrs (WFDB Toolbox)

%    - Pan-Tompkins (FECGSYN)

%    - Maxima Search (OSET/FECGSYN)

%    - matched filtering

%

% --

% 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/>.

qrs = cell(1,5); % CHANGE ME! Using 6 levels of detection

WIN_ALG = round(0.08*fs); % 100ms allow re-alignment of 80ms for gqrs

REF_WIN = round(0.17*fs);    % refractory window PT algorithm, no detection should occur here.

WIN_BLOCK = round(5*fs);    % window for PT algorithm

OLAP = round(1*fs);

MIRR = round(2*fs);           % length to mirror signal at beginning and end to avoid border effects

LEN = length(signal);

signal = [signal(1:MIRR); signal ; signal(LEN-MIRR+1:LEN)];

LEN = LEN + 2*MIRR;

%% gQRS

try

    recordName = ['gd_' fname];

    tm = 1/fs:1/fs:LEN/fs;

    wrsamp(tm',signal,recordName,fs,1)

    disp(['Running gqrs ' fname '...'])

    if isunix                 

        system(['gqrs -r ' recordName ' -f 0 -s 0']);

        disp(['gqrs ran ' fname '...'])

    else

       gqrs(recordName)

    end

    qrs{1} = rdann(recordName,'qrs');        

    disp(['read gqrs output ' fname '...'])

    [~,lag]=arrayfun(@(x) max(abs(signal(x-WIN_ALG:x+WIN_ALG))),qrs{1}(2:end-1));

    qrs{1}(2:end-1) = qrs{1}(2:end-1) + lag - WIN_ALG - 1;

    delete([recordName '.*'])

catch

    warning('gqrs failed.')

    delete([recordName '.*'])

    qrs{1} = [];

end

clear lag tm recordName

%% P&T algorithm (in windows of 5 seconds - 1 sec overlap)

try

    disp(['Running jqrs ' fname '...'])

    qrs{2} = jqrs(signal,REF_WIN,[],fs,[],[],0);

catch

    warning('Pan Tompkins failed.')

    qrs{2} = [];

end

clear qrs_tmp startp endp count

%% MaxSearch (mQRS)

qrsmax=qrs{double(length(qrs{2})>length(qrs{1}))+1}; % selecting as reference detector with most detections

try

    disp(['Running maxsearch ' fname '...'])

    if length(qrsmax)<5

        HR_PARAM = 1.3; % default value

    else

        HR_PARAM = fs./nanmedian(diff(qrsmax))+0.1;  % estimate rate (in Hz) for MaxSearch

    end

    qrs{3} = OSET_MaxSearch(signal,HR_PARAM/fs)';

catch

    warning('MaxSearch failed.')

    qrs{3} = [];

end

% Put all detectors on same sign (max or min) - shift by median lag

x = decimate(signal,5);    % just to reduce calculus

y = sort(x);    

flag = mean(abs(y(round(end-0.05*length(y)):end)))>mean(abs(y(1:round(0.05*length(y)))));

for i = 1:length(qrs)

    if flag

        [~,lag]=arrayfun(@(x) max(signal(x-WIN_ALG:x+WIN_ALG)),qrs{i}(3:end-2));

    else

        [~,lag]=arrayfun(@(x) min(signal(x-WIN_ALG:x+WIN_ALG)),qrs{i}(3:end-2));

    end

    qrs{i}(3:end-2) = qrs{i}(3:end-2) + lag - WIN_ALG - 1;

end

clear HR_PARAM qrsmax

%% Matched filtering QRS detection

try

    disp(['Running matchedfilter ' fname '...'])

    [btype,templates,tempuncut] = beat_class(signal,qrs{3},round(WIN_ALG));

    [qrs{4},ect_qrs]=matchedQRS(signal,templates,qrs{3},btype);

    matchedf = 1; % status signal used as feature

catch

    warning('Matched filter failed.')

    disp('Matched filter failed.')

    qrs{4} = [];

    ect_qrs = [];

    matchedf = 0;

end

%% Consensus detection

try

    consqrs = kde_fusion(cell2mat(qrs(1:4)'),fs,length(signal)); % 2x gqrs

    consqrs(diff(consqrs)<REF_WIN) = []; % removing double detections

    [~,lag]=arrayfun(@(x) max(abs(signal(x-WIN_ALG:x+WIN_ALG))),consqrs(2:end-2));

    consqrs(2:end-2) = consqrs(2:end-2) + lag - WIN_ALG - 1;

    % Put all detectors on same sign (max or min) - shift by median lag

    if flag

        [~,lag]=arrayfun(@(x) max(signal(x-WIN_ALG:x+WIN_ALG)),consqrs(3:end-2));

    else

        [~,lag]=arrayfun(@(x) min(signal(x-WIN_ALG:x+WIN_ALG)),consqrs(3:end-2));

    end

    qrs{5} = consqrs(3:end-2) + lag - WIN_ALG - 1;

catch

    disp('Consensus failed.')

    qrs{5} = qrs{3};

end

% Remove border detections

qrs = cellfun(@(x) x(x>MIRR&x<LEN-MIRR)-MIRR,qrs,'UniformOutput',false);

signal = signal(MIRR+1:end-MIRR); % just for plotting

ect_qrs = ect_qrs(ect_qrs>MIRR&ect_qrs<(LEN-MIRR)) - MIRR; 

clear lag flag consqrs

%% In case ectopic beats are present (consensus QRS locations is ignored, matched filter assumes)

% try

%     if ~isempty(ect_qrs)

%         % Remove ectopic beats

%         cons = qrs{4};

%         [~,dist] = dsearchn(cons,ect_qrs);

%         insbeats = ect_qrs(dist>REF_WIN); % inserting these beats

%         

%         cons = sort([cons; insbeats]);

%         cons(arrayfun(@(x) find(cons==x),insbeats)) = NaN; % insert NaN on missing beats

%         cons = round(inpaint_nans(cons)); % Interpolate with simple spline over missing beats

%         cons(cons<1|cons>length(signal)) = [];

%         cons = sort(cons);

%         cons(diff(cons)<REF_WIN) = [];

%         qrs{6} = cons;

%         

%         %    % visualise

% %             plot(signal)

% %             hold on

% %             plot(qrs{4},2.*ones(size(qrs{4})),'ob')

% %             plot(ect_qrs,2.*ones(size(ect_qrs)),'xr')

% %             plot(qrs{6},3.*ones(size(qrs{6})),'md')

% %             legend('signal','matchedf','ectopic','interpolated')

%         %     close

%         

%     else

%         qrs{6} = qrs{5};

%     end

% catch

%     warning('Figuring out ectopic failed.')

%     qrs{6} = qrs{5};

%     ect_qrs = [];

% end

%% Generating some features

% Amplitude around QRS complex

normsig = signal./median(signal(qrs{end})); % normalized amplitude

feat_amp = var(normsig(qrs{end})); % QRS variations in amplitude

feat_amp2 = std(normsig(qrs{end})); % QRS variations in amplitude

feat_amp3 = mean(diff(normsig(qrs{end}))); % QRS variations in amplitude

feats = [feat_amp feat_amp2 feat_amp3];

%% Plots for sanity check

% subplot(2,1,1)

% plot(signal,'color',[0.7 0.7 0.7])

% hold on

% plot(qrs{1},signal(qrs{1}),'sg')

% plot(qrs{2},signal(qrs{2}),'xb')

% plot(qrs{3},signal(qrs{3}),'or')

% plot(qrs{4},signal(qrs{4}),'dy')

% plot(qrs{5},1.3*median(signal(qrs{5}))*ones(size(qrs{5})),'dm')

% legend('signal','gqrs','jqrs','maxsearch','matchedfilt','kde consensus')

% subplot(2,1,2)

% plot(qrs{6}(2:end),diff(qrs{6})*1000/fs)

% ylabel('RR intervals (ms)')

% ylim([250 1300])

% close

end

function det=kde_fusion(qrs,fs,dlength)

% Function uses Kernel density estimation on fusing multiple detections

%

% Input

%    qrs      Array with all detections

%    fs       Sampling frequency

% dlength     Signal length

%

%

qrs = sort(qrs);

w_std = 0.10*fs;    % standard deviation of gaussian kernels [ms]

pt_kernel = round(fs/2);    % half window for kernel function [samples]

%% Calculating Kernel Density Estimation

% preparing annotations

peaks = hist(qrs,1:dlength);

% kde (adding gaussian kernels around detections)

kernel = exp(-(-pt_kernel:pt_kernel).^2./(2*w_std^2));

% calculating kernel density estimate

kde = conv(peaks,kernel,'same');

%% Decision

% Parameters

min_dist = round(0.2*fs);    % minimal distance between consecutive peaks [ms]

th = max(kde)/3;      % threshold for true positives (good measure is number_of_channels/3)

% Finding candidate peaks

wpoints = 1+min_dist:min_dist:dlength; % start points of windows (50% overlap)

wpoints(wpoints>dlength-2*min_dist) = []; % cutting edges

M = arrayfun(@(x) kde(x:x+2*min_dist-1)',wpoints(1:end), ...

    'UniformOutput',false);  % windowing signal (50% overlap)

M = cell2mat(M);

% adding first segment

head = [kde(1:min_dist) zeros(1,min_dist)]';

M = [head M];

% adding last segment

tail = kde((wpoints(end)+2*min_dist):dlength)';

tail = [tail; zeros(2*min_dist-length(tail),1)];

M = [M tail];

[~,idx] = max(M);   % finding maxima

idx = idx+[1 wpoints wpoints(end)+2*min_dist]; % absolute locations

idx(idx < 0) = [];

idx(idx > length(kde)) = [];

i = 1;

while i<=length(idx)

    doubled = (abs(idx-idx(i))<min_dist);

    if sum(doubled) > 1

        [~,idxmax]=max(kde(idx(doubled)));

        idxmax = idx(idxmax + find(doubled,1,'first') - 1);

        idx(doubled) = [];

        idx = [idxmax idx];

        clear doubled idxmax

    end

    i = i+1;

end

% threshold check

idx(kde(idx)<th) = [];

det = sort(idx)';

end