%%Written by Erin Boykin, 2011
%————————————————————————————————-
%Inputs: Nt – number of time points in signal
% NC – number of channels or time series
% w – intercept vector from ARfit
% A – coefficient matrix from ARfit
% Sig – noise covariance matrix from ARfit
% freq – vector of discrete frequencies at which gc is calculated
%
%Outputs: gc – gc values in form frequency x channel i x channel j
% where causality is from channel i to channel j
%————————————————————————–
%Note that this function requires one to first fit an AR model to the time
%series. It is assumed that this is done using the ARfit package available
%here: http://www.gps.caltech.edu/~tapio/arfit/
%
function [gc] = pairwiseGC(Nt,Nc,w,A,Sig,freq)
gc = zeros(length(freq),Nc,Nc);
for xx = 1:Nc
for yy = xx+1:Nc
if (yy ~= xx)
[S,H] = sig2spect_arfitin(Nt,Nc,w,A,Sig,freq); %Signals to AR auto&cross-spectra
for ff = 1:length(freq)
%%If Geweke’s transform is done earlier then:
%causality(ff,yy,xx) = log(abs(S(1,1,ff))/abs((H(1,1,ff)*Sig(1,1)*conj(H(1,1,ff)))));
%causality(ff,xx,yy) = log(abs(S(2,2,ff))/abs((H(2,2,ff)*Sig(2,2)*conj(H(2,2,ff)))));
%%If Geweke’s Transform has not been done:
denom = Sig(2,2)-Sig(1,2)^2/Sig(1,1);
gc(ff,yy,xx) = log(abs(S(1,1,ff))/(abs(S(1,1,ff))-denom*(abs(H(1,2,ff)))^2));
denom = Sig(1,1)-Sig(2,1)^2/Sig(2,2);
gc(ff,xx,yy) = log(abs(S(2,2,ff))/(abs(S(2,2,ff))-denom*(abs(H(2,1,ff)))^2));
end
end
end
end
function [S,H]= sig2spect_arfitin(Nt,m,w,B,Sig,freq)
popt = size(B,2) / m;
A0 = eye(m);
Bf = zeros(m,m,length(freq));
for g=1:length(freq)
for ch=1:m
Bj = B(:,[ch:m:m*popt]);
for j = 1:popt
Bf(:,ch,g) = Bf(:,ch,g) + Bj(:,j).*exp(-2*pi*1i*j*freq(g));
end
end
Bf(:,:,g) = A0 – Bf(:,:,g);
end
%if m == 2 %Do Geweke’s transform
% Bf(2,1,:) = Bf(2,1,:) – Sig(1,2)/Sig(1,1)*Bf(1,1,:);
% Bf(2,2,:) = Bf(2,2,:) – Sig(1,2)/Sig(1,1)*Bf(1,2,:);
% Sig(2,2) = Sig(2,2) – Sig(1,2)^2/Sig(1,1);
%end
for g = 1:length(freq)
H(:,:,g) = inv(Bf(:,:,g));
S(:,:,g) = H(:,:,g)*Sig*H(:,:,g)’;
end
