compute time frequency in label in source space with MNE

doc/source/getting_started.rst

@@ -52,26 +52,6 @@ If you get a new prompt with no error messages, you should be good to go.
 Learning Python
 ---------------
 
-If you are new to Python here are some online courses:
+If you are new to Python here is a very good place to get started:
 
-    * http://docs.python.org/
-
-    * http://learnpythonthehardway.org/book/
-
-    * http://code.google.com/intl/fr/edu/languages/google-python-class/
-
-    * http://homepage.mac.com/s_lott/books/python.html
-
-More specific tutorials on scientific computing with Python using Numpy and Scipy see:
-
-    * http://showmedo.com/videotutorials/numpy (video)
-
-    * http://www.scipy.org/NumPy_for_Matlab_Users
-
-    * http://mathesaurus.sourceforge.net/matlab-numpy.html
-
-Some discussions online are:
-
-    * http://stackoverflow.com/questions/111857/what-did-you-use-to-teach-yourself-python
-
-    * http://stackoverflow.com/questions/17988/how-to-learn-python
+    * http://scipy-lectures.github.com

examples/time_frequency/plot_source_label_time_frequency.py

@@ -0,0 +1,84 @@
+"""
+=========================================================
+Compute power and phase lock in label of the source space
+=========================================================
+
+Returns time-frequency maps of induced power and phase lock
+in the source space. The inverse method is linear based on dSPM inverse operator.
+
+"""
+# Authors: Alexandre Gramfort <gramfort@nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import numpy as np
+
+import mne
+from mne import fiff
+from mne.datasets import sample
+from mne.minimum_norm import read_inverse_operator, source_induced_power
+
+###############################################################################
+# Set parameters
+data_path = sample.data_path('..')
+raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
+label_name = 'Aud-lh'
+fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name
+
+tmin, tmax, event_id = -0.2, 0.5, 1
+
+# Setup for reading the raw data
+raw = fiff.Raw(raw_fname)
+events = mne.find_events(raw)
+inverse_operator = read_inverse_operator(fname_inv)
+
+include = []
+exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053']  # bads + 2 more
+
+# picks MEG gradiometers
+picks = fiff.pick_types(raw.info, meg=True, eeg=False, eog=True,
+                                stim=False, include=include, exclude=exclude)
+
+# Load condition 1
+event_id = 1
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
+                    preload=True)
+
+# Compute a source estimate per frequency band
+frequencies = np.arange(7, 30, 2)  # define frequencies of interest
+label = mne.read_label(fname_label)
+power, phase_lock = source_induced_power(epochs, inverse_operator, frequencies,
+                            label, baseline=(-0.1, 0), baseline_mode='percent',
+                            n_cycles=2, n_jobs=1)
+
+power = np.mean(power, axis=0)  # average over sources
+phase_lock = np.mean(phase_lock, axis=0)  # average over sources
+times = epochs.times
+
+###############################################################################
+# View time-frequency plots
+import pylab as pl
+pl.clf()
+pl.subplots_adjust(0.1, 0.08, 0.96, 0.94, 0.2, 0.43)
+pl.subplot(2, 1, 1)
+pl.imshow(20 * power, extent=[times[0], times[-1],
+                                      frequencies[0], frequencies[-1]],
+          aspect='auto', origin='lower')
+pl.xlabel('Time (s)')
+pl.ylabel('Frequency (Hz)')
+pl.title('Induced power in %s' % label_name)
+pl.colorbar()
+
+pl.subplot(2, 1, 2)
+pl.imshow(phase_lock, extent=[times[0], times[-1],
+                              frequencies[0], frequencies[-1]],
+          aspect='auto', origin='lower')
+pl.xlabel('Time (s)')
+pl.ylabel('Frequency (Hz)')
+pl.title('Phase-lock in %s' % label_name)
+pl.colorbar()
+pl.show()

examples/time_frequency/plot_source_space_time_frequency.py

@@ -17,7 +17,7 @@
 import mne
 from mne import fiff
 from mne.datasets import sample
-from mne.minimum_norm import read_inverse_operator, source_induced_power
+from mne.minimum_norm import read_inverse_operator, source_band_induced_power
 
 ###############################################################################
 # Set parameters
@@ -48,8 +48,8 @@
 # Compute a source estimate per frequency band
 bands = dict(alpha=[9, 11], beta=[18, 22])
 
-stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
-                            use_fft=False, n_jobs=-1)
+stcs = source_band_induced_power(epochs, inverse_operator, bands, n_cycles=2,
+                                 use_fft=False, n_jobs=1)
 
 for b, stc in stcs.iteritems():
     stc.save('induced_power_%s' % b)

mne/baseline.py

@@ -29,7 +29,7 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True):
         If baseline is equal ot (None, None) all the time
         interval is used.
 
-    mode: 'logratio' | 'ratio' | 'zscore' | 'mean'
+    mode: 'logratio' | 'ratio' | 'zscore' | 'mean' | 'percent'
         Do baseline correction with ratio (power is divided by mean
         power during baseline) or zscore (power is divided by standard
         deviatio of power during baseline after substracting the mean,
@@ -44,7 +44,7 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True):
     if copy:
         data = data.copy()
 
-    valid_modes = ['logratio', 'ratio', 'zscore', 'mean']
+    valid_modes = ['logratio', 'ratio', 'zscore', 'mean', 'percent']
     if mode not in valid_modes:
         raise Exception('mode should be any of : %s' % valid_modes)
 
@@ -73,6 +73,9 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True):
             std = np.std(data[..., imin:imax], axis=-1)[..., None]
             data -= mean
             data /= std
+        elif mode == 'percent':
+            data -= mean
+            data /= mean
 
     elif verbose:
         print "No baseline correction applied..."

mne/minimum_norm/__init__.py

@@ -1,4 +1,4 @@
 from .inverse import read_inverse_operator, apply_inverse, \
                      apply_inverse_raw, make_inverse_operator, \
                      apply_inverse_epochs
-from .time_frequency import source_induced_power
+from .time_frequency import source_band_induced_power, source_induced_power

mne/minimum_norm/tests/test_time_frequency.py

@@ -1,13 +1,14 @@
 import os.path as op
 
 import numpy as np
-from numpy.testing import assert_array_almost_equal, assert_equal
+from numpy.testing import assert_array_almost_equal
 from nose.tools import assert_true
 
 from ...datasets import sample
 from ... import fiff, find_events, Epochs
+from ...label import read_label
 from ..inverse import read_inverse_operator
-from ..time_frequency import source_induced_power
+from ..time_frequency import source_band_induced_power, source_induced_power
 
 
 examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples')
@@ -16,10 +17,11 @@
                                         'sample_audvis-meg-oct-6-meg-inv.fif')
 fname_data = op.join(data_path, 'MEG', 'sample',
                                         'sample_audvis_raw.fif')
+fname_label = op.join(data_path, 'MEG', 'sample', 'labels', 'Aud-lh.label')
 
 
 def test_tfr_with_inverse_operator():
-    """Test MNE inverse computation"""
+    """Test time freq with MNE inverse computation"""
 
     tmin, tmax, event_id = -0.2, 0.5, 1
 
@@ -33,27 +35,41 @@ def test_tfr_with_inverse_operator():
 
     # picks MEG gradiometers
     picks = fiff.pick_types(raw.info, meg=True, eeg=False, eog=True,
-                                    stim=False, include=include, exclude=exclude)
+                                stim=False, include=include, exclude=exclude)
 
     # Load condition 1
     event_id = 1
-    events = events[:3]  # take 3 events to keep the computation time low
-    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
-                        baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
-                        preload=True)
+    events3 = events[:3]  # take 3 events to keep the computation time low
+    epochs = Epochs(raw, events3, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
+                    preload=True)
 
     # Compute a source estimate per frequency band
     bands = dict(alpha=[10, 10])
+    label = read_label(fname_label)
 
-    stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
-                                use_fft=False, pca=True)
+    stcs = source_band_induced_power(epochs, inverse_operator, bands,
+                            n_cycles=2, use_fft=False, pca=True, label=label)
 
     stc = stcs['alpha']
     assert_true(len(stcs) == len(bands.keys()))
     assert_true(np.all(stc.data > 0))
     assert_array_almost_equal(stc.times, epochs.times)
 
-    stcs_no_pca = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
-                                use_fft=False, pca=False)
+    stcs_no_pca = source_band_induced_power(epochs, inverse_operator, bands,
+                            n_cycles=2, use_fft=False, pca=False, label=label)
 
     assert_array_almost_equal(stcs['alpha'].data, stcs_no_pca['alpha'].data)
+
+    # Compute a source estimate per frequency band
+    epochs = Epochs(raw, events[:10], event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
+                    preload=True)
+
+    frequencies = np.arange(7, 30, 2)  # define frequencies of interest
+    power, phase_lock = source_induced_power(epochs, inverse_operator,
+                            frequencies, label, baseline=(-0.1, 0),
+                            baseline_mode='percent', n_cycles=2, n_jobs=1)
+    assert_true(np.all(phase_lock > 0))
+    assert_true(np.all(phase_lock < 1))
+    assert_true(np.max(power) > 10)