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 <[email protected]> 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
