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

