In [1]:

    

import numpy as np

import scipy.io
import scipy.signal

from matplotlib import rcParams, animation
import matplotlib.pyplot as plt
%matplotlib inline

from IPython.display import HTML


    

In [2]:

    

# Global configuration.
FILENAME = "../../data/megicann_train_v2.mat"
DATASET_X = 'train_day2'
FREQ_DOMAINS = [3, 4, 5]
WINDOW_LENGTH = 40


    

In [3]:

    

# This function preprocesses data for connectivity analysis:
# * Original data is in format train[frequency][sample, channel, time]. The function makes different
#   frequencies appear as channels and swap channel and time.
# * Compute hilbert translation for the data to get the phase data
# * Normalize the data so that you can easily compute e.g. PLV, PLI, wPLI, entropy, etc.

def preprocess_x(raw_data):
    # Move different frequencies to channels
    samples = raw_data[0, FREQ_DOMAINS[0]][:]
    for freq_domain in FREQ_DOMAINS[1:]:
        samples = np.concatenate([samples, raw_data[0, freq_domain][:]], axis=1)
        
    # Normalize the data within each sample (per channel) and ensure that type is ok
    for i, sample in enumerate(samples):
        sample_mean = sample.mean(axis=1);
        sample_std = sample.std(axis=1);
        samples[i] = ((sample.transpose() - sample_mean) / sample_std).transpose()

    samples = scipy.signal.hilbert(samples)
    samples = samples / abs(samples)
    
    # Transpose the data to [sample, time, channel] format
    samples = np.transpose(samples, axes=(0, 2, 1))
    
    return samples


    

In [4]:

    

# Load samples and preprocess them
train_dict = scipy.io.loadmat(FILENAME)
dataset = preprocess_x(train_dict[DATASET_X])


    

In [5]:

    

# Compute windowed PLV for the data and return connectivity chart in format [sample, time, channel, channel]
def wplv(samples, window_length):
    plvs = []
    for sample in samples:
        sliding_plv = []
        for i in range(sample.shape[0] - window_length):
            _sample = sample[i:i+window_length]
            plv = np.abs(np.dot(_sample.conj().T, _sample)) / window_length
            np.fill_diagonal(plv, 0)
            sliding_plv.append(plv)
        plvs.append(sliding_plv)
        
    return np.asarray(plvs)


    

In [6]:

    

# Compute connectivity for a single sample
sample = wplv(dataset[0:1], WINDOW_LENGTH)

# Default is ffmpeg - and I couldn't get it working right away. Change writer to avconv.
plt.rcParams['animation.writer'] = 'avconv'

# Create figure and plot the first frame - but don't show it.
fig = plt.figure(figsize=(10,10))
im = plt.imshow(sample[0, 0], animated=True, cmap=plt.get_cmap('hot'))
plt.colorbar()
plt.tight_layout()
plt.close()

# Function to plot a single frame
def updatefig(i):
    im.set_array(sample[0, i])
    return im,

# Create animation
ani = animation.FuncAnimation(fig, updatefig, interval=50, blit=True, frames=sample.shape[1])
HTML(ani.to_html5_video())


    

    
Out[6]:

In [ ]: