Created
December 29, 2019 22:21
-
-
Save danielmk/7eafbfad7700eaf2784c6882e13c1861 to your computer and use it in GitHub Desktop.
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# -*- 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] |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment