function [phi1,phi2,omeg1,omeg2,q1,q2,Qcoef1,Qcoef2,norm1,norm2,dirin] =...
    co_fbtransf2(theta1,theta2,or,safr,ngrid)
% DAMOCO Toolbox, function CO_FBTRANSF2, version 17.01.11
% This high-level function performs all steps of the phase transformation at once,
% using the Fourier-based method.
% Given two protophases theta1 and theta2, the order of the Fourier expansion, or, 
% and the sampling frequency of the data, safr, this function computes:
%      the true phases phi1 and phi2, 
%      the true coupling functions q1 and q2 on a grid of size ngrid, 
%      the Fourier coefficients of the coupling functions, Qcoef1 and Qcoef2,
%      the frequencies (constant terms of coupling functions) omeg1 and omeg2,
%      the norms of the coupling functions, norm1 and norm2, 
%      and the directionality index dirin.
% Form of call: 
%      [phi1,phi2,omeg1,omeg2,q1,q2,Qcoef1,Qcoef2,norm1,norm2,dirin] = co_fbtransf2(theta1,theta2,or,safr,ngrid)
%      [phi1,phi2,omeg1,omeg2,q1,q2,Qcoef1,Qcoef2,norm1,norm2,dirin] = co_fbtransf2(theta1,theta2,or,safr)
% Input:
%       theta1:     protophase of the 1st oscillator
%       theta2:     protophase of the 2nd oscillator
%       or:         order of the Fourier expansion: or = 10 is a good first choice.
%                   To check and optimize this parameter: plot the  matrix of abs(Qcoef) 
%                   using co_plotcoef and verify that the coefficients at the boundary of the matrix are small. 
%                   If not, increase the parameter or. If the region, where the coefficients are nearly zero, 
%                   is large, reduce or.    
%       safr:       sampling frequency
%       ngrid:      size of the grid to compute q1, q2; default value is ngrid=100
% Output: 
%       phi1:       true phase of the 1st oscillator
%       phi2:       true phase of the 2nd oscillator
%       omeg1:      autonomous frequency (constant term of the coupling function) of the 1st oscillator
%       omeg2:      autonomous frequency (constant term of the coupling function) of the 2nd oscillator
%       q1:         coupling function of the 1st oscillator, computed on a grid 
%                   Note: It does not contain the constant term, i.e. the
%                   phase dynamics is given by dphi1/dt = omeg1 + q1(phi1,phi2).
%       q2:         coupling function of the 2nd oscillator, computed on a
%                   grid, without constant term.
%       Qcoef1:     the coefficients of the Fourier expansion of q1.
%       Qcoef1:     the coefficients of the Fourier expansion of q2.
%       norm1:      the norm of q1.
%       norm2:      the norm of q2.
%       dirin:      the index of directionality of interaction between both oscillators: 
%                   dirin = 1: unidirectional coupling from 1 to 2. 
%                   dirin = -1: unidirectional coupling from 2 to 1.
%                   dirin = 0: equally strong bidirectional coupling  
%
co_testproto(theta1, theta2);        % Testing protophases
pr_theta1 = co_fbtransf1(theta1);    % Preprocessing with univariate Fourier based transform
pr_theta2 = co_fbtransf1(theta2);    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%       WARNING,  if protophases are too strongly correlated
[dummy, maxind]=co_maxsync(pr_theta1, pr_theta2, 2*or); clear dummy;
if maxind > 0.6;
    disp(['Synchronisation index is ' num2str(maxind)]);
    disp('Protophases are strongly correlated, the results may not be reliable!');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%           FOURIER BASED METHOD
[Fcoef1, Fcoef2] = co_fexp2(pr_theta1, pr_theta2, or, safr);        % Computing coeffs of Fourier expansion of the protophase dynamics
[sigfc1, sigfc2] = co_fbsolv(Fcoef1, Fcoef2);                       % Computing coefficients of phase transformation
[phi1, phi2] = co_fbth2phi(pr_theta1, pr_theta2, sigfc1, sigfc2);   % Performing phase transformation
if nargin == 4; ngrid = 100; end;                                   % Default value of grid szize
[Qcoef1, Qcoef2, q1, q2] = co_fexp2(phi1, phi2, or, safr, ngrid);   % Computing coupling functions and their Fourier coeffs
omeg1 = real(Qcoef1(or+1,or+1));                    % Autonomous frequency is extracted from matrix Qcoef1 
omeg2 = real(Qcoef2(or+1,or+1));                    % Autonomous frequency is extracted from matrix Qcoef2
q1 = q1 - omeg1;                                    % Coupling function of the phase does not contain the autonomous frequency
q2 = q2 - omeg2;                                    % Coupling function of the phase does not contain the autonomous frequency
Qcoef1(or+1,or+1) = 0;                              % The coefficients describe the coupling function only, 
Qcoef2(or+1,or+1) = 0;                              % they do not contain the autonomous frequency omeg1
norm1 = co_fbnorm(Qcoef1);                          % Norms of the coupling functions
norm2 = co_fbnorm(Qcoef2);                      
dirin = co_dirin(norm1, norm2, omeg1,omeg2);        % Direction of couplinmg
end