[med-svn] [python-mne] 234/353: ENH: raw.add_proj, updated example

[Yaroslav Halchenko]

This is an automated email from the git hooks/post-receive script.

yoh pushed a commit to tag 0.4
in repository python-mne.

commit fd0d62ba1ad13e40bd41f8456019152a48ea3cf3
Author: mluessi at nmr.mgh.harvard.edu <mluessi at nmr.mgh.harvard.edu>
Date:   Wed Jun 27 03:00:56 2012 -0400

    ENH: raw.add_proj, updated example
---
 .../plot_compute_raw_data_spectrum.py              | 49 +++++++++++++++-------
 mne/fiff/raw.py                                    | 25 +++++++++++
 2 files changed, 58 insertions(+), 16 deletions(-)

diff --git a/examples/time_frequency/plot_compute_raw_data_spectrum.py b/examples/time_frequency/plot_compute_raw_data_spectrum.py
index 20c6bd6..cfbcc2a 100644
--- a/examples/time_frequency/plot_compute_raw_data_spectrum.py
+++ b/examples/time_frequency/plot_compute_raw_data_spectrum.py
@@ -1,21 +1,21 @@
 """
 ==================================================
-Compute the spectrum of raw data using multi-taper
+Compute the power spectral density of raw data
 ==================================================
 
 This script shows how to compute the power spectral density (PSD)
-of measurements on a raw dataset.
-
+of measurements on a raw dataset. It also show the effect of applying SSP
+to the data to reduce ECG and EOG artifacts.
 """
 # Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
-#
+#          Martin Luessi <mluessi at nmr.mgh.harvard.edu>
 # License: BSD (3-clause)
 
 print __doc__
 
 import numpy as np
 
-from mne import fiff
+from mne import fiff, read_proj, read_selection
 from mne.time_frequency import compute_raw_psd
 from mne.datasets import sample
 
@@ -23,37 +23,54 @@ from mne.datasets import sample
 # Set parameters
 data_path = sample.data_path('..')
 raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+proj_fname = data_path + '/MEG/sample/sample_audvis_eog_proj.fif'
 
 # Setup for reading the raw data
 raw = fiff.Raw(raw_fname)
 exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053']  # bads + 2 more
 
-# picks MEG gradiometers
-picks = fiff.pick_types(raw.info, meg='grad', eeg=False, eog=False,
-                        stim=False, exclude=exclude)
+# Add SSP projection vectors to reduce EOG and ECG artifacts
+projs = read_proj(proj_fname)
+raw.add_proj(projs, activate=True, remove_existing=True)
+
+# Pick MEG magnetometers in the Left-temporal region
+selection = read_selection('Left-temporal')
+picks = fiff.pick_types(raw.info, meg='mag', eeg=False, eog=False,
+                        stim=False, exclude=exclude, selection=selection)
 
 tmin, tmax = 0, 60  # use the first 60s of data
 fmin, fmax = 2, 300  # look at frequencies between 5 and 70Hz
-NFFT = 2048 # the FFT size (NFFT). Ideally a power of 2
+NFFT = 2048  # the FFT size (NFFT). Ideally a power of 2
 psds, freqs = compute_raw_psd(raw, tmin=tmin, tmax=tmax, picks=picks,
                               fmin=fmin, fmax=fmax, NFFT=NFFT, n_jobs=1,
-                              plot=False)
+                              plot=False, proj=False)
+
+# And now do the same with SSP applied
+psds_ssp, freqs = compute_raw_psd(raw, tmin=tmin, tmax=tmax, picks=picks,
+                                  fmin=fmin, fmax=fmax, NFFT=NFFT, n_jobs=1,
+                                  plot=False, proj=True)
 
 # Convert PSDs to dB
 psds = 10 * np.log10(psds)
+psds_ssp = 10 * np.log10(psds_ssp)
 
 ###############################################################################
 # Compute mean and standard deviation accross channels and then plot
-psd_mean = np.mean(psds, axis=0)
-psd_std = np.std(psds, axis=0)
+def plot_psds(freqs, psds, fill_color):
+    psd_mean = np.mean(psds, axis=0)
+    psd_std = np.std(psds, axis=0)
+    hyp_limits = (psd_mean - psd_std, psd_mean + psd_std)
 
-hyp_limits = (psd_mean - psd_std, psd_mean + psd_std)
+    pl.plot(freqs, psd_mean)
+    pl.fill_between(freqs, hyp_limits[0], y2=hyp_limits[1], color=fill_color,
+                    alpha=0.5)
 
 import pylab as pl
 pl.figure()
-pl.plot(freqs, psd_mean)
-pl.fill_between(freqs, hyp_limits[0], y2=hyp_limits[1], color=(1, 0, 0, .3),
-                alpha=0.5)
+plot_psds(freqs, psds, (0, 0, 1, .3))
+plot_psds(freqs, psds_ssp, (0, 1, 0, .3))
 pl.xlabel('Freq (Hz)')
 pl.ylabel('Power Spectral Density (dB/Hz)')
+pl.legend(['Without SSP', 'With SSP'])
 pl.show()
+
diff --git a/mne/fiff/raw.py b/mne/fiff/raw.py
index 5c5ee3e..7be18ca 100644
--- a/mne/fiff/raw.py
+++ b/mne/fiff/raw.py
@@ -17,6 +17,7 @@ from .meas_info import read_meas_info, write_meas_info
 from .tree import dir_tree_find
 from .tag import read_tag
 from .pick import pick_types
+from .proj import activate_proj
 
 from ..filter import low_pass_filter, high_pass_filter, band_pass_filter
 from ..parallel import parallel_func
@@ -535,6 +536,30 @@ class Raw(object):
         self.filter(None, freq, picks, n_jobs=n_jobs, verbose=verbose,
                     filter_length=filter_length)
 
+    def add_proj(self, projs, activate=True, remove_existing=False):
+        """Add SSP projection vectors
+
+        Parameters
+        ----------
+        projs : list
+            List with projection vectors
+
+        activate : bool
+            Activate the projection vectors being added
+
+        remove_existing : bool
+            Remove the projection vectors currently in the file
+        """
+        if activate:
+            projs = activate_proj(projs, copy=True)
+        else:
+            projs = copy.deepcopy(projs)
+
+        if remove_existing:
+            self.info['projs'] = projs
+        else:
+            self.info['projs'].extend(projs)
+
     def save(self, fname, picks=None, tmin=0, tmax=None, buffer_size_sec=10,
              drop_small_buffer=False):
         """Save raw data to file

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git