function [phi, Cav] = co_avcyc(x, theta, N, PL, alpha)
% DAMOCO Toolbox, function CO_AVCYC, version 18.05.14
%
% Given an oscillatory periodic time series, x, and a first estimate 
% of its protophase, theta, the function yields a new protophase, phi, 
% based on the projetion of the embedded time series on the average cycle.
%
% Form of call: 
%                   [phi, Cav] = co_avcyc(x, theta, N, alpha, PL)
% Input:            x:          oscillatory time series
%                   theta:      first estimate of the protophase
%                   N:          maximal order of the Fourier expansion
%                   alpha:      defines the relative weight in the 
%                               minimization procedure: 
%                               (1) alpha=0: the new protophase, phi, is
%                                   based on minimizing the distance to the
%                                   average cycle only. alpha =0 is default
%                                   value.
%                               (2) weight alpha quantifies how the 
%                                   difference of phi to the first estimate, 
%                                   theta, is used as a second constraint.
%                                   Hence, for very large alpha,  phi=theta.
%                                   This is useful in case of time series
%                                   possessing flat regions where the change
%                                   of the signal is almost on the noise level. 
%                                   In these regions the signal does not 
%                                   contain information on the evolution of 
%                                   the phase. A method to optimize alpha is
%                                   described in the supplementary material 
%                                   in Ref. [5]. 
%                   PL:         If PL==1, some results are plotted:
%                                     (1) average cycle
%                                     (2) the coefficients Cav
%                                     (3) the derivative of new protophase, phi (green), 
%                                         and of the first estimate (red).
% Output:           phi:        the new protophase, based on projection on the average cycle. 
%                   Cav:        the coefficients of Fourier expansion of the average cycle.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin <5; alpha = 0; end;                      % default value  of alpha
if nargin < 4; PL=0; end;                          % default value of PL
S = size(x); if S(1)> S(2); x=x.'; end             % Check for correct format
S = size(theta); if S(1)> S(2); theta=theta.'; end % Check for correct format
x = x / max(x);                                    % Normalize amplitude of the time series
H = hilbert(x);                                    % Computing 2-dimensional embedding 
%                                                    by means of Hilbert transform 
H = H(1000:end-1000);
theta = theta(1000:end-1000);
[Cav] = av_comp(H, theta, N, PL);               % Computing coefficients of the average 
%                                                    cycle, Cav.
S=size(Cav); if S(1)> S(2); Cav=Cav.'; end         % Check for correct format
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% Computation of new protophase, phi, by means of
%%%%%%%%%%%%%%%%%%% projection on the average cycle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
phi = zeros(1,length(theta));
DIS=mean(gradient(theta));                          % defining size of search region
for n = 1:length(theta)                      % Estimating new protophase, phi, by minimizing error
    phi(n) = fminbnd(@ERav, theta(n)-DIS, theta(n) + DIS,[], Cav, H(n), theta(n), alpha);   
end
if PL ==1;
    figure;
        plot(mod(phi,2*pi), gradient(phi),'g.','MarkerSize',4);
        hold on
            plot(mod(theta,2*pi), gradient(theta),'r.','MarkerSize',4);
        hold off
end        
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%  SUBFUNCTION   [Cav] = co_av_comp(H, theta, N, PL)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Cav] = av_comp(H, theta, N, PL)
% DAMOCO Toolbox, function CO_AV_COMP, version 06.05.14
% Given a 2-dimensional complex embedding by means of Hilbert transform, H, and an estimate of the
% protophase, theta, the function yields the coefficients of the average cycle, Cav.
% Form of call: 
%          [Cav] = co_av_comp(H, theta, N, PL)
% Input:   H:       2-dimensional embedding by means of the Hilbert transform
%          theta:   estimate of the protophase
%          N:       maximal order of Fourier expansion,
%          PL:      if PL ==1 the results are plotted
% Output:  Cav      Fourier coefficients of the average cycle
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=size(theta); if S(1) > S(2); theta = theta'; end  % Check for correct format
S=size(H); if S(1) > S(2);  H = H.'; end
Cav = zeros(N+1,1);
d_theta = gradient(theta);                          % Computing derivative of theta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%     Computing the coefficients Cav
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n = 0 : N 
    Cav(n+1) =  sum(H.*exp(-1i*n*theta).*d_theta)/(theta(end)-theta(1)); 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%  Plotting results if PL ==1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if PL ==1;
    Y= zeros(1,10000);
    p2=2*pi*(0:9999)/10000;
    for n = 0 : N
        Y = Y + Cav(n +1)*exp(1i*n*p2);
    end
    figure;
    plot(H,'.','MarkerSize', 4);
    hold on;
        plot(Y,'.r','MarkerSize',4);
    hold off
    figure; plot(abs(Cav))
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% ERRORFUNCTION  ERav(phi, Cav, H, theta, alpha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ER] = ERav(phi, Cav, H, theta, alpha)
N = (0:length(Cav)-1).';
Z = Cav * exp(1i*N*phi);
ER = (abs(H-Z).^2 + alpha* abs(exp(1i*phi)-exp(1i*theta) ).^2);
end