MNE : Demo of time-frequency analysis

In [1]: 
# add plot inline in the page
%matplotlib inline

First, load the mne package:

In [2]: 
import mne

We set the log-level to 'WARNING' so the output is less verbose

In [3]: 
mne.set_log_level('WARNING')

Access raw data

Now we import the sample dataset. If you don't already have it, it will be downloaded automatically (but be patient approx. 2GB)

In [4]: 
data_path = '/Users/alex/Sync/karolinska_teaching/'
raw_fname = data_path + '/meg/workshop_visual_sss.fif'

print raw_fname

/Users/alex/Sync/karolinska_teaching//meg/workshop_visual_sss.fif

Read data from file:

In [5]: 
raw = mne.fiff.Raw(raw_fname)
print raw

<Raw  |  n_channels x n_times : 322 x 988000>

!!! we need to fix now tiny acquisition problems

In [6]: 
def fix_info(raw):
    raw.info['chs'][raw.ch_names.index('BIO001')]['kind'] = mne.fiff.constants.FIFF.FIFFV_EOG_CH
    raw.info['chs'][raw.ch_names.index('BIO002')]['kind'] = mne.fiff.constants.FIFF.FIFFV_EOG_CH
    raw.info['chs'][raw.ch_names.index('BIO003')]['kind'] = mne.fiff.constants.FIFF.FIFFV_ECG_CH

fix_info(raw)
In [7]: 
print raw.info

<Info | 24 non-empty fields
    acq_pars : unicode | 25251 items
    bads : list | 0 items
    buffer_size_sec : numpy.float64 | 1.0
    ch_names : list | BIO001, BIO002, BIO003, IASX+, IASX-, IASY+, IASY-, ...
    chs : list | 322 items
    comps : list | 0 items
    description : unicode | 17 items
    dev_head_t : dict | 3 items
    dig : list | 207 items
    experimenter : unicode | 15 items
    file_id : dict | 4 items
    filename : str | /Users/ale.../workshop_visual_sss.fif
    filenames : list | 1 items
    highpass : float | 0.10000000149
    lowpass : float | 330.0
    meas_date : numpy.ndarray | 2013-12-16 09:58:50
    meas_id : dict | 4 items
    nchan : int | 322
    orig_blocks : list | 5 items
    orig_fid_str : instance | <StringIO.StringIO instance at 0x1105c1878>
    proj_id : numpy.ndarray | 1 items
    proj_name : unicode | 8 items
    projs : list | 0 items
    sfreq : float | 1000.0
    acq_stim : NoneType
    ctf_head_t : NoneType
    dev_ctf_t : NoneType
>
In [8]: 
print raw.info['bads']

[]

Define and read epochs

First extract events:

In [9]: 
events = mne.find_events(raw, stim_channel='STI101', verbose=True)

Reading 0 ... 987999  =      0.000 ...   987.999 secs...
[done]
1473 events found
Events id: [10 11 12 13 14 16 17 20 21 22 23 24 26 27]
In [10]: 
print events[:5]  # events is a 2d array

[[172624      0     10]
 [173142      0     11]
 [173892      0     12]
 [175201      0     13]
 [175460      0     14]]

Look at the design in a graphical way:

In [11]: 
plt.plot((events[:, 0] - raw.first_samp) / raw.info['sfreq'], events[:, 2], '.');
plt.xlabel('time (s)')
plt.axis([0.1168, 41.4003, 9.7295, 17.3887])
Out[11]:
[0.1168, 41.4003, 9.7295, 17.3887]

From raw to epochs

Define epochs parameters:

In [12]: 
event_id = dict(cue_left=10, cue_right=20)  # event trigger and conditions
tmin = -0.3  # start of each epoch
tmax = 1.6  # end of each epoch
In [13]: 
event_id
Out[13]:
{'cue_left': 10, 'cue_right': 20}

Pick the good channels:

In [14]: 
picks = mne.fiff.pick_types(raw.info, meg=True, eeg=False, eog=True,
                            stim=False, exclude='bads')

Define the baseline period:

In [15]: 
baseline = (None, 0)  # means from the first instant to t = 0

Define peak-to-peak rejection parameters for gradiometers, magnetometers and EOG:

In [16]: 
reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)

Read epochs:

In [17]: 
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=picks, baseline=baseline, reject=reject,
                    preload=True)  # with preload
print epochs

<Epochs  |  n_events : 126 (all good), tmin : -0.3 (s), tmax : 1.6 (s), baseline : (None, 0),
 'cue_left': 68, 'cue_right': 58>
In [18]: 
epochs.plot_drop_log()
Out[18]:
91.446028513238289

Create trialinfo array like Fieldtrip

In [19]: 
def make_trialinfo(epochs, events):
    trialinfo = []
    n_epochs = len(epochs.selection)
    sfreq = epochs.info['sfreq']
    for ii in epochs.selection:
        this_trialinfo = []
        for jj in range(ii + 1, ii + 5):
            if jj >= len(events):
                continue
            dt = (events[jj, 0] - events[ii, 0]) / sfreq
            if dt < epochs.tmax:
                this_trialinfo.append(dt)
        trialinfo.append(this_trialinfo)
    return np.array(trialinfo)

trialinfo = make_trialinfo(epochs, events)

Plot ERF on gradiometers and average latencies

In [20]: 
fig = mne.fiff.pick_types_evoked(epochs.average(), meg='grad').plot(show=False)
ax = plt.gca()
ax.axvline(0., color='b')
ax.axvline(1000. * np.mean(trialinfo[:, 0]), color='g')
ax.axvline(1000. * np.mean(trialinfo[:, 1]), color='r')
plt.show()
In [21]: 
evoked_left = epochs['cue_left'].average()
print evoked_left

evoked_right = epochs['cue_right'].average()
print evoked_right

<Evoked  |  comment : 'cue_left', time : [-0.300000, 1.600000], n_epochs : 68, n_channels x n_times : 306 x 1901>
<Evoked  |  comment : 'cue_right', time : [-0.300000, 1.600000], n_epochs : 58, n_channels x n_times : 306 x 1901>

Plot ERFs for both and then each conditions:

In [22]: 
fig = (evoked_right + evoked_left).plot()  # left and right
In [23]: 
fig = evoked_right.plot()
In [24]: 
fig = evoked_left.plot()
In [25]: 
contrast = evoked_left - evoked_right
fig = contrast.plot()
In [26]: 
import numpy as np
fig = evoked_left.plot_topomap(times=np.linspace(0.07, 0.25, 5), ch_type='mag')
fig = evoked_right.plot_topomap(times=np.linspace(0.07, 0.25, 5), ch_type='mag')
In [27]: 
fig = evoked_left.plot_topomap(times=np.linspace(0.07, 0.25, 5), ch_type='grad')
fig = evoked_right.plot_topomap(times=np.linspace(0.07, 0.25, 5), ch_type='grad')
In [30]: 
from mne.viz import plot_topo
layout = mne.find_layout(evoked_left.info)
plot_topo([evoked_left, evoked_right], layout=layout, color=['yellow', 'red']);

Now look at the response

In [31]: 
event_id = dict(speedup_left=13, speedup_right=23)  # event trigger and conditions
tmin = -0.1  # start of each epoch
tmax = 0.3  # end of each epoch
baseline = (None, 0)
In [32]: 
epochs_speedup = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                            picks=picks, baseline=baseline, reject=reject,
                            preload=True)  # with preload
print epochs_speedup

<Epochs  |  n_events : 280 (all good), tmin : -0.1 (s), tmax : 0.3 (s), baseline : (None, 0),
 'speedup_left': 149, 'speedup_right': 131>
In [33]: 
trialinfo = make_trialinfo(epochs_speedup, events)
In [34]: 
plt.plot(np.sort(1000. * trialinfo[:, 0]))
plt.xlim(0, 140.)
Out[34]:
(0, 140.0)

Let's look at the ERF on gradiometers only

In [102]: 
fig = mne.fiff.pick_types_evoked(epochs_speedup.average(), meg='grad').plot(show=False)
ax = plt.gca()
ax.axvline(0., color='b')
ax.axvline(1000. * np.mean(trialinfo[:, 0]), color='g')
plt.show()

Let's epoch on the response

In [103]: 
event_id = dict(resp_left=14, resp_right=24)  # event trigger and conditions
tmin = -0.75  # start of each epoch
tmax = 1.5  # end of each epoch
baseline = (-0.3, -0.15)

if 'mag' in reject:
    del reject['mag']  # get rid of rejection value for magnetometers

# Restrict the analysis to occipital sensors for speed
selection = mne.read_selection('occipital')
selection = [s.replace(' ', '') for s in selection]  # remove empty spaces

picks = mne.fiff.pick_types(raw.info, meg='grad', eeg=False, eog=True,
                            stim=False, exclude='bads',
                            selection=selection)

epochs_resp = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                         picks=picks, baseline=baseline, reject=reject,
                         preload=True)  # with preload
print epochs_resp

<Epochs  |  n_events : 60 (all good), tmin : -0.75 (s), tmax : 1.5 (s), baseline : (-0.3, -0.15),
 'resp_left': 33, 'resp_right': 27>

Let's average and look at the ERF

In [104]: 
evoked_resp = epochs_resp.average()
In [105]: 
fig = evoked_resp.plot()

Run time-frequency analysis on response

Subtract evoked response for frequency analysis

In [87]: 
epochs_resp = epochs_resp.subtract_evoked()

Let's not run TF on EOG channels:

In [106]: 
print epochs_resp.ch_names[:5]

['BIO001', 'BIO002', 'MEG1642', 'MEG1643', 'MEG1712']
In [107]: 
epochs_resp.drop_picks([0, 1])
In [108]: 
print epochs_resp.ch_names[:5]

['MEG1642', 'MEG1643', 'MEG1712', 'MEG1713', 'MEG1722']

Define paramters for time-frequency analysis

In [109]: 
frequencies = np.arange(5., 90., 2.5)  # define frequencies of interest
n_cycles = frequencies / 7.  # different number of cycle per frequency
Fs = raw.info['sfreq']  # sampling in Hz
decim = 5
data = epochs_resp.get_data()
print data.shape

(60, 48, 2251)

Let's compute

In [110]: 
from mne.time_frequency import induced_power
power, phase_lock = induced_power(data, Fs=Fs, frequencies=frequencies,
                                  n_cycles=n_cycles, n_jobs=4, use_fft=False,
                                  decim=decim, zero_mean=True)

Let's look at some shapes to understand what we manipulate

In [111]: 
power.shape, data.shape, frequencies.shape
Out[111]:
((48, 34, 451), (60, 48, 2251), (34,))

Let's plot the power

In [112]: 
# %matplotlib qt4
%matplotlib inline

from mne.viz import plot_topo_power

baseline = (-0.2, 0)  # set the baseline for induced power
mode = 'percent'  # set mode for baseline rescaling
layout = mne.find_layout(epochs_resp.info)

title = 'Induced power - Response'
plot_topo_power(epochs_resp, power, frequencies, layout, baseline=baseline,
                mode=mode, decim=decim, vmin=-3., vmax=3, title=title,
                dB=False)
plt.show()
In []: