function [ECGmean,ECGsd,meanPhase] = MeanECGExtraction(x,phase,bins,flag)
%
% [ECGmean,ECGsd,meanPhase] = MeanECGExtraction(x,phase,bins,flag)
% Calculation of the mean and SD of ECG waveforms in different beats
%
% inputs:
% x: input ECG signal
% phase: ECG phase
% bins: number of desired phase bins
% flag
%     1: aligns the baseline on zero, by using the mean of the first 10%
%     segment of the calculated mean ECG beat
%     0: no baseline alignment
%
% outputs:
% ECGmean: mean ECG beat
% ECGsd: standard deviation of ECG beats
% meanPhase: the corresponding phase for one ECG beat
%
%
% Open Source ECG Toolbox, version 2.0, March 2008
% Released under the GNU General Public License
% Copyright (C) 2008  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- LIS-INPG, Grenoble, France
% reza.sameni@gmail.com

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

meanPhase = zeros(1,bins);
ECGmean = zeros(1,bins);
ECGsd = zeros(1,bins);

I = find( phase>=(pi-pi/bins)  | phase<(-pi+pi/bins) );
if(~isempty(I))
    meanPhase(1) = -pi;
    ECGmean(1) = mean(x(I));
    ECGsd(1) = std(x(I));
else
    meanPhase(1) = 0;
    ECGmean(1) =0;
    ECGsd(1) = -1;
end
for i = 1 : bins-1;
    I = find( phase >= 2*pi*(i-0.5)/bins-pi & phase < 2*pi*(i+0.5)/bins-pi );
    if(~isempty(I))
        meanPhase(i + 1) = mean(phase(I));
        ECGmean(i + 1) = mean(x(I));
        ECGsd(i + 1) = std(x(I));
    else
        meanPhase(i + 1) = 0;
        ECGmean(i + 1) = 0;
        ECGsd(i + 1) = -1;
    end
end
K = find(ECGsd==-1);
for i = 1:length(K),
    switch K(i)
        case 1
            meanPhase(K(i)) = -pi;
            ECGmean(K(i)) = ECGmean(K(i)+1);
            ECGsd(K(i)) = ECGsd(K(i)+1);
        case bins
            meanPhase(K(i)) = pi;
            ECGmean(K(i)) = ECGmean(K(i)-1);
            ECGsd(K(i)) = ECGsd(K(i)-1);
        otherwise
            meanPhase(K(i)) = mean([meanPhase(K(i)-1),meanPhase(K(i)+1)]);
            ECGmean(K(i)) = mean([ECGmean(K(i)-1),ECGmean(K(i)+1)]);
            ECGsd(K(i)) = mean([ECGsd(K(i)-1),ECGsd(K(i)+1)]);
    end
end

if(flag==1)
    ECGmean = ECGmean - mean(ECGmean(1:ceil(length(ECGmean)/10)));
end