function filtdir
% This function filters two signals, estimates phases and computes
% two directionality indices
%
% One algorithm can be commented out because it is rather time consuming
% (see below)
% For description of algorithms see:
%  www.agnld.uni-potsdam.de/~mros/pdf/PRE41909.pdf

%%%%%% load the signals %%%%%%%%%%%%
%%%%%% you have to adjust this part to input your data 
load input_file;
t=M(:,1); x=M(:,2); y=M(:,4);
clear M;
%% t is time, x is signal 1, y is signal 2
%%%%%%%%%%%%%%
Fsamp=500; fn=Fsamp/2.; % sampling and Nyquist frequency
pi2=pi+pi; npt=length(x);
ncut=100;   % 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 = <oscillation period>

%%%%% filter design %%%%%%%%%
%%%%% if your data do not require filtering then just 
%%%%% skip this and copy x to xf, y to yf
%%%%%
%%%%% if your data requires more complex preprocessing then 
%%%%% you may insert it here
%%%%%
%%%%%  !!!!!!!! IMPORTANT: often you really need preprocessing to 
%%%%%  obtain well-defined (or, at least, reasonable :-)) phase!
%%%%%  Otherwise the following computations may be completely 
%%%%%  senseless! 
%%%%%
Wn=[4.5 5.5]/fn;     % determines the bandpass
b = fir1(1024,Wn);   % determines the filter length

xf=filtfilt(b,1,x);  % filtering data
yf=filtfilt(b,1,y);
%%%%%%%%%%%%%%%%%%%
%%%%%% Hilbert transform and computation of phases
%%%%%% Phase plots for check
figure(1);
ht=hilbert(xf); phi1=angle(ht); phi1=phi1(ncut:npt-ncut);
subplot(1,3,1); plot(ht);
xlabel('x'), ylabel('x_H'), axis square;
ht=hilbert(yf); phi2=angle(ht); phi2=phi2(ncut:npt-ncut);
subplot(1,3,2); plot(ht);
xlabel('y'), ylabel('y_H'), axis square;
subplot(1,3,3); plot(phi1,phi2,'r.','MarkerSize',1);
xlabel('\phi_1'), ylabel('\phi_2'), axis square; axis([-pi pi -pi pi]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
phi1=unwrap(phi1);       phi2=unwrap(phi2);% unwrap phases
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);
subplot(2,1,1);   plot(t,xf,t,yf,'r'); % plot filtered data
title('filtered data');
subplot(2,1,2);   plot(t(ncut:npt-ncut),(phi1-phi2)/pi2);
title('phase difference');

%%%%%% Check synchronization index
synchroindex=(mean(cos(phi1-phi2)))^2 + (mean(sin(phi1-phi2)))^2;
if synchroindex > 0.6
  disp('WARNING: synchronization index is rather high!!!');
  synchroindex
  disp('The results might be not significant!');
  disp('Comparison of different algorithms is recommended');
end
%%%%%%%% directionality indices
%%%%%%%% EMA %%%%%%%%%%%%%%%%%
results=dirc(phi1,phi2,npforw); % inputs are unwrapped phases
disp('       c1       c2       d12');
disp(results);
%%%%%%% IPA : rather time consuming
% comment out next 3 lines if you do not want second method
results=prdirc(phi1,phi2);      % inputs are unwrapped phases
disp('       c1       c2       r12');
disp(results);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% EPA algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [results]=dirc(phi1,phi2,npforw) %%% inputs are unwrapped phases
% Version for LFP-EMG data
pi2=pi+pi;

npt=length(phi1)-npforw;
dphi1=phi1(npforw+1:npforw+npt);
dphi2=phi2(npforw+1:npforw+npt);
phi1=phi1(1:npt);   phi2=phi2(1:npt);
dphi1=dphi1-phi1;   dphi2=dphi2-phi2;
phi1=mod(phi1,pi2); phi2=mod(phi2,pi2);

% design matrix is common for both fittings
X = [ones(size(phi1)) sin(phi1) cos(phi1) sin(phi2) cos(phi2) ...
     sin(phi1-phi2) cos(phi1-phi2) sin(phi1+phi2) cos(phi1+phi2)...
     sin(2*phi1) cos(2*phi1) sin(2*phi2) cos(2*phi2)...
     sin(3*phi1) cos(3*phi1) sin(3*phi2) cos(3*phi2)];
%  Delta\phi_1
a=X\dphi1;  c1=sqrt(dzdy2(a));
a=X\dphi2;  c2=sqrt(dzdx2(a));

if (c1<0.01) & (c2<0.01)
 disp('WARNING: coefficients c1 and c2 are too small!');
 disp('Be careful: probably the signals are not correlated');
end;
d12=(c2-c1)/(c1+c2);
results=[c1 c2 d12];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r=dzdx2(a)
 r= 4*a(10)^2 + 4*a(11)^2 + 9*a(14)^2 + 9*a(15)^2 + a(2)^2 ...
    + a(3)^2 + a(6)^2 + a(7)^2 + a(8)^2 + a(9)^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r=dzdy2(a)
 r= 4*a(12)^2 + 4*a(13)^2 + 9*a(16)^2 + 9*a(17)^2 + a(4)^2 ...
    + a(5)^2 + a(6)^2 + a(7)^2 + a(8)^2 + a(9)^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% IPA algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [results]=prdirc(phi1,phi2)       %%% inputs are unwrapped phases
% Dependence of instantaneous period on the phase of the second oscillator

[rt1,npt1]=instper(phi1);  % compute series of instantaneous periods
[rt2,npt2]=instper(phi2);

npt=min(npt1,npt2);

phi1=phi1(1:npt); phi2=phi2(1:npt); rt1=rt1(1:npt); rt2=rt2(1:npt);
% design matrix for fitting
X = [ones(size(phi1)) sin(phi1) cos(phi1) sin(phi2) cos(phi2) ...
     sin(phi1-phi2) cos(phi1-phi2) sin(phi1+phi2) cos(phi1+phi2)...
     sin(2*phi1) cos(2*phi1) sin(2*phi2) cos(2*phi2)...
     sin(3*phi1) cos(3*phi1) sin(3*phi2) cos(3*phi2)];
% approximate inst periods for 1st signal
a=X\rt1';   c1=sqrt(dzdy2(a));
% approximate inst periods for 2nd signal
a=X\rt2';   c2=sqrt(dzdx2(a));

results=[c1 c2 (c2-c1)/(c1+c2)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rt,nptnew]=instper(phi) % input: unwrapped phases

pi2=pi+pi;
npspl=5;  % how many points to take for spline approximation

npt=length(phi); tim=1:npt;
finphi=phi(npt)-pi2;
nptnew=npt;
while phi(nptnew) > finphi
  nptnew=nptnew-1;
end;

ipi2=1;
for i=1:nptnew
 nextphi=phi(i)+pi2;
 while phi(ipi2) < nextphi
   ipi2=ipi2+1;
 end;
%%%%% remove comments if you want spline interpolation %%%%%%%%%
%bef=ipi2-5; if bef < 1 bef=1; end;
%aft=ipi2+5; if aft npt aft=npt; end;
%tspl=tim(bef:aft);
%phispl=phi(bef:aft);
%tintpl=spline(phispl,tspl,nextphi);
%rt(i)=tintpl-i;
%
%%%%%  next lines are for faster version: linear interpolation
 timfrac=(nextphi-phi(ipi2-1))/(phi(ipi2)-phi(ipi2-1));
 rt(i)=ipi2-1-i+timfrac;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%