function [COUP, NORM, OMEGA]=co_tricplfan(Qcoef1, Qcoef2, Qcoef3, meth, thresh)
% DAMOCO Toolbox, function CO_TRICPLFAN, version 12.03.14
%
% This function is a tool to analyze the coupling structure of
% a three oscillator network. Given the Fourier coefficients
% Qcoef1, Qcoef2, Qcoef3 of the coupling functions, obtained by co_fcpltri,
% co_tricplfan returns the coupling structure in the form of matrix COUP.
%
% The n-th row contains measures of the effect ON the n-th oscillator.
% The m-th column contains measures of the effect CAUSED BY the m-th oscillator.
% For example, COUP(1,2) quantifies effect of node 2 on node 1, 2-->1
% 0 COUP(1,2) COUP(1,3)
% COUP(2,1) 0 COUP(2,3)
% COUP(3,1) COUP(3,2) 0
% The diagonal is by definition empty, since it would present the effects of an
% oscillator on itself; this place is however used as discussed below.
%
% COUP is computed by two different techniques, determined by the parameter 'meth'
%
% meth = 1 or meth = 2 : PARTIAL NORM
%
% A non-diagonal element COUP(n,m) is partial norm, i.e. norm of the Fourier terms of the coupling
% function Qcoefn, depending only on phases n,m (see Ref [4] for details).
% A diagonal value COUP(n,n) is the norm of the terms Qcoefn(n,n~0,k~0), i.e. of the terms,
% depending on both external phases.
% If diagonal terms are not small, the coupling cannot be explained by pairwise interaction
% only (see Ref [4]).
%
%
% meth = 3 or meth = 4 : PARTIAL DERIVATIVE
%
% The coefficients COUP are computed by means of the norm of the partial
% derivative with respect to one of external phases. For example
% the effect of the 2nd oscillator on the first is quantified as:
% COUP = sqrt(sum( m^2 Qcoef(n,m,k)^2) ).
% The drawback of this method is that noise is enhanced for large order n,m,k.
% The diagonal terms here are zero.
%
% co_tricplfan also returns the norms of the all coupling functions in the 3-dimensional
% array NORM and estimates of the autonomous frequencies of the oscillators in the 3-dimensional
% array OMEGA.
%
% To reduce numerical / noise effects, a threshold 'thresh' can be chosen in terms of percent
% of the maximal absolute value of the coefficients Qcoef:
% all coefficients Qcoef with an absolute value smaller than thresh/100*max(abs(Qcoef))
% are set to zero. The default value is thresh=2%.
%
% From COUP, which quantifies the interaction of oscillators in absolute values, a relative
% quantity can be defined by dividing COUP by the autonomous frequency OMEGA of the affected
% oscillator:
% COUP_rel = COUP / OMEGA
% This measure quantify the strength of coupling in relation to autonomous frequency;
% hence, it quantify the tendency to synchronization.
% Output COUP_rel can be chosen by setting meth = 2 or meth = 4.
%
% Form of call: [COUP, NORM, OMEGA]=co_3osccoup(Qcoef1, Qcoef2, Qcoef3) [all parameters = default]
% [COUP, NORM, OMEGA]=co_3osccoup(Qcoef1, Qcoef2, Qcoef3, meth): specifying method
% [COUP, NORM, OMEGA]=co_3osccoup(Qcoef1, Qcoef2, Qcoef3, meth, thresh): specifying
% both method and threshold
%
% INPUT:
% Qcoef1: matrix of Fourier coefficients of the coupling function
% of the 1st oscillator, produced by co_fcpltri
% Qcoef2: .. of the 2nd oscillator..
% Qcoef3: .. of the 3rd oscillator..
% meth: (1) (default)= partial norm, (2)= partial norm / omega, (3) = partial derivative,
% (4) = partial derivative / omega
% thresh: threshold for reduction of noise effects (also of computational noise) in terms
% of percent of the maximal absolute value of the coefficients Qcoef:
% thresh = 10 is good choice for dominant components only, the default
% value is thresh = 2%; typically this removes numerical artifacts, while keeping
% all other terms. Feel free to vary this parameter [0,..100], to remove noise, etc..
%
% OUTPUT:
% COUP matrix describing the coupling structure using the chosen method meth
% NORM: norms of all coupling functions
% OMEGA: the autonomous frequencies of the 3 oscillators
%
if nargin == 4;
thresh = 2; % defining default value thresh=2%
end;
if nargin == 3;
meth = 1; % default method =1: partial norm
thresh = 2; % defining default value thresh=2;
end;
N = size(Qcoef1); N = (N(1)-1)/2;
COUP = zeros(3,3);
NORM = zeros(1,3);
OMEGA(1) = abs( Qcoef1(N+1,N+1,N+1)); Qcoef1(N+1,N+1,N+1)=0; % extracting autonomous frequencies
OMEGA(2) = abs( Qcoef2(N+1,N+1,N+1)); Qcoef2(N+1,N+1,N+1)=0;
OMEGA(3) = abs( Qcoef3(N+1,N+1,N+1)); Qcoef3(N+1,N+1,N+1)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if meth<3; % choosing method
[COUP(1,2), COUP(1,3), COUP(1,1), ~, ~, NORM(1)] = co_3to2(Qcoef1,N,thresh);
else
[~, ~, ~, COUP(1,2), COUP(1,3), NORM(1)] = co_3to2(Qcoef1,N,thresh);
end
if meth==2 || meth==4;
COUP(1,2) = COUP(1,2) / OMEGA(1);
COUP(1,3) = COUP(1,3) / OMEGA(1);
COUP(1,1) = COUP(1,1) / OMEGA(1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if meth<3;
[COUP(2,3), COUP(2,1), COUP(2,2), ~, ~, NORM(2)] = co_3to2(Qcoef2,N,thresh);
else
[~, ~, ~, COUP(2,3), COUP(2,1), NORM(2)] = co_3to2(Qcoef2,N,thresh);
end
if meth==2 || meth==4;
COUP(2,3) = COUP(2,3) / OMEGA(2);
COUP(2,1) = COUP(2,1) / OMEGA(2);
COUP(2,2) = COUP(2,2) / OMEGA(2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if meth<3;
[COUP(3,1), COUP(3,2), COUP(3,3), ~, ~, NORM(3)] = co_3to2(Qcoef3,N,thresh);
else
[~, ~, ~, COUP(3,1), COUP(3,2), NORM(3)] = co_3to2(Qcoef3,N,thresh);
end
if meth==2 || meth==4;
COUP(3,1) = COUP(3,1) / OMEGA(3);
COUP(3,2) = COUP(3,2) / OMEGA(3);
COUP(3,3) = COUP(3,3) / OMEGA(3);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [C12, C13, C123, PART12, PART13, TOT] = co_3to2(Q, N, thresh) % this function sorts the coefficients and computes the norms
thresh=max(abs(Q(:)))*thresh/100;
ind= abs(Q)