A python port of a method described by Adhikari et al. 2010 to estimate the directionality and the lag between the eeg recordings from two brain areas.

amp_crosscorr

# -*- coding: utf-8 -*-
"""
@author: Daniel Müller-Komorowska
"""

import numpy as np
import scipy.signal as signal

def amp_crosscorr(eeg1, eeg2, samp_freq, low_freq, high_freq):
    """
    """
    
    if len(eeg1) != len(eeg2):
        raise ValueError("eeg1 and eeg2 are not equal in length")
    if len(eeg1.shape) != 1:
        raise ValueError("eeg1 must be one dimensional array")
    if len(eeg2.shape) != 1:
        raise ValueError("eeg2 must be one dimensional array")
    
    order = round(samp_freq)
    if order % 2 != 0:
        order = order - 1
    
    nyquist = np.floor(samp_freq/2)
    
    # signal.firwin() is the scipy equivalent of matlabs fir1()
    # Original: MyFilt=fir1(order,[low_freq high_freq]/Nyquist)
    # We have to change some defaults of firwin to get the same behavior
    cutoff_freqs = np.array([low_freq, high_freq])
    filt_coeff = signal.firwin(order, cutoff_freqs / nyquist)
    
    # Publication uses Filter0()
    # Original: filtered1 = Filter0(MyFilt,eeg1)
    # Found no matlab function Filter0() so I pretend they used filtfilt()
    # scipy.signal has filtfilt() implemented
    filtered1 = signal.filtfilt(filt_coeff, eeg1)
    filtered2 = signal.filtfilt(filt_coeff, eeg2)
    
    # The hilbert transform splits the signal into a real and an imaginary part
    # The imaginary part is the hilbert transform.
    # Original: filt_hilb1 = hilbert(filtered1)
    # scipy.signal.hilbert should behave similarly.
    filt_hilb1 = signal.hilbert(filtered1)
    filt_hilb2 = signal.hilbert(filtered2)
    
    # In MATLAB the abs() function returns the imaginary.
    # We use the numpy.absolute function for this.
    amp1 = np.absolute(filt_hilb1)
    amp2 = np.absolute(filt_hilb2)
    
    # Now we simply subtract the mean from amp1 and amp2
    amp1 = amp1 - amp1.mean()
    amp2 = amp2 - amp2.mean()
    
    # Original: [crosscorr,lags]=xcorr(amp1, amp2,round(samp_freq/10),'coeff');
    # We use scipy.signal.correlate(). Unlike MATLABs xcorr, we do not get lag
    # We need to calculate the lag for ourselves
    # I am not sure if method = 'direct' is the best choice
    ycorr = signal.correlate(amp1, amp2, method = "direct")
    
    # Now we calculate the lag in ms
    xcorr = np.arange(ycorr.size)
    
    lags = xcorr - (amp1.size-1)
    
    lags = -lags * (1/samp_freq)
    
    idx_max_corr = ycorr.argmax()
    
    max_corr_lag = lags[idx_max_corr]