function co_example4
% DAMOCO Toolbox, function CO_EXAMPLE4, version 16.01.11
%
% This function illustrates the bivariate analysis.
% It reads from the file co_vdp2uni.mat the input data,
% reconstructs the phase model and quantifies the coupling.
% The data are from two unidirectionally coupled van der Pol 
% oscillators, see manual for parameters.
% This function differs from co_example2 by the method it uses:
% here the nonlinear equation system for unknown transformation 
% functions is solved by iterations.
%
pi2=pi+pi;
load co_vdp2uni.mat;
% The previous line loads to wokspace three variables:
%    signals x1 and x2, 
%    sampling frequency fsamp. 
disp(' ');disp(' ');disp(' ');disp(' ');
disp('--------  Starting co_example4  -----------');
disp('--  Plotting a piece of bivariate data  --');
figure(1); time=(1:1000)/fsamp;     % plotting a piece of data
plot(time,x1(1:1000),time,x2(1:1000),'-r'); 
xlabel('time'); ylabel('x_1,x_2'); 
title('Data from 2 coupled van der Pol oscillators (not all shown)');
xlim([0 time(end)]);
%
%%%%%%%%%%%%%%   step 1: computation of protophases  %%%%%%%%%%%%%%%%%%%%%
disp('--  Computung and plotting protophases via Hilbert transform  --');
%
% In simple cases like this one, the protophases can be computed via
% the Hilbert transform
[theta1,minampl]= co_hilbproto(x1,2); % Embedding is plotted as figure 2
% We use the default values for other parameters
disp(['Minimal amplitude over average amplitude is ' num2str(minampl)])
% Graphic output and computation of the minimal instantaneous amplitude
% check whether the trajectory does not come close to the origin,
% i.e. that the protophase is well-defined
theta2= co_hilbproto(x2); % Second protophase: without plot and 
% explicit check. Note that the program automatically checks the minamp 
% parameter, even if it is skipped in the call, and issues a warning, 
% if the data are not good enough for the Hilbert transform.
figure(3); plot(theta1,theta2,'r.','MarkerSize',1);
xlabel('\theta_1'); ylabel('\theta_2'); axis square; 
title('Raw protophases'); axis tight;
%
% For further comparison we plot the coupling functins, constructed 
% from raw protophases
%
Forder=10;      % order of Fourier series is 10
ngrid=50;       % grid size is 50
[Fcoef1,Fcoef2,f1,f2]=co_fexp2(theta1,theta2,Forder,fsamp,ngrid); 
co_plot2cplf(f1,f2,0,4,'\theta_1','\theta_2','F_1(\theta_1,\theta_2)',...
    'F_2(\theta_2,\theta_1)',...
    'Coupling function from raw protophases');
%
%%%%%%   step 2: univariate transformation as preprocessing   %%%%%%
%
[theta1_1,arg,sigma1] = co_fbtransf1(theta1);% default parameters are used
[theta2_1,arg,sigma2] = co_fbtransf1(theta2);
figure(5); plot(arg,sigma1,arg,sigma2,'r-'); % transformation functions are 
xlabel('\theta');ylabel('\sigma_{1,2}');     % plotted for illustration  
title('Univariate transformation functions'); xlim([0 pi2]);
figure(6); plot(theta1_1,theta2_1,'r.','MarkerSize',1);
xlabel('\phi_1'); ylabel('\phi_2'); axis square; 
title('Phases after univariate transformations'); axis tight;
%
%%%%%%%%%%%%     step 3: bivariate transformation           %%%%%%%%%%%%      
%
disp('--  Reconstructing phase dynamics from 1D transformed protophases  --');
[Fcoef1,Fcoef2,f1,f2]=co_fexp2(theta1_1,theta2_1,Forder,fsamp,ngrid); 
% These are coupling functions constructed from 1D-transformed protophases,
% They are plotted for further comparison.
co_plot2cplf(f1,f2,0,7,'\theta_1','\theta_2','F_1(\theta_1,\theta_2)',...
    'F_2(\theta_2,\theta_1)',...
    'Coupling function from 1D-transformed protophases');
%
% Now comes the bivariate transformation itself, this time iteration technique
% is used to solve nonlinear equations.
% Number of iterations is 3, results of iterations are shown
niter=3;   showiter=1; % intermediate results are show if showiter >0
[q1,q2,omeg1,omeg2,Nrmq1,Nrmq2,sigma1,sigma2]=co_itersolv(f1,f2,niter,showiter);
%
figure(8); plot(arg,sigma1,arg,sigma2,'-r'); % transformation functions are 
xlabel('\theta');ylabel('\sigma_{1,2}');     % plotted for illustration  
title('Bivariate transformation functions');  xlim([0 pi2]);
%
disp('Constant terms of the coupling functions');
disp('   (estimates of natural frequencies) are:');
disp(['omega1 = ' num2str(omeg1) ', omega2 = ' num2str(omeg2)]);
co_plot2cplf(q1,q2,0,9,'\phi_1','\phi_2','\omega_1+Q_1(\phi_1,\phi_2)',...
    '\omega_2+Q_2(\phi_2,\phi_1)',...
    'True coupling functions after bivariate transformation');
%
phi1=co_gth2phi(theta1_1,sigma1);   % final transform of 1D-transformed protophases
phi2=co_gth2phi(theta2_1,sigma2);
%
%%%%%%%%%%    step 4: quantification of coupling  %%%%%%%%%%%%%%%%%%%%
%
siproto=co_sync(theta1,theta2,1,1);
disp(['Synchronization index from protophases ' num2str(siproto)]); 
siphase=co_sync(phi1,phi2,1,1);
disp(['Synchronization index from phases ' num2str(siphase)]); 
%
dirin=co_dirin(Nrmq1,Nrmq2,omeg1,omeg2);
disp(['Directionality index ' num2str(dirin)]); 
%
disp('--------  End of co_example2  --------');
end