%%Updated by Erin Boykin, 2011
%%Adapted from code by Rosenblum and Pikovsky:
%%http://www.agnld.uni-potsdam.de/~mros/damoco.html
%————————————————————————————————-
%Inputs: y – time series for analysis (rows are time points, columns are
% channels)
% This function is for trivariate series with *3* channels
% Fsamp – sampling frequency of signal
%
%Outputs: results contains the following [c12 c13 c21 c23 c31 c32]
% where cij is pdm from i to j
%————————————————————————–
%This function assumes that the signals in y are narrow-band enough that
%they have a well-defined phase. Since this project deals with neuronal
%oscillators which are narrow-band, this assumption is upheld.
function [results] = PDM_3(y,Fsamp)
if size(y,2) ~= 3
error(‘Must have trivarate time series’);
end
npt=length(y);
ncut=1; % number of points to be thrown away in the beginning
% and int the end of signals
npforw = 100; % this is the only parameter for dirc routine
% this parameter corresponds to tau in the paper
%%%% (shift in points to determine delta phi)
% reasonable choice is npforw =
%%%%%%%%%%%%%%%%%%%
%%%%%% Hilbert transform and computation of phases
%%%%%% Phase plots for check
ht=hilbert(y(:,1)); phi1=angle(ht); phi1=phi1(ncut:npt-ncut);
ht=hilbert(y(:,2)); phi2=angle(ht); phi2=phi2(ncut:npt-ncut);
ht=hilbert(y(:,3)); phi3=angle(ht); phi3=phi3(ncut:npt-ncut);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
phi1=unwrap(phi1); phi2=unwrap(phi2); phi3=unwrap(phi3);% unwrap phases
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[results]=dirc(phi1,phi2,phi3,npforw); % inputs are unwrapped phases
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [results]=dirc(phi1,phi2,phi3,npforw) %%% inputs are unwrapped phases
pi2=pi+pi;
npt=length(phi1)-npforw;
dphi1=phi1(npforw+1:npforw+npt);
dphi2=phi2(npforw+1:npforw+npt);
dphi3=phi3(npforw+1:npforw+npt);
phi1=phi1(1:npt); phi2=phi2(1:npt); phi3=phi3(1:npt);
dphi1=dphi1-phi1; dphi2=dphi2-phi2; dphi3=dphi3-phi3;
phi1=mod(phi1,pi2); phi2=mod(phi2,pi2); phi3=mod(phi3,pi2);
% Design coefficient matrix for each channel
X1 = [ones(size(phi1)) sin(phi1) cos(phi1) sin(phi2) cos(phi2) sin(phi3) cos(phi3)...
sin(phi1-phi2) cos(phi1-phi2) sin(phi1-phi3) cos(phi1-phi3) sin(phi1+phi2) cos(phi1+phi2) sin(phi1+phi3) cos(phi1+phi3)...
sin(2*phi1) cos(2*phi1) sin(2*phi2) cos(2*phi2) sin(2*phi3) cos(2*phi3)...
sin(3*phi1) cos(3*phi1) sin(3*phi2) cos(3*phi2) sin(3*phi3) cos(3*phi3)];
X2 = [ones(size(phi2)) sin(phi2) cos(phi2) sin(phi1) cos(phi1) sin(phi3) cos(phi3)...
sin(phi2-phi1) cos(phi2-phi1) sin(phi2-phi3) cos(phi2-phi3) sin(phi1+phi2) cos(phi1+phi2) sin(phi2+phi3) cos(phi2+phi3)...
sin(2*phi2) cos(2*phi2) sin(2*phi1) cos(2*phi1) sin(2*phi3) cos(2*phi3)...
sin(3*phi2) cos(3*phi2) sin(3*phi1) cos(3*phi1) sin(3*phi3) cos(3*phi3)];
X3 = [ones(size(phi3)) sin(phi3) cos(phi3) sin(phi1) cos(phi1) sin(phi2) cos(phi2)...
sin(phi3-phi1) cos(phi3-phi1) sin(phi3-phi2) cos(phi3-phi2) sin(phi1+phi3) cos(phi1+phi3) sin(phi2+phi3) cos(phi2+phi3)...
sin(2*phi3) cos(2*phi3) sin(2*phi1) cos(2*phi1) sin(2*phi2) cos(2*phi2)...
sin(3*phi3) cos(3*phi3) sin(3*phi1) cos(3*phi1) sin(3*phi2) cos(3*phi2)];
% Delta\phi_1
a1=X1\dphi1;
c21=sqrt(a1(4)^2 + a1(5)^2 + a1(8)^2 + a1(9)^2 + a1(12)^2 + a1(13)^2 + 4*a1(18)^2 + 4*a1(19)^2 + 9*a1(24)^2 + 9*a1(25)^2);
c31=sqrt(a1(6)^2 + a1(7)^2 + a1(10)^2 + a1(11)^2 + a1(14)^2 + a1(15)^2 + 4*a1(20)^2 + 4*a1(21)^2 + 9*a1(26)^2 + 9*a1(27)^2);
a2=X2\dphi2;
c12=sqrt(a2(4)^2 + a2(5)^2 + a2(8)^2 + a2(9)^2 + a2(12)^2 + a2(13)^2 + 4*a2(18)^2 + 4*a2(19)^2 + 9*a2(24)^2 + 9*a2(25)^2);
c32=sqrt(a2(6)^2 + a2(7)^2 + a2(10)^2 + a2(11)^2 + a2(14)^2 + a2(15)^2 + 4*a2(20)^2 + 4*a2(21)^2 + 9*a2(26)^2 + 9*a2(27)^2);
a3=X3\dphi3;
c13=sqrt(a3(4)^2 + a3(5)^2 + a3(8)^2 + a3(9)^2 + a3(12)^2 + a3(13)^2 + 4*a3(18)^2 + 4*a3(19)^2 + 9*a3(24)^2 + 9*a3(25)^2);
c23=sqrt(a3(6)^2 + a3(7)^2 + a3(10)^2 + a3(11)^2 + a3(14)^2 + a3(15)^2 + 4*a3(20)^2 + 4*a3(21)^2 + 9*a3(26)^2 + 9*a3(27)^2);
if (c21<0.01) && (c31<0.01) && (c12<0.01) && (c32<0.01) && (c13<0.01) && (c23<0.01)
disp('WARNING: coefficients c1, c2, c3 are too small!');
disp('Be careful: probably the signals are not correlated');
end;
results=[c12 c13 c21 c23 c31 c32];
