ENH: multi-taper PSD estimation

doc/source/whats_new.rst

@@ -9,6 +9,8 @@ Current
 Changelog
 ~~~~~~~~~
 
+   - Multi-taper PSD estimation for single epochs in source space using minimum norm by `Martin Luessi`_
+
    - Read and visualize .dip files obtained with xfit or mne_dipole_fit by `Alex Gramfort`_
 
 .. _changes_0_4:

examples/time_frequency/plot_compute_source_psd_epochs.py

@@ -0,0 +1,88 @@
+"""
+=====================================================================
+Compute Power Spectral Density of inverse solution from single epochs
+=====================================================================
+
+Compute PSD of dSPM inverse solution on single trial epochs restricted
+to a brain label. The PSD is computed using a multi-taper method with
+Discrete Prolate Spheroidal Sequence (DPSS) windows.
+
+"""
+
+# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import numpy as np
+import pylab as pl
+import mne
+from mne.datasets import sample
+from mne.fiff import Raw, pick_types
+from mne.minimum_norm import read_inverse_operator, compute_source_psd_epochs
+
+
+data_path = sample.data_path('..')
+fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
+fname_raw = data_path + '/MEG/sample/sample_audvis_raw.fif'
+fname_event = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'
+label_name = 'Aud-lh'
+fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name
+
+event_id, tmin, tmax = 1, -0.2, 0.5
+snr = 1.0  # use smaller SNR for raw data
+lambda2 = 1.0 / snr ** 2
+method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)
+
+# Load data
+inverse_operator = read_inverse_operator(fname_inv)
+label = mne.read_label(fname_label)
+raw = Raw(fname_raw)
+events = mne.read_events(fname_event)
+
+# Set up pick list
+include = []
+exclude = raw.info['bads'] + ['EEG 053']  # bads + 1 more
+
+# pick MEG channels
+picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,
+                                            include=include, exclude=exclude)
+# Read epochs
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,
+                                                    eog=150e-6))
+
+# define frequencies of interest
+fmin, fmax = 0., 70.
+bandwidth = 4.  # bandwidth of the windows in Hz
+
+# compute source space psd in label
+
+# Note: By using "return_generator=True" stcs will be a generator object
+# instead of a list. This allows us so to iterate without having to
+# keep everything in memory.
+
+stcs = compute_source_psd_epochs(epochs, inverse_operator, lambda2=lambda2,
+                                 method=method, fmin=fmin, fmax=fmax,
+                                 bandwidth=bandwidth, return_generator=True)
+
+# compute average PSD over the first 10 epochs
+n_epochs = 10
+for i, stc in enumerate(stcs):
+    if i >= n_epochs:
+        break
+
+    if i == 0:
+        psd_avg = np.mean(stc.data, axis=0)
+    else:
+        psd_avg += np.mean(stc.data, axis=0)
+
+psd_avg /= n_epochs
+freqs = stc.times  # the frequencies are stored here
+
+pl.figure()
+pl.plot(freqs, psd_avg)
+pl.xlabel('Freq (Hz)')
+pl.ylabel('Power Spectral Density')
+pl.show()

mne/minimum_norm/__init__.py

@@ -4,4 +4,4 @@
                      apply_inverse_raw, make_inverse_operator, \
                      apply_inverse_epochs, write_inverse_operator
 from .time_frequency import source_band_induced_power, source_induced_power, \
-                            compute_source_psd
+                            compute_source_psd, compute_source_psd_epochs

mne/minimum_norm/tests/test_time_frequency.py

@@ -7,11 +7,15 @@
 from mne.datasets import sample
 from mne import fiff, find_events, Epochs
 from mne.label import read_label
-from mne.minimum_norm.inverse import read_inverse_operator
+from mne.minimum_norm.inverse import read_inverse_operator, \
+                                     apply_inverse_epochs
 from mne.minimum_norm.time_frequency import source_band_induced_power, \
-                            source_induced_power, compute_source_psd
+                            source_induced_power, compute_source_psd, \
+                            compute_source_psd_epochs
 
 
+from mne.time_frequency import multitaper_psd
+
 examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples')
 data_path = sample.data_path(examples_folder)
 fname_inv = op.join(data_path, 'MEG', 'sample',
@@ -93,3 +97,61 @@ def test_source_psd():
     # Time max at line frequency (60 Hz in US)
     assert_true(59e-3 <= stc.times[np.argmax(np.sum(stc.data, axis=0))]
                       <= 61e-3)
+
+
+def test_source_psd_epochs():
+    """Test multi-taper source PSD computation in label from epochs"""
+
+    raw = fiff.Raw(fname_data)
+    inverse_operator = read_inverse_operator(fname_inv)
+    label = read_label(fname_label)
+
+    event_id, tmin, tmax = 1, -0.2, 0.5
+    lambda2, method = 1. / 9., 'dSPM'
+    bandwidth = 8.
+    fmin, fmax = 0, 100
+
+    picks = fiff.pick_types(raw.info, meg=True, eeg=False, stim=True,
+                            ecg=True, eog=True, include=['STI 014'],
+                            exclude=raw.info['bads'])
+    reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)
+
+    events = find_events(raw)
+    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=reject)
+
+    # only look at one epoch
+    epochs.drop_bad_epochs()
+    one_epochs = epochs[:1]
+
+    # return list
+    stc_psd = compute_source_psd_epochs(one_epochs, inverse_operator,
+                                        lambda2=lambda2, method=method,
+                                        pick_normal=True, label=label,
+                                        bandwidth=bandwidth,
+                                        fmin=fmin, fmax=fmax)[0]
+
+    # return generator
+    stcs = compute_source_psd_epochs(one_epochs, inverse_operator,
+                                     lambda2=lambda2, method=method,
+                                     pick_normal=True, label=label,
+                                     bandwidth=bandwidth,
+                                     fmin=fmin, fmax=fmax,
+                                     return_generator=True)
+
+    for stc in stcs:
+        stc_psd_gen = stc
+
+    assert_array_almost_equal(stc_psd.data, stc_psd_gen.data)
+
+    # compare with direct computation
+    stc = apply_inverse_epochs(one_epochs, inverse_operator,
+                               lambda2=lambda2, method=method,
+                               pick_normal=True, label=label)[0]
+
+    sfreq = epochs.info['sfreq']
+    psd, freqs = multitaper_psd(stc.data, sfreq=sfreq, bandwidth=bandwidth,
+                                fmin=fmin, fmax=fmax)
+
+    assert_array_almost_equal(psd, stc_psd.data)
+    assert_array_almost_equal(freqs, stc_psd.times)