[med-svn] [python-mne] 182/376: ENH : finish artefact rejection in cov estimation
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:32 UTC 2015
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 ccd178e2cc4c134362ab838696162b2071202601
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Mon Apr 4 11:13:26 2011 -0400
ENH : finish artefact rejection in cov estimation
---
examples/plot_estimate_covariance_matrix.py | 2 +-
mne/cov.py | 123 ++++++++++++++++++++++++----
mne/epochs.py | 89 ++++++++++----------
mne/fiff/pick.py | 11 +++
mne/fiff/tests/data/test-cov.fif | Bin 541025 -> 541025 bytes
mne/tests/test_cov.py | 13 +--
mne/time_frequency/tests/test_tfr.py | 4 +-
7 files changed, 170 insertions(+), 72 deletions(-)
diff --git a/examples/plot_estimate_covariance_matrix.py b/examples/plot_estimate_covariance_matrix.py
index 5b41bb5..7f595db 100644
--- a/examples/plot_estimate_covariance_matrix.py
+++ b/examples/plot_estimate_covariance_matrix.py
@@ -20,7 +20,7 @@ fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
raw = fiff.Raw(fname)
# Set up pick list: MEG + STI 014 - bad channels
-cov = mne.noise_covariance_segment(raw)
+cov = mne.noise_covariance_segment(raw, reject=dict(eeg=40e-6, eog=150e-6))
print cov
###############################################################################
diff --git a/mne/cov.py b/mne/cov.py
index 2873a5c..24d8a8b 100644
--- a/mne/cov.py
+++ b/mne/cov.py
@@ -3,6 +3,7 @@
#
# License: BSD (3-clause)
+import copy
import os
from math import floor, ceil
import numpy as np
@@ -18,8 +19,9 @@ from .fiff.write import start_block, end_block, write_int, write_name_list, \
write_double, write_float_matrix, start_file, end_file
from .fiff.proj import write_proj, make_projector
from .fiff import fiff_open
-from .fiff.pick import pick_types
-from .epochs import Epochs
+from .fiff.pick import pick_types, channel_indices_by_type
+from .epochs import Epochs, _is_good
+
def rank(A, tol=1e-8):
s = linalg.svd(A, compute_uv=0)
@@ -315,11 +317,9 @@ def read_cov(fid, node, cov_kind):
###############################################################################
# Estimate from data
-def _estimate_noise_covariance_from_epochs(epochs, bmin, bmax,
+def _estimate_noise_covariance_from_epochs(epochs, bmin, bmax, reject, flat,
keep_sample_mean):
"""Estimate noise covariance matrix from epochs
-
- XXX : doc missing
"""
picks_no_eog = pick_types(epochs.info, meg=True, eeg=True, eog=False)
n_channels = len(picks_no_eog)
@@ -331,7 +331,15 @@ def _estimate_noise_covariance_from_epochs(epochs, bmin, bmax,
if bmax is None:
bmax = epochs.times[-1]
bmask = (epochs.times >= bmin) & (epochs.times <= bmax)
+
+ idx_by_type = channel_indices_by_type(epochs.info)
+
for e in epochs:
+
+ if not _is_good(e, epochs.ch_names, idx_by_type, reject, flat):
+ print "Artefact detected in epoch"
+ continue
+
e = e[picks_no_eog]
mu = e[:,bmask].mean(axis=1)
e -= mu[:,None]
@@ -339,21 +347,52 @@ def _estimate_noise_covariance_from_epochs(epochs, bmin, bmax,
e -= np.mean(e, axis=0)
data += np.dot(e, e.T)
n_samples += e.shape[1]
+
print "Number of samples used : %d" % n_samples
print '[done]'
return data, n_samples, ch_names
-def noise_covariance_segment(raw, reject_params=None, keep_sample_mean=True,
- tmin=None, tmax=None, tstep=0.2):
+def noise_covariance_segment(raw, tmin=None, tmax=None, tstep=0.2,
+ reject=None, flat=None, keep_sample_mean=True):
"""Estimate noise covariance matrix from a continuous segment of raw data
- XXX : doc missing
+ It is typically useful to estimate a noise covariance
+ from empty room data or time intervals before starting
+ the stimulation.
Parameters
----------
+ raw : instance of Raw
+ Raw data
+ tmin : float
+ Beginning of time interval in seconds
+ tmax : float
+ End of time interval in seconds
+ tstep : float
+ Size of data chunks for artefact rejection.
+ reject : dict
+ Rejection parameters based on peak to peak amplitude.
+ Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg'.
+ If reject is None then no rejection is done.
+ Values are float. Example:
+ reject = dict(grad=4000e-13, # T / m (gradiometers)
+ mag=4e-12, # T (magnetometers)
+ eeg=40e-6, # uV (EEG channels)
+ eog=250e-6 # uV (EOG channels)
+ )
+ flat : dict
+ Rejection parameters based on flatness of signal
+ Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg'
+ If flat is None then no rejection is done.
+ keep_sample_mean : bool
+ If False data are centered at each instant before computing
+ the covariance.
+
Returns
-------
+ cov : instance of Covariance
+ Noise covariance matrix.
"""
sfreq = raw.info['sfreq']
@@ -370,16 +409,24 @@ def noise_covariance_segment(raw, reject_params=None, keep_sample_mean=True,
n_samples = 0
mu = 0
+ info = copy.copy(raw.info)
+ info['chs'] = [info['chs'][k] for k in picks]
+ info['ch_names'] = [info['ch_names'][k] for k in picks]
+ info['nchan'] = len(picks)
+ idx_by_type = channel_indices_by_type(info)
+
# Read data in chuncks
for first in range(start, stop, step):
last = first + step
if last >= stop:
last = stop
raw_segment, times = raw[picks, first:last]
- # XXX : check for artefacts
- mu += raw_segment[idx].sum(axis=1)
- data += np.dot(raw_segment[idx], raw_segment[idx].T)
- n_samples += raw_segment.shape[1]
+ if _is_good(raw_segment, info['ch_names'], idx_by_type, reject, flat):
+ mu += raw_segment[idx].sum(axis=1)
+ data += np.dot(raw_segment[idx], raw_segment[idx].T)
+ n_samples += raw_segment.shape[1]
+ else:
+ print "Artefact detected in [%d, %d]" % (first, last)
mu /= n_samples
data -= n_samples * mu[:,None] * mu[None,:]
@@ -394,16 +441,53 @@ def noise_covariance_segment(raw, reject_params=None, keep_sample_mean=True,
def noise_covariance(raw, events, event_ids, tmin, tmax,
- bmin=None, bmax=None, reject_params=None,
+ bmin=None, bmax=None, reject=None, flat=None,
keep_sample_mean=True):
- """Estimate noise covariance matrix from raw file
+ """Estimate noise covariance matrix from raw file and events.
- XXX : doc missing
+ The noise covariance is typically estimated on pre-stim periods
+ when the stim onset if defined from events.
Parameters
----------
+ raw : Raw instance
+ The raw data
+ events : array
+ The events a.k.a. the triggers
+ event_ids : array-like of int
+ The valid events to consider
+ tmin : float
+ Initial time in (s) around trigger. Ex: if tmin=0.2
+ then the beginning is 200ms before trigger onset.
+ tmax : float
+ Final time in (s) around trigger. Ex: if tmax=0.5
+ then the end is 500ms after trigger onset.
+ bmin : float
+ Used to specify a specific baseline for the offset.
+ bmin is the init of baseline.
+ bmax : float
+ bmax is the end of baseline.
+ reject : dict
+ Rejection parameters based on peak to peak amplitude.
+ Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg'.
+ If reject is None then no rejection is done.
+ Values are float. Example:
+ reject = dict(grad=4000e-13, # T / m (gradiometers)
+ mag=4e-12, # T (magnetometers)
+ eeg=40e-6, # uV (EEG channels)
+ eog=250e-6 # uV (EOG channels)
+ )
+ flat : dict
+ Rejection parameters based on flatness of signal
+ Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg'
+ If flat is None then no rejection is done.
+ keep_sample_mean : bool
+ If False data are centered at each instant before computing
+ the covariance.
Returns
-------
+ cov : instance of Covariance
+ The computed covariance.
"""
# Pick all channels
picks = pick_types(raw.info, meg=True, eeg=True, eog=True)
@@ -417,12 +501,15 @@ def noise_covariance(raw, events, event_ids, tmin, tmax,
epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
baseline=(None, 0))
d, n, ch_names = _estimate_noise_covariance_from_epochs(epochs,
- bmin=bmin, bmax=bmax, keep_sample_mean=keep_sample_mean)
-
- # XXX : check artefacts
+ bmin=bmin, bmax=bmax, reject=reject, flat=flat,
+ keep_sample_mean=keep_sample_mean)
data += d
n_samples += n
+ if n_samples == 0:
+ raise ValueError('Not enough samples to compute the noise covariance'
+ ' matrix : %d samples' % n_samples)
+
data /= n_samples - 1
cov = Covariance(None)
cov.data = data
diff --git a/mne/epochs.py b/mne/epochs.py
index 916ef57..562ac95 100644
--- a/mne/epochs.py
+++ b/mne/epochs.py
@@ -7,7 +7,7 @@ import copy
import numpy as np
import fiff
from .fiff import Evoked
-from .fiff.pick import channel_type, pick_types
+from .fiff.pick import pick_types, channel_indices_by_type
class Epochs(object):
@@ -223,49 +223,16 @@ class Epochs(object):
n_reject = 0
for k in range(n_events):
e = self._get_epoch_from_disk(k)
- if self._is_good_epoch(e):
+ if ((self.reject is not None or self.flat is not None) and
+ not _is_good(e, self.ch_names, self._channel_type_idx,
+ self.reject, self.flat)):
+ n_reject += 1
+ else:
data[cnt] = self._get_epoch_from_disk(k)
cnt += 1
- else:
- n_reject += 1
print "Rejecting %d epochs." % n_reject
return data[:cnt]
- def _is_good_epoch(self, e):
- """Test is epoch e is good
- """
- if self.reject is not None:
- for key, thresh in self.reject.iteritems():
- idx = self._channel_type_idx[key]
- name = key.upper()
- if len(idx) > 0:
- e_idx = e[idx]
- deltas = np.max(e_idx, axis=1) - np.min(e_idx, axis=1)
- idx_max_delta = np.argmax(deltas)
- delta = deltas[idx_max_delta]
- if delta > thresh:
- ch_name = self.ch_names[idx[idx_max_delta]]
- print '\tRejecting epoch based on %s : %s (%s > %s).' \
- % (name, ch_name, delta, thresh)
- return False
- if self.flat is not None:
- for key, thresh in self.flat.iteritems():
- idx = self._channel_type_idx[key]
- name = key.upper()
- if len(idx) > 0:
- e_idx = e[idx]
- deltas = np.max(e_idx, axis=1) - np.min(e_idx, axis=1)
- idx_min_delta = np.argmin(deltas)
- delta = deltas[idx_min_delta]
- if delta < thresh:
- ch_name = self.ch_names[idx[idx_min_delta]]
- print '\tRejecting flat epoch based on %s : %s (%s < %s).' \
- % (name, ch_name, delta, thresh)
- return False
-
- return True
-
-
def get_data(self):
"""Get all epochs as a 3D array
@@ -285,11 +252,7 @@ class Epochs(object):
if self.reject is None and self.flat is None:
return
- idx = dict(grad=[], mag=[], eeg=[], eog=[], ecg=[])
- for k, ch in enumerate(self.ch_names):
- for key in idx.keys():
- if channel_type(self.info, k) == key:
- idx[key].append(k)
+ idx = channel_indices_by_type(self.info)
for key in idx.keys():
if (self.reject is not None and key in self.reject) \
@@ -347,7 +310,7 @@ class Epochs(object):
evoked.data = data
evoked.times = self.times.copy()
evoked.comment = self.name
- evoked.aspect_kind = np.array([100]) # XXX
+ evoked.aspect_kind = np.array([100]) # for standard average file
evoked.nave = n_events
evoked.first = - np.sum(self.times < 0)
evoked.last = np.sum(self.times > 0)
@@ -365,3 +328,39 @@ class Epochs(object):
evoked.info['nchan'] = len(data_picks)
evoked.data = evoked.data[data_picks]
return evoked
+
+
+def _is_good(e, ch_names, channel_type_idx, reject, flat):
+ """Test if data segment e is good according to the criteria
+ defined in reject and flat.
+ """
+ if reject is not None:
+ for key, thresh in reject.iteritems():
+ idx = channel_type_idx[key]
+ name = key.upper()
+ if len(idx) > 0:
+ e_idx = e[idx]
+ deltas = np.max(e_idx, axis=1) - np.min(e_idx, axis=1)
+ idx_max_delta = np.argmax(deltas)
+ delta = deltas[idx_max_delta]
+ if delta > thresh:
+ ch_name = ch_names[idx[idx_max_delta]]
+ print '\tRejecting epoch based on %s : %s (%s > %s).' \
+ % (name, ch_name, delta, thresh)
+ return False
+ if flat is not None:
+ for key, thresh in flat.iteritems():
+ idx = channel_type_idx[key]
+ name = key.upper()
+ if len(idx) > 0:
+ e_idx = e[idx]
+ deltas = np.max(e_idx, axis=1) - np.min(e_idx, axis=1)
+ idx_min_delta = np.argmin(deltas)
+ delta = deltas[idx_min_delta]
+ if delta < thresh:
+ ch_name = ch_names[idx[idx_min_delta]]
+ print '\tRejecting flat epoch based on %s : %s (%s < %s).' \
+ % (name, ch_name, delta, thresh)
+ return False
+
+ return True
diff --git a/mne/fiff/pick.py b/mne/fiff/pick.py
index 31f2ef0..22a7dc8 100644
--- a/mne/fiff/pick.py
+++ b/mne/fiff/pick.py
@@ -261,3 +261,14 @@ def pick_channels_forward(orig, include=[], exclude=[]):
for k in sel]
return fwd
+
+def channel_indices_by_type(info):
+ """Get indices of channels by type
+ """
+ idx = dict(grad=[], mag=[], eeg=[], eog=[], ecg=[])
+ for k, ch in enumerate(info['chs']):
+ for key in idx.keys():
+ if channel_type(info, k) == key:
+ idx[key].append(k)
+
+ return idx
diff --git a/mne/fiff/tests/data/test-cov.fif b/mne/fiff/tests/data/test-cov.fif
index 574ac2a..b84054b 100644
Binary files a/mne/fiff/tests/data/test-cov.fif and b/mne/fiff/tests/data/test-cov.fif differ
diff --git a/mne/tests/test_cov.py b/mne/tests/test_cov.py
index 6ba2e31..a833db8 100644
--- a/mne/tests/test_cov.py
+++ b/mne/tests/test_cov.py
@@ -39,7 +39,7 @@ def test_cov_estimation_on_raw_segment():
cov = mne.noise_covariance_segment(raw)
cov_mne = mne.Covariance(erm_cov_fname)
assert cov_mne.ch_names == cov.ch_names
- assert (linalg.norm(cov.data - cov_mne.data, ord='fro')
+ assert (linalg.norm(cov.data - cov_mne.data, ord='fro')
/ linalg.norm(cov.data, ord='fro')) < 1e-6
@@ -49,13 +49,14 @@ def test_cov_estimation_with_triggers():
raw = Raw(raw_fname)
events = mne.find_events(raw)
event_ids = [1, 2, 3, 4]
- cov = mne.noise_covariance(raw, events, event_ids,
- tmin=-0.2, tmax=0, keep_sample_mean=True)
+ cov = mne.noise_covariance(raw, events, event_ids, tmin=-0.2, tmax=0,
+ reject=dict(grad=10000e-13, mag=4e-12,
+ eeg=80e-6, eog=150e-6),
+ keep_sample_mean=True)
cov_mne = mne.Covariance(cov_fname)
assert cov_mne.ch_names == cov.ch_names
- # XXX : check something
- # assert (linalg.norm(cov.data - cov_mne.data, ord='fro')
- # / linalg.norm(cov.data, ord='fro')) < 0.1 # XXX : fix
+ assert (linalg.norm(cov.data - cov_mne.data, ord='fro')
+ / linalg.norm(cov.data, ord='fro')) < 0.05
def test_whitening_cov():
diff --git a/mne/time_frequency/tests/test_tfr.py b/mne/time_frequency/tests/test_tfr.py
index 145d8fb..d0c5594 100644
--- a/mne/time_frequency/tests/test_tfr.py
+++ b/mne/time_frequency/tests/test_tfr.py
@@ -32,8 +32,8 @@ def test_time_frequency():
stim=False, include=include, exclude=exclude)
picks = picks[:2]
- epochs = mne.Epochs(raw, events, event_id,
- tmin, tmax, picks=picks, baseline=(None, 0))
+ epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+ baseline=(None, 0))
data = epochs.get_data()
times = epochs.times
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git
More information about the debian-med-commit
mailing list