--- a
+++ b/featurebased-approach/subfunctions/matchedQRS.m
@@ -0,0 +1,114 @@
+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
+
+