MNE : Demo of time-frequency analysis

Author : Alexandre Gramfort

date : Jan. 2014

DISCLAIMER : This tutorial is for illustration purposes. The function to perform the TFR analysis is likely to change over the next few months. Objective is to simplify the code required.

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 0x110464d40>
    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]: 
import numpy as np
times = np.linspace(0.07, 0.25, 5)
fig = evoked_left.plot_topomap(times=times, ch_type='mag')
fig = evoked_right.plot_topomap(times=times, ch_type='mag')
In [24]: 
fig = evoked_left.plot_topomap(times=times, ch_type='grad')
fig = evoked_right.plot_topomap(times=times, ch_type='grad')

Now look at the stimulus period

In [25]: 
event_id = dict(stim_left=12, stim_right=22)  # event trigger and conditions
tmin = -0.75  # start of each epoch
tmax = 0.8  # end of each epoch
baseline = (-0.5, 0)
In [37]: 
# 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=True, eeg=False, eog=True,
                            stim=False, exclude='bads', selection=None)
In [38]: 
epochs_stim = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                         picks=picks, baseline=baseline, reject=reject,
                         preload=True)  # with preload
print epochs_stim

<Epochs  |  n_events : 228 (all good), tmin : -0.75 (s), tmax : 0.8 (s), baseline : None,
 'stim_left': 120, 'stim_right': 108>

Let's look at the ERF on gradiometers only

In [39]: 
fig = epochs_stim.average().plot()

Let's try to find some visual gamma... time-frequency analysis

Let's not run TF on EOG channels:

In [40]: 
print epochs_stim.ch_names[:5]

['BIO001', 'BIO002', 'MEG0111', 'MEG0112', 'MEG0113']
In [41]: 
epochs_stim.drop_channels(['BIO001', 'BIO002'])
In [42]: 
print epochs_stim.ch_names[:5]  # EOGs are gone...

['MEG0111', 'MEG0112', 'MEG0113', 'MEG0121', 'MEG0122']

Define paramters for time-frequency analysis

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

(228, 306, 1551)

Let's compute

In [44]: 
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 [45]: 
power.shape, data.shape, frequencies.shape
Out[45]:
((306, 46, 311), (228, 306, 1551), (46,))

Let's visualize

In [46]: 
from mne.viz import plot_topo_power

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

power_viz = power.copy()
baseline = None
mode = None
mask = (epochs_stim.times[::decim] < -0.1) & (epochs_stim.times[::decim] > -0.5)
power_viz = power_viz / np.mean(power_viz[:,:,mask], axis=-1)[:, :, None] - 1.

title = 'Induced power - Stim ON'
plot_topo_power(epochs_stim, power_viz, frequencies, layout, baseline=baseline,
                mode=mode, decim=decim, title=title,
                dB=False, vmin=-.3, vmax=0.3)
plt.show()
In [49]: 
ch_idx = epochs_stim.ch_names.index('MEG1922')

times = 1e3 * epochs_stim.times[::decim]
extent = (times[0], times[-1], frequencies[0], frequencies[-1])
plt.imshow(power_viz[ch_idx], extent=extent, aspect="auto", origin="lower",
          vmin=-0.3, vmax=0.3, picker=True)

plt.xlim([-300., 700.])
plt.xlabel('Time (ms)')
plt.ylabel('Freq. (Hz)')
plt.title('Power ' + epochs.ch_names[ch_idx])
plt.colorbar()
Out[49]:
<matplotlib.colorbar.Colorbar instance at 0x1547c5878>
In [51]: 
%matplotlib inline

plt.imshow(phase_lock[ch_idx], extent=extent, aspect="auto", origin="lower",
          vmin=0, vmax=0.8, picker=True)

plt.xlim([-300., 700.])
plt.xlabel('Time (ms)')
plt.ylabel('Freq. (Hz)')
plt.title('Phase Lock ' + epochs.ch_names[ch_idx])
plt.colorbar()
Out[51]:
<matplotlib.colorbar.Colorbar instance at 0x11411b3f8>
In []: