""" This file contains all the function definitions necessary for running spectral Granger causality. It is based on Mingzhou Ding's Matlab code package BSMART, available from www.brain-smart.org. Typical usage is as follows: from bsmart import pwcausalr F,pp,cohe,Fx2y,Fy2x,Fxy=pwcausalr(x,ntrls,npts,p,fs,freq); Outputs: F is the frequency vector for the remaining quantities pp is the spectral power cohe is the coherence Fx2y is the causality of channel X to channel Y Fy2x is the causality of channel Y to channel X Fxy is the "instantaneous" causality (cohe-Fx2y-Fy2x I think) Inputs: x is the data for at least two channels, e.g. a 2x8000 array consisting of two LFP time series ntrls is the number of trials (whatever that means -- just leave it at 1) npts is the number of points in the data (in this example, 8000) p is the order of the polynomial fit (e.g. 10 for a smooth fit, 20 for a less smooth fit) fs is the sampling rate (e.g. 200 Hz) freq is the maximum frequency to calculate (e.g. fs/2=100, which will return 0:100 Hz) The other two functions (armorf and spectrum_AR) can also be called directly, but more typically they are used by pwcausalr in intermediate calculations. Note that the sampling rate of the returned quantities is calculated as fs/2. To calculate the power spectrum powspec of a single time series x over the frequency range 0:freq, use the following (NB: now accessible via "from spectrum import ar") from bsmart import armorf, spectrum_AR [A,Z,tmp]=armorf(x,ntrls,npts,p) # Calculate autoregressive fit for i in range(freq+1): # Loop over frequencies [S,H]=spectrum_AR(A,Z,p,i,fs) # Calculate spectrum powspec[i]=abs(S**2) # Calculate and store power In either case (pwcausalr or spectrum_AR), the smoothness of the spectra is determined by the polynomial order p. Larger values of p give less-smooth spectra. Version: 2012dec10 by Cliff Kerr (cliffk@neurosim.downstate.edu) """ # ARMORF -- AR parameter estimation via LWR method modified by Morf. # # X is a matrix whose every row is one variable's time series # ntrls is the number of realizations, npts is the length of every realization # If the time series are stationary long, just let ntrls=1, npts=length(x) # # A = ARMORF(X,NR,NL,ORDER) returns the polynomial coefficients A corresponding to # the AR model estimate of matrix X using Morf's method. # ORDER is the order of the AR model. # # [A,E] = ARMORF(...) returns the final prediction error E (the variance # estimate of the white noise input to the AR model). # # [A,E,K] = ARMORF(...) returns the vector K of reflection coefficients (parcor coefficients). # # Ref: M. Morf, etal, Recursive Multichannel Maximum Entropy Spectral Estimation, # IEEE trans. GeoSci. Elec., 1978, Vol.GE-16, No.2, pp85-94. # S. Haykin, Nonlinear Methods of Spectral Analysis, 2nd Ed. # Springer-Verlag, 1983, Chapter 2 def timefreq(x,fs=200): """ TIMEFREQ This function takes the time series and the sampling rate and calculates the total number of points, the maximum frequency, the minimum (or change in) frequency, and the vector of frequency points F. Version: 2011may04 """ from numpy import size, arange, append maxfreq=float(fs)/2.0 # Maximum frequency minfreq=float(fs)/float(size(x,0)) # Minimum and delta frequency -- simply the inverse of the length of the recording in seconds F=arange(minfreq,maxfreq+minfreq,minfreq) # Create frequencies evenly spaced from 0:minfreq:maxfreq F=append(0,F) # Add zero-frequency component return F def ckchol(M): """ CKCHOL This function computes the Cholesky decomposition of the matrix if it's positive-definite; else it returns the identity matrix. It was written to handle the "matrix must be positive definite" error in linalg.cholesky. Version: 2011may03 """ from numpy import linalg, matrix, eye, size try: # First, try the Cholesky decomposition output=linalg.cholesky(M) except: # If not, just return garbage print 'WARNING: Cholesky failed, so returning (invalid) identity matrix!' output=matrix(eye(size(M,0))) return output def armorf(x,ntrls,npts,p): from scipy import shape, array, matrix, zeros, concatenate, eye, dstack from numpy import linalg # for inverse and Cholesky factorization; inv=linalg.inv; # Make name consistent with Matlab # Initialization x=matrix(x) [L,N]=shape(x); # L is the number of channels, N is the npts*ntrls pf=pb=pfb=ap=bp=En=matrix(zeros((L,L,1))); # covariance matrix at 0, # calculate the covariance matrix? for i in range(ntrls): En=En+x[:,i*npts:(i+1)*npts]*x[:,i*npts:(i+1)*npts].H; ap=ap+x[:,i*npts+1:(i+1)*npts]*x[:,i*npts+1:(i+1)*npts].H; bp=bp+x[:,i*npts:(i+1)*npts-1]*x[:,i*npts:(i+1)*npts-1].H; ap = inv((ckchol(ap/ntrls*(npts-1)).T).H); bp = inv((ckchol(bp/ntrls*(npts-1)).T).H); for i in range(ntrls): efp = ap*x[:,i*npts+1:(i+1)*npts]; ebp = bp*x[:,i*npts:(i+1)*npts-1]; pf = pf + efp*efp.H; pb = pb + ebp*ebp.H; pfb = pfb + efp*ebp.H; En = (ckchol(En/N).T).H; # Covariance of the noise # Initial output variables tmp=[] for i in range(L): tmp.append([]) # In Matlab, coeff=[], and anything can be appended to that. coeff = matrix(tmp);# Coefficient matrices of the AR model kr = matrix(tmp); # reflection coefficients aparr=array(ap) # Convert AP matrix to an array, so it can be dstacked bparr=array(bp) for m in range(p): # Calculate the next order reflection (parcor) coefficient ck = inv((ckchol(pf).T).H)*pfb*inv(ckchol(pb).T); kr=concatenate((kr,ck),1); # Update the forward and backward prediction errors ef = eye(L)- ck*ck.H; eb = eye(L)- ck.H*ck; # Update the prediction error En = En*(ckchol(ef).T).H; # Update the coefficients of the forward and backward prediction errors Z=zeros((L,L)) # Make it easier to define this aparr=dstack((aparr,Z)) bparr=dstack((bparr,Z)) pf = pb = pfb = Z # Do some variable juggling to handle Python's array/matrix limitations a=b=zeros((L,L,0)) for i in range(m+2): tmpap1=matrix(aparr[:,:,i]) # Need to convert back to matrix to perform operations tmpbp1=matrix(bparr[:,:,i]) tmpap2=matrix(aparr[:,:,m+1-i]) tmpbp2=matrix(bparr[:,:,m+1-i]) tmpa = inv((ckchol(ef).T).H)*(tmpap1-ck*tmpbp2); tmpb = inv((ckchol(eb).T).H)*(tmpbp1-ck.H*tmpap2); a=dstack((a,array(tmpa))) b=dstack((b,array(tmpb))) for k in range(ntrls): efp = zeros((L,npts-m-2)); ebp = zeros((L,npts-m-2)); for i in range(m+2): k1=m+2-i+k*npts; k2=npts-i+k*npts; efp = efp+matrix(a[:,:,i])*matrix(x[:,k1:k2]); ebp = ebp+matrix(b[:,:,m+1-i])*matrix(x[:,k1-1:k2-1]); pf = pf + efp*efp.H; pb = pb + ebp*ebp.H; pfb = pfb + efp*ebp.H; aparr = a; bparr = b; for j in range(p): coeff = concatenate((coeff,inv(matrix(a[:,:,0]))*matrix(a[:,:,j+1])),1); return coeff, En*En.H, kr #Port of spectrum_AR.m # Version: 2010jan18 def spectrum_AR(A,Z,M,f,fs): # Get the spectrum in one specific frequency-f from scipy import eye, size, exp, pi from numpy import linalg; inv=linalg.inv N = size(Z,0); H = eye(N,N); # identity matrix for m in range(M): H = H + A[:,m*N:(m+1)*N]*exp(-1j*(m+1)*2*pi*f/fs); # Multiply f in the exponent by sampling interval (=1/fs). See Richard Shiavi H = inv(H); S = H*Z*H.H/fs; return S,H # Using Geweke's method to compute the causality between any two channels # # x is a two dimentional matrix whose each row is one variable's time series # Nr is the number of realizations, # Nl is the length of every realization # If the time series have one ralization and are stationary long, just let Nr=1, Nl=length(x) # porder is the order of AR model # fs is sampling frequency # freq is a vector of frequencies of interest, usually freq=0:fs/2 # CK: WRONG!! freq must be a scalar, else the for loop doesn't work. # # Fx2y is the causality measure from x to y # Fy2x is causality from y to x # Fxy is instantaneous causality between x and y # the order of Fx2y/Fy2x is 1 to 2:L, 2 to 3:L,....,L-1 to L. That is, # 1st column: 1&2; 2nd: 1&3; ...; (L-1)th: 1&L; ...; (L(L-1))th: (L-1)&L. # revised Jan. 2006 by Yonghong Chen # Note: remove the ensemble mean before using this code def pwcausalr(x,Nr,Nl,porder,fs,freq=0): # Note: freq determines whether the frequency points are calculated or chosen from pylab import size, shape, real, log, conj, zeros, arange, array from numpy import linalg; det=linalg.det import numpy as np # Just for "sum"; can't remember what's wrong with pylab's sum [L,N] = shape(x); #L is the number of channels, N is the total points in every channel if freq==0: F=timefreq(x[0,:],fs) # Define the frequency points else: F=array(range(0,freq+1)) # Or just pick them npts=size(F,0) # Initialize arrays maxindex=np.sum(arange(1,L)) pp=zeros((L,npts)) # Had these all defined on one line, and stupidly they STAY linked!! cohe=zeros((maxindex,npts)) Fy2x=zeros((maxindex,npts)) Fx2y=zeros((maxindex,npts)) Fxy=zeros((maxindex,npts)) index = 0; for i in range(1,L): for j in range(i+1,L+1): y=zeros((2,N)) # Initialize y index = index + 1; y[0,:] = x[i-1,:]; y[1,:] = x[j-1,:]; A2,Z2,tmp = armorf(y,Nr,Nl,porder); #fitting a model on every possible pair eyx = Z2[1,1] - Z2[0,1]**2/Z2[0,0]; #corrected covariance exy = Z2[0,0] - Z2[1,0]**2/Z2[1,1]; f_ind = 0; for f in F: f_ind = f_ind + 1; S2,H2 = spectrum_AR(A2,Z2,porder,f,fs); pp[i-1,f_ind-1] = abs(S2[0,0]*2); # revised if (i==L-1) & (j==L): pp[j-1,f_ind-1] = abs(S2[1,1]*2); # revised cohe[index-1,f_ind-1] = real(abs(S2[0,1])**2 / S2[0,0]/S2[1,1]); Fy2x[index-1,f_ind-1] = log(abs(S2[0,0])/abs(S2[0,0]-(H2[0,1]*eyx*conj(H2[0,1]))/fs)); #Geweke's original measure Fx2y[index-1,f_ind-1] = log(abs(S2[1,1])/abs(S2[1,1]-(H2[1,0]*exy*conj(H2[1,0]))/fs)); Fxy[index-1,f_ind-1] = log(abs(S2[0,0]-(H2[0,1]*eyx*conj(H2[0,1]))/fs)*abs(S2[1,1]-(H2[1,0]*exy*conj(H2[1,0]))/fs)/abs(det(S2))); return F,pp,cohe,Fx2y,Fy2x,Fxy def granger(vec1,vec2,order=10,rate=200,maxfreq=0): """ GRANGER Provide a simple way of calculating the key quantities. Usage: F,pp,cohe,Fx2y,Fy2x,Fxy=granger(vec1,vec2,order,rate,maxfreq) where: F is a 1xN vector of frequencies pp is a 2xN array of power spectra cohe is the coherence between vec1 and vec2 Fx2y is the causality from vec1->vec2 Fy2x is the causality from vec2->vec1 Fxy is non-directional causality (cohe-Fx2y-Fy2x) vec1 is a time series of length N vec2 is another time series of length N rate is the sampling rate, in Hz maxfreq is the maximum frequency to be returned, in Hz Version: 2011jul18 """ from bsmart import timefreq, pwcausalr from scipy import array, size if maxfreq==0: F=timefreq(vec1,rate) # Define the frequency points else: F=array(range(0,maxfreq+1)) # Or just pick them npts=size(F,0) data=array([vec1,vec2]) F,pp,cohe,Fx2y,Fy2x,Fxy=pwcausalr(data,1,npts,order,rate,maxfreq) return F,pp[0,:],cohe[0,:],Fx2y[0,:],Fy2x[0,:],Fxy[0,:]