Module with functions for spectral Granger causality

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

Typical usage is as follows: from bsmart import pwcausalr F,pp,cohe,Fx2y,Fy2x,Fxy=pwcausalr(x,ntrls,npts,p,fs,freq);


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)


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: 2019jun17 by Cliff Kerr (


timefreq(x[, fs])




armorf(x, ntrls, npts, p)

spectrum_AR(A, Z, M, f, fs)

pwcausalr(x, Nr, Nl, porder, fs[, freq])

granger(vec1, vec2[, order, rate, maxfreq])

GRANGER, fs=200)[source]


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: 2019jun17[source]


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: 2019jun17, ntrls, npts, p)[source], Z, M, f, fs)[source], Nr, Nl, porder, fs, freq=0)[source], vec2, order=10, rate=200, maxfreq=0)[source]


Provide a simple way of calculating the key quantities.




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: 2019jun17