[med-svn] [python-mne] 52/376: adding way to estimate noise cov matrix in python

[Yaroslav Halchenko]

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

yoh pushed a commit to annotated tag v0.1
in repository python-mne.

commit f35a9ba95855c7d508d9fd66ed309e69b3e5153e
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date:   Tue Jan 25 15:32:09 2011 -0500

    adding way to estimate noise cov matrix in python
---
 examples/plot_estimate_covariance_matrix.py | 52 +++++++++++++++++++++++++++++
 examples/read_write_raw.py                  | 16 ++++-----
 mne/cov.py                                  | 33 +++++++++++++++++-
 3 files changed, 91 insertions(+), 10 deletions(-)

diff --git a/examples/plot_estimate_covariance_matrix.py b/examples/plot_estimate_covariance_matrix.py
new file mode 100644
index 0000000..63fe3e3
--- /dev/null
+++ b/examples/plot_estimate_covariance_matrix.py
@@ -0,0 +1,52 @@
+"""
+==============================================
+Estimate covariance matrix from a raw FIF file
+==============================================
+
+"""
+# Author: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import os
+import mne
+from mne import fiff
+
+fname = os.environ['MNE_SAMPLE_DATASET_PATH']
+fname += '/MEG/sample/sample_audvis_raw.fif'
+
+raw = fiff.setup_read_raw(fname)
+
+# Set up pick list: MEG + STI 014 - bad channels
+want_meg = True
+want_eeg = False
+want_stim = False
+
+picks = fiff.pick_types(raw['info'], meg=want_meg, eeg=want_eeg,
+                        stim=want_stim, exclude=raw['info']['bads'])
+
+print "Number of picked channels : %d" % len(picks)
+
+full_cov = mne.Covariance(kind='full')
+full_cov.estimate_from_raw(raw, picks=picks)
+print full_cov
+
+diagonal_cov = mne.Covariance(kind='diagonal')
+diagonal_cov.estimate_from_raw(raw, picks=picks)
+print diagonal_cov
+
+###############################################################################
+# Show covariance
+import pylab as pl
+pl.figure(figsize=(10, 5))
+pl.subplot(1, 2, 1)
+pl.imshow(full_cov.data, interpolation="nearest")
+pl.title('Full covariance matrix (%d channels)' % full_cov.data.shape[0])
+pl.subplot(1, 2, 2)
+pl.imshow(diagonal_cov.data, interpolation="nearest")
+pl.title('Diagonal covariance matrix (%d channels)' % \
+                                              diagonal_cov.data.shape[0])
+pl.subplots_adjust(0.06, 0.02, 0.98, 0.94, 0.16, 0.26)
+pl.show()
diff --git a/examples/read_write_raw.py b/examples/read_write_raw.py
index 59456cb..09293ee 100644
--- a/examples/read_write_raw.py
+++ b/examples/read_write_raw.py
@@ -11,10 +11,12 @@ Read and write raw data in 60-sec blocks
 
 print __doc__
 
+import os
 from math import ceil
 from mne import fiff
 
-infile = 'MNE-sample-data/MEG/sample/sample_audvis_raw.fif'
+infile = os.environ['MNE_SAMPLE_DATASET_PATH']
+infile += '/MEG/sample/sample_audvis_raw.fif'
 outfile = 'sample_audvis_small_raw.fif'
 
 raw = fiff.setup_read_raw(infile)
@@ -51,15 +53,11 @@ quantum = int(ceil(quantum_sec * raw['info']['sfreq']))
 #
 first_buffer = True
 for first in range(start, stop, quantum):
-    last = start + quantum
+    last = first + quantum
     if last >= stop:
-        last = stop
-    try:
-        data, times = raw[picks, first:last]
-    except Exception as inst:
-        raw['fid'].close()
-        outfid.close()
-        print inst
+        last = stop+1
+
+    data, times = raw[picks, first:last]
 
     print 'Writing ... ',
     fiff.write_raw_buffer(outfid, data, cals)
diff --git a/mne/cov.py b/mne/cov.py
index f08c808..4439aae 100644
--- a/mne/cov.py
+++ b/mne/cov.py
@@ -24,7 +24,7 @@ class Covariance(object):
     _kind_to_id = dict(full=1, sparse=2, diagonal=3) # XXX : check
     _id_to_kind = {1: 'full', 2: 'sparse', 3: 'diagonal'} # XXX : check
 
-    def __init__(self, kind):
+    def __init__(self, kind='full'):
         self.kind = kind
 
     def load(self, fname):
@@ -48,6 +48,37 @@ class Covariance(object):
         """save covariance matrix in a FIF file"""
         write_cov_file(fname, self._cov)
 
+    def estimate_from_raw(self, raw, picks=None, quantum_sec=10):
+        """Estimate noise covariance matrix from a raw FIF file
+        """
+        #   Set up the reading parameters
+        start = raw['first_samp']
+        stop = raw['last_samp'] + 1
+        quantum = int(quantum_sec * raw['info']['sfreq'])
+
+        cov = 0
+        n_samples = 0
+
+        # Read data
+        for first in range(start, stop, quantum):
+            last = first + quantum
+            if last >= stop:
+                last = stop
+
+            data, times = raw[picks, first:last]
+
+            if self.kind is 'full':
+                cov += np.dot(data, data.T)
+            elif self.kind is 'diagonal':
+                cov += np.diag(np.sum(data ** 2, axis=1))
+            else:
+                raise ValueError, "Unsupported covariance kind"
+
+            n_samples += data.shape[1]
+
+        self.data = cov / n_samples # XXX : check
+        print '[done]'
+
     def __repr__(self):
         s = "kind : %s" % self.kind
         s += ", size : %s x %s" % self.data.shape

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