--- a
+++ b/featurebased-approach/subfunctions/stackbeats.m
@@ -0,0 +1,72 @@
+function [M,RES] = stackbeats(signal,qrs,T_LENGTH,NB_BINS)
+% Stack and wraps beats onto a matrix M
+%  RES = time reference to resample from bins to time
+%
+% --
+% 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/>.
+
+% Parameter
+if size(signal,1) < size(signal,2), signal = signal'; end
+
+phase = FECGx_kf_PhaseCalc(qrs,length(signal));
+ini_cycles = find(phase(2:end)<0&phase(1:end-1)>0)+1; % start of cycles
+cycle_len = diff(ini_cycles); % distance between cycles
+end_cycles = ini_cycles(1:end-1)+cycle_len-1; % start of cycles
+meanPhase = linspace(-pi,pi,NB_BINS);
+RES =  2.*round((median(cycle_len)+1)/2)-1;
+% stacking cycles
+cycle = arrayfun(@(x) interp1(phase(ini_cycles(x):end_cycles(x)),...
+    signal(ini_cycles(x):end_cycles(x)),meanPhase,'spline'),...
+    1:length(ini_cycles)-1,'UniformOutput',0);
+M = cell2mat(cycle')';
+
+% re-aligning beats
+i = 1;
+while i<=size(M,2)
+    lag = finddelay(mean(M'),M(:,i)); % finds lag on the beat
+    if abs(lag) > T_LENGTH    % only allows 20% shift
+        M(:,i) = [];            % otherwise ignore beat
+        continue
+    end
+    if size(M,2)>2 && lag < 0
+        M(end-abs(lag):end,i) = 0;
+    elseif size(M,2)>2
+        M(1:lag,i) = 0;
+    end
+    M(:,i)=circshift(M(:,i),[-lag,1]);
+    i = i+1;
+end
+
+
+end