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