Settings (These should be changed to fit your system)



In [ ]:

    
data_folder = "/home/mje/Projects/agency_connectivity/data/test/"  #Set the path to the data
scripts_folder  = "/home/mje/Projects/agency_connectivity"  #Set the path to the scripts



In [ ]:

    
import mne
import matplotlib.pyplot as plt 
import numpy as np
import seaborn as sns
import os

%matplotlib qt
# change "qt" to "inline" for the figures to be place in the notebook



In [ ]:

    
os.chdir(scripts_folder)

Home made functions



In [ ]:

    
from preprocessing import convert_bdf2fif, filter_raw, compute_ica, epoch_data, hilbert_process
from tf_functions import morlet_analysis, single_trial_tf, calc_spatial_resolution, calc_wavelet_duration
from connectivity_functions import ISPC_over_trials

Preprocessing data

Convert from bdf to fiff file (the mne format)



In [ ]:

    
convert_bdf2fif("P2", data_folder)

Filter data

The filter_raw function apply a bandpass from 1Hz to 98Hz



In [ ]:

    
filter_raw("P2", data_folder)

ICA

Compute ICA and automactically remove the component with the highest correlation with the EOG channel (EXGF4)



In [ ]:

    
compute_ica("P2", data_folder)

Epocing



In [ ]:

    
epoch_data("P2", data_folder)



In [ ]:

    
epochs = mne.read_epochs(data_folder + "P2_ds_bp_ica-epo.fif")

Evoked



In [ ]:

    
evoked_vol.plot_joint(title="voluntary", topomap_args={"cmap": "viridis"})



In [ ]:

    
evoked_inv.plot_joint(title="involuntary", topomap_args={"cmap": "viridis"})

Time-frequency analysis



In [ ]:

    
# load epochs
epochs = mne.read_epochs(data_folder + "P2_ds_bp_ica-epo.fif")



In [ ]:

    
tfr_vol, itc_vol = morlet_analysis(epochs["voluntary"], n_cycles=4)
tfr_invol, itc_invol = morlet_analysis(epochs["involuntary"], n_cycles=4)



In [ ]:

    
print("vol max: %s \ninvol max: %s") % (tfr_vol.data.max(), tfr_invol.data.max())
print("vol min: %s \ninvol min: %s") % (tfr_vol.data.min(), tfr_invol.data.min())

Plot results

Average power



In [ ]:

    
tfr_vol.plot_topo(baseline=(None, -1.8), mode="zscore", cmap="viridis", vmax=30, vmin=-30)



In [ ]:

    
tfr_invol.plot_topo(baseline=(None, -1.8), mode="zscore", cmap="viridis", vmax=30, vmin=-30)

Intertrial coherence



In [ ]:

    
itc_vol.plot_topo(cmap="viridis", vmax=1, vmin=0)



In [ ]:

    
itc_invol.plot_topo(cmap="viridis", vmax=1, vmin=0)

Make single trial Morlet transform of epochs

The "frequencies", "n_cycles", and "times" are not used but for convinience when plotting etc.



In [ ]:

    
frequencies = np.arange(6., 30., 1.)
n_cycles = 5.
times = epochs.times

tfr_vol = single_trial_tf(epochs["voluntary"], n_cycles=4)
tfr_invol = single_trial_tf(epochs["involuntary"], n_cycles=4)

# convert to numpy array
tfr_vol = np.asarray(tfr_vol)
tfr_invol = np.asarray(tfr_invol)

Calculate spectral and temporal resolution of the Morlet wavelets

From http://www.fieldtriptoolbox.org/tutorial/timefrequencyanalysis (Time-frequency analysis IV -- Morlet Wavelets). We get the spectral bandwidth is:

$$\frac{\text{Frequency}}{\text{number of cycles}}\times2 = \text{Spectral resolution Hz}$$



In [ ]:

    
frequencies = np.arange(6., 31., 1.)



In [ ]:

    
spectral_resolution_4 = calc_spatial_resolution(frequencies, n_cycles=4)
spectral_resolution_5 = calc_spatial_resolution(frequencies, n_cycles=5)

plt.figure()
plt.plot(frequencies, spectral_resolution_4, label="4 cycles")
plt.plot(frequencies, spectral_resolution_5, label="5 cycles")
plt.ylabel("Spectral bandwidth (Hz)")
plt.xlabel("Frequencies")
plt.legend()
plt.show()

plt.figure()
plt.plot(frequencies, spectral_resolution_4 - spectral_resolution_5)
plt.ylabel("Spectral bandwidth Hz ")
plt.xlabel("Frequencies")
plt.title("Difference: 4 cycles - 5 cycles")
plt.show()

From http://www.fieldtriptoolbox.org/tutorial/timefrequencyanalysis (Time-frequency analysis IV -- Morlet Wavelets). We get the wavelet duration is:

$$\Bigg(\frac{\frac{\text{Number of cycles}}{\text{Frequency}}}{\pi}\Bigg)\times 1000 = \text{Wavelet duration (in miliseconds)}$$



In [ ]:

    
Wavelet_duration_4 = calc_wavelet_duration(frequencies, n_cycles=4)
Wavelet_duration_5 = calc_wavelet_duration(frequencies, n_cycles=5)

plt.figure()
plt.plot(frequencies, Wavelet_duration_4, label="4 cycles")
plt.plot(frequencies, Wavelet_duration_5, label="5 cycles")
plt.ylabel("Wavelet duration (in ms)")
plt.xlabel("Frequencies")
plt.legend()
plt.show()

plt.figure()
plt.plot(frequencies, Wavelet_duration_5 - Wavelet_duration_4, label="difference")
plt.ylabel("Wavelet duration (in ms)")
plt.xlabel("Frequencies")
plt.title("Difference: 5 cycles - 4 cycles")
plt.legend()
plt.show()

Connectivity: ISPC

Calculate ISPC over trials



In [ ]:

    
chan_A, chan_B = 37, 47  # change for different channels
freq_idx = 14  # change for another frequciency

ispc_vol = ISPC_over_trials(tfr_vol, freqs=[2, 7], channels=[chan_A, chan_B], faverage=True)
ispc_vol = ISPC_over_trials(tfr_invol, freqs=[2, 7], channels=[chan_A, chan_B], faverage=True)



In [ ]:

    
plt.figure()
plt.plot(times, ispc_vol, 'b', label="Voluntary")
plt.plot(times, ispc_vol, 'r', label="Involuntary")
plt.legend()
plt.title("ISPC over time between %s and %s for freq: %s:%s" % (epochs.ch_names[chan_A],
                                                  epochs.ch_names[chan_B],
                                                  frequencies[2], frequencies[7]))
plt.show()

Calculate ISPC over time



In [ ]:

    
ispc_vol = np.empty(tfr_vol.shape[-1])
ispc_invol = np.empty(tfr_invol.shape[-1])

chan_A, chan_B = 37, 47  # change for different channels
freq_idx = 14  # change for another frequciency

for i in range(len(ispc_vol)):
    ispc_vol[i] = np.abs(np.mean(np.exp(
        1j*(np.angle(tfr_vol[:, chan_A, freq_idx, i]) -
            np.angle(tfr_vol[:, chan_B, freq_idx, i])))))

for i in range(len(ispc_invol)):
    ispc_invol[i] = np.abs(np.mean(np.exp(
        1j*(np.angle(tfr_invol[:, chan_A, freq_idx, i]) -
            np.angle(tfr_invol[:, chan_B, freq_idx, i])))))



In [ ]:

    
plt.figure()
plt.plot(times, ispc_vol, 'b', label="Voluntary")
plt.plot(times, ispc_invol, 'r', label="Involuntary")
plt.legend()
plt.title("ISPC over time between %s and %s for freq: %s:%s" % (epochs.ch_names[chan_A],
                                                  epochs.ch_names[chan_B],
                                                  frequencies[2], frequencies[7]))
plt.show()



In [ ]:



In [ ]:



In [ ]:



In [ ]:



In [ ]:



In [ ]:



In [ ]: