diff --git a/doc/source/python_reference.rst b/doc/source/python_reference.rst index 6b4a9fb45dc..3df27b8e3c0 100644 --- a/doc/source/python_reference.rst +++ b/doc/source/python_reference.rst @@ -358,6 +358,17 @@ Manipulate channels and set sensors locations for processing and plotting: run_ica infomax +EEG referencing: + +.. currentmodule:: mne.io + +.. autosummary:: + :toctree: generated/ + :template: function.rst + + set_bipolar_reference + set_eeg_reference + :py:mod:`mne.filter`: .. automodule:: mne.filter diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst index fdb00d27e9d..862203ee024 100644 --- a/doc/source/whats_new.rst +++ b/doc/source/whats_new.rst @@ -49,6 +49,7 @@ Changelog - Add parameter to `whiten_evoked`, `compute_whitener` and `prepare_noise_cov` to set the exact rank by `Martin Luessi`_ and `Denis Engemann`_ + - Add support for bipolar referencing by `Marijn van Vliet`_ BUG ~~~ @@ -80,6 +81,8 @@ BUG - Fix bug in `mne.viz.plot_topo` if ylim was passed for single sensor layouts by `Denis Engemann`_ + - Average reference projections will no longer by automatically added after applying a custom EEG reference by `Marijn van Vliet`_ + API ~~~ diff --git a/examples/plot_rereference_eeg.py b/examples/plot_rereference_eeg.py new file mode 100644 index 00000000000..b9e18d3dd6d --- /dev/null +++ b/examples/plot_rereference_eeg.py @@ -0,0 +1,61 @@ +""" +============================= +Re-referencing the EEG signal +============================= + +Load raw data and apply some EEG referencing schemes. +""" +# Author: Marijn van Vliet +# +# License: BSD (3-clause) + +import mne +from mne.datasets import sample +from matplotlib import pyplot as plt + +# Setup for reading the raw data +data_path = sample.data_path() +raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif' +event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif' +event_id, tmin, tmax = 1, -0.2, 0.5 + +# Read the raw data +raw = mne.io.Raw(raw_fname, preload=True) +events = mne.read_events(event_fname) + +# The EEG channels will be plotted to visualize the difference in referencing +# schemes. +picks = mne.pick_types(raw.info, meg=False, eeg=True, eog=True, exclude='bads') + +############################################################################### +# Apply different EEG referencing schemes and plot the resulting evokeds + +reject = dict(eeg=180e-6, eog=150e-6) +epochs_params = dict(events=events, event_id=event_id, tmin=tmin, tmax=tmax, + picks=picks) + +# No reference. This assumes that the EEG has already been referenced properly. +# This explicitly prevents MNE from adding a default EEG reference. +raw_no_ref, _ = mne.io.set_eeg_reference(raw, []) +evoved_no_ref = mne.Epochs(raw_no_ref, **epochs_params).average() +del raw_no_ref # Free memory + +fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1, sharex=True) +evoved_no_ref.plot(axes=ax1, titles=dict(eeg='EEG Original reference')) + +# Average reference. This is normally added by default, but can also be added +# explicitly. +raw_car, _ = mne.io.set_eeg_reference(raw) +evoked_car = mne.Epochs(raw_car, **epochs_params).average() +del raw_car + +evoked_car.plot(axes=ax2, titles=dict(eeg='EEG Average reference')) + +# Use the mean of channels EEG 001 and EEG 002 as a reference +raw_custom, _ = mne.io.set_eeg_reference(raw, ['EEG 001', 'EEG 002']) +evoked_custom = mne.Epochs(raw_custom, **epochs_params).average() +del raw_custom + +evoked_custom.plot(axes=ax3, titles=dict(eeg='EEG Custom reference')) + +mne.viz.tight_layout() diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index 84aabff85d4..9fd7a726e82 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -14,7 +14,7 @@ from ..utils import logger, verbose from ..io.pick import pick_types from ..forward import _subject_from_forward -from ..minimum_norm.inverse import combine_xyz +from ..minimum_norm.inverse import combine_xyz, _check_reference from ..source_estimate import SourceEstimate from ..time_frequency import CrossSpectralDensity, compute_epochs_csd from ._lcmv import _prepare_beamformer_input @@ -186,6 +186,7 @@ def dics(evoked, forward, noise_csd, data_csd, reg=0.01, label=None, Gross et al. Dynamic imaging of coherent sources: Studying neural interactions in the human brain. PNAS (2001) vol. 98 (2) pp. 694-699 """ + _check_reference(evoked) info = evoked.info data = evoked.data tmin = evoked.times[0] @@ -244,6 +245,7 @@ def dics_epochs(epochs, forward, noise_csd, data_csd, reg=0.01, label=None, Gross et al. Dynamic imaging of coherent sources: Studying neural interactions in the human brain. PNAS (2001) vol. 98 (2) pp. 694-699 """ + _check_reference(epochs) info = epochs.info tmin = epochs.times[0] @@ -481,6 +483,7 @@ def tf_dics(epochs, forward, noise_csds, tmin, tmax, tstep, win_lengths, NOTE : Dalal et al. used a synthetic aperture magnetometry beamformer (SAM) in each time-frequency window instead of DICS. """ + _check_reference(epochs) if pick_ori not in [None, 'normal']: raise ValueError('Unrecognized orientation option in pick_ori, ' diff --git a/mne/beamformer/_lcmv.py b/mne/beamformer/_lcmv.py index 878645af6d2..c9bf0e9a716 100644 --- a/mne/beamformer/_lcmv.py +++ b/mne/beamformer/_lcmv.py @@ -15,7 +15,7 @@ from ..io.proj import make_projector from ..io.pick import pick_types, pick_channels_forward, pick_channels_cov from ..forward import _subject_from_forward -from ..minimum_norm.inverse import _get_vertno, combine_xyz +from ..minimum_norm.inverse import _get_vertno, combine_xyz, _check_reference from ..cov import compute_whitener, compute_covariance from ..source_estimate import _make_stc, SourceEstimate from ..source_space import label_src_vertno_sel @@ -285,6 +285,7 @@ def lcmv(evoked, forward, noise_cov, data_cov, reg=0.01, label=None, beamformers for neuromagnetic source reconstruction. Biomedical Engineering (2004) vol. 51 (10) pp. 1726--34 """ + _check_reference(evoked) info = evoked.info data = evoked.data @@ -348,6 +349,7 @@ def lcmv_epochs(epochs, forward, noise_cov, data_cov, reg=0.01, label=None, beamformers for neuromagnetic source reconstruction. Biomedical Engineering (2004) vol. 51 (10) pp. 1726--34 """ + _check_reference(epochs) info = epochs.info tmin = epochs.times[0] @@ -422,6 +424,7 @@ def lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=None, beamformers for neuromagnetic source reconstruction. Biomedical Engineering (2004) vol. 51 (10) pp. 1726--34 """ + _check_reference(raw) info = raw.info @@ -610,6 +613,7 @@ def tf_lcmv(epochs, forward, noise_covs, tmin, tmax, tstep, win_lengths, time-frequency dynamics of cortical activity. NeuroImage (2008) vol. 40 (4) pp. 1686-1700 """ + _check_reference(epochs) if pick_ori not in [None, 'normal']: raise ValueError('Unrecognized orientation option in pick_ori, ' diff --git a/mne/forward/forward.py b/mne/forward/forward.py index 68e6d409d32..81566160f7b 100644 --- a/mne/forward/forward.py +++ b/mne/forward/forward.py @@ -346,6 +346,11 @@ def _read_forward_meas_info(tree, fid): raise ValueError('MEG/head coordinate transformation not found') info['bads'] = read_bad_channels(fid, parent_meg) + + # Check if a custom reference has been applied + tag = find_tag(fid, parent_mri, FIFF.FIFF_CUSTOM_REF) + info['custom_ref_applied'] = bool(tag.data) if tag is not None else False + return info diff --git a/mne/inverse_sparse/_gamma_map.py b/mne/inverse_sparse/_gamma_map.py index 01d11e58e7f..eab2e220615 100644 --- a/mne/inverse_sparse/_gamma_map.py +++ b/mne/inverse_sparse/_gamma_map.py @@ -8,7 +8,7 @@ from ..forward import is_fixed_orient, _to_fixed_ori from ..io.pick import pick_channels_evoked -from ..minimum_norm.inverse import _prepare_forward +from ..minimum_norm.inverse import _prepare_forward, _check_reference from ..utils import logger, verbose from .mxne_inverse import _make_sparse_stc, _prepare_gain @@ -231,6 +231,8 @@ def gamma_map(evoked, forward, noise_cov, alpha, loose=0.2, depth=0.8, Wipf et al. A unified Bayesian framework for MEG/EEG source imaging, NeuroImage, vol. 44, no. 3, pp. 947-66, Mar. 2009. """ + _check_reference(evoked) + # make forward solution in fixed orientation if necessary if loose is None and not is_fixed_orient(forward): forward = deepcopy(forward) diff --git a/mne/inverse_sparse/mxne_inverse.py b/mne/inverse_sparse/mxne_inverse.py index 7223096c161..bb9dc7a79e9 100644 --- a/mne/inverse_sparse/mxne_inverse.py +++ b/mne/inverse_sparse/mxne_inverse.py @@ -8,6 +8,7 @@ from ..source_estimate import SourceEstimate from ..minimum_norm.inverse import combine_xyz, _prepare_forward +from ..minimum_norm.inverse import _check_reference from ..forward import compute_orient_prior, is_fixed_orient, _to_fixed_ori from ..io.pick import pick_channels_evoked from .mxne_optim import mixed_norm_solver, norm_l2inf, tf_mixed_norm_solver @@ -159,6 +160,8 @@ def mixed_norm(evoked, forward, noise_cov, alpha, loose=0.2, depth=0.8, if not isinstance(evoked, list): evoked = [evoked] + _check_reference(evoked[0]) + all_ch_names = evoked[0].ch_names if not all(all_ch_names == evoked[i].ch_names for i in range(1, len(evoked))): @@ -361,6 +364,8 @@ def tf_mixed_norm(evoked, forward, noise_cov, alpha_space, alpha_time, The residual a.k.a. data not explained by the sources. Only returned if return_residual is True. """ + _check_reference(evoked) + all_ch_names = evoked.ch_names info = evoked.info diff --git a/mne/io/__init__.py b/mne/io/__init__.py index 0d74a626a65..0dedf4d78f5 100644 --- a/mne/io/__init__.py +++ b/mne/io/__init__.py @@ -30,4 +30,5 @@ # for backward compatibility from .fiff import RawFIFF from .fiff import RawFIFF as Raw -from .base import concatenate_raws, get_chpi_positions, set_eeg_reference +from .base import concatenate_raws, get_chpi_positions +from .reference import set_eeg_reference, set_bipolar_reference diff --git a/mne/io/base.py b/mne/io/base.py index 885343dea37..f5fde7cea14 100644 --- a/mne/io/base.py +++ b/mne/io/base.py @@ -3,6 +3,7 @@ # Martin Luessi # Denis Engemann # Teon Brooks +# Marijn van Vliet # # License: BSD (3-clause) @@ -20,8 +21,7 @@ from .constants import FIFF from .pick import pick_types, channel_type, pick_channels from .meas_info import write_meas_info -from .proj import (setup_proj, activate_proj, proj_equal, ProjMixin, - _has_eeg_average_ref_proj, make_eeg_average_ref_proj) +from .proj import setup_proj, activate_proj, proj_equal, ProjMixin from ..channels.channels import (ContainsMixin, PickDropChannelsMixin, SetChannelsMixin, InterpolationMixin) from ..channels.layout import read_montage, apply_montage, Montage @@ -77,15 +77,6 @@ def __hash__(self): raise RuntimeError('Cannot hash raw unless preloaded') return object_hash(dict(info=self.info, data=self._data)) - def _add_eeg_ref(self, add_eeg_ref): - """Helper to add an average EEG reference""" - if add_eeg_ref: - eegs = pick_types(self.info, meg=False, eeg=True, ref_meg=False) - projs = self.info['projs'] - if len(eegs) > 0 and not _has_eeg_average_ref_proj(projs): - eeg_ref = make_eeg_average_ref_proj(self.info, activate=False) - projs.append(eeg_ref) - def _parse_get_set_params(self, item): # make sure item is a tuple if not isinstance(item, tuple): # only channel selection passed @@ -397,7 +388,8 @@ def filter(self, l_freq, h_freq, picks=None, filter_length='10s', if l_freq < h_freq: logger.info('Band-pass filtering from %0.2g - %0.2g Hz' % (l_freq, h_freq)) - self._data = band_pass_filter(self._data, fs, l_freq, h_freq, + self._data = band_pass_filter( + self._data, fs, l_freq, h_freq, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, h_trans_bandwidth=h_trans_bandwidth, @@ -406,7 +398,8 @@ def filter(self, l_freq, h_freq, picks=None, filter_length='10s', else: logger.info('Band-stop filtering from %0.2g - %0.2g Hz' % (h_freq, l_freq)) - self._data = band_stop_filter(self._data, fs, h_freq, l_freq, + self._data = band_stop_filter( + self._data, fs, h_freq, l_freq, filter_length=filter_length, l_trans_bandwidth=h_trans_bandwidth, h_trans_bandwidth=l_trans_bandwidth, method=method, @@ -741,7 +734,7 @@ def save(self, fname, picks=None, tmin=0, tmax=None, buffer_size_sec=10, int=FIFF.FIFFT_INT, single=FIFF.FIFFT_FLOAT, double=FIFF.FIFFT_DOUBLE) - if not format in type_dict.keys(): + if format not in type_dict.keys(): raise ValueError('format must be "short", "int", "single", ' 'or "double"') reset_dict = dict(short=False, int=False, single=True, double=True) @@ -1330,60 +1323,6 @@ def add_events(self, events, stim_channel=None): self._data[pick, idx - self.first_samp] += events[:, 2] -def set_eeg_reference(raw, ref_channels, copy=True): - """Rereference eeg channels to new reference channel(s). - - If multiple reference channels are specified, they will be averaged. - - Parameters - ---------- - raw : instance of Raw - Instance of Raw with eeg channels and reference channel(s). - - ref_channels : list of str - The name(s) of the reference channel(s). - - copy : bool - Specifies whether instance of Raw will be copied or modified in place. - - Returns - ------- - raw : instance of Raw - Instance of Raw with eeg channels rereferenced. - - ref_data : array - Array of reference data subtracted from eeg channels. - """ - # Check to see that raw data is preloaded - if not raw.preload: - raise RuntimeError('Raw data needs to be preloaded. Use ' - 'preload=True (or string) in the constructor.') - # Make sure that reference channels are loaded as list of string - if not isinstance(ref_channels, list): - raise IOError('Reference channel(s) must be a list of string. ' - 'If using a single reference channel, enter as ' - 'a list with one element.') - # Find the indices to the reference electrodes - ref_idx = [raw.ch_names.index(c) for c in ref_channels] - - # Get the reference array - ref_data = raw._data[ref_idx].mean(0) - - # Get the indices to the eeg channels using the pick_types function - eeg_idx = pick_types(raw.info, exclude="bads", eeg=True, meg=False, - ref_meg=False) - - # Copy raw data or modify raw data in place - if copy: # copy data - raw = raw.copy() - - # Rereference the eeg channels - raw._data[eeg_idx] -= ref_data - - # Return rereferenced data and reference array - return raw, ref_data - - def _allocate_data(data, data_buffer, data_shape, dtype): if data is None: # if not already done, allocate array with right type @@ -1526,7 +1465,8 @@ def _write_raw(fname, raw, info, picks, format, data_type, reset_range, start, # Split files if necessary, leave some space for next file info if pos >= split_size - this_buff_size_bytes - 2 ** 20: - next_fname, next_idx = _write_raw(fname, raw, info, picks, format, + next_fname, next_idx = _write_raw( + fname, raw, info, picks, format, data_type, reset_range, first + buffer_size, stop, buffer_size, projector, inv_comp, drop_small_buffer, split_size, part_idx + 1, use_fname) @@ -1646,7 +1586,7 @@ def _write_raw_buffer(fid, buf, cals, format, inv_comp): if buf.shape[0] != len(cals): raise ValueError('buffer and calibration sizes do not match') - if not format in ['short', 'int', 'single', 'double']: + if format not in ['short', 'int', 'single', 'double']: raise ValueError('format must be "short", "single", or "double"') if np.isrealobj(buf): @@ -1864,7 +1804,7 @@ def _check_update_montage(info, montage): # raise error if positions are missing if missing_positions: err = ("The following positions are missing from the montage " - "definitions: %s. If those channels lack positions " - "because they are EOG channels use the eog parameter." - % str(missing_positions)) + "definitions: %s. If those channels lack positions " + "because they are EOG channels use the eog parameter." + % str(missing_positions)) raise KeyError(err) diff --git a/mne/io/brainvision/brainvision.py b/mne/io/brainvision/brainvision.py index b1857e14355..1bf3a870440 100644 --- a/mne/io/brainvision/brainvision.py +++ b/mne/io/brainvision/brainvision.py @@ -436,6 +436,7 @@ def _get_eeg_info(vhdr_fname, reference, eog, misc): info['orig_blocks'] = None info['line_freq'] = None info['subject_info'] = None + info['custom_ref_applied'] = False eeg_info = {} diff --git a/mne/io/bti/bti.py b/mne/io/bti/bti.py index df528939167..7b1db0bf2b6 100644 --- a/mne/io/bti/bti.py +++ b/mne/io/bti/bti.py @@ -781,22 +781,22 @@ def _read_bti_header(pdf_fname, config_fname): fid.seek(1, 1) info.update({'data_format': read_int16(fid), - 'acq_mode': read_int16(fid), - 'total_epochs': read_int32(fid), - 'input_epochs': read_int32(fid), - 'total_events': read_int32(fid), - 'total_fixed_events': read_int32(fid), - 'sample_period': read_float(fid), - 'xaxis_label': read_str(fid, 16), - 'total_processes': read_int32(fid), - 'total_chans': read_int16(fid)}) + 'acq_mode': read_int16(fid), + 'total_epochs': read_int32(fid), + 'input_epochs': read_int32(fid), + 'total_events': read_int32(fid), + 'total_fixed_events': read_int32(fid), + 'sample_period': read_float(fid), + 'xaxis_label': read_str(fid, 16), + 'total_processes': read_int32(fid), + 'total_chans': read_int16(fid)}) fid.seek(2, 1) info.update({'checksum': read_int32(fid), - 'total_ed_classes': read_int32(fid), - 'total_associated_files': read_int16(fid), - 'last_file_index': read_int16(fid), - 'timestamp': read_int32(fid)}) + 'total_ed_classes': read_int32(fid), + 'total_associated_files': read_int16(fid), + 'last_file_index': read_int16(fid), + 'timestamp': read_int32(fid)}) fid.seek(20, 1) _correct_offset(fid) @@ -988,7 +988,7 @@ def __init__(self, pdf_fname, config_fname='config', logger.info('Reading 4D PDF file %s...' % pdf_fname) bti_info = _read_bti_header(pdf_fname, config_fname) - # XXX indx is informed guess. Normally only one transform is stored. + # XXX indx is informed guess. Normally only one transform is stored. dev_ctf_t = bti_info['bti_transform'][0].astype('>f8') bti_to_nm = bti_to_vv_trans(adjust=rotation_x, translation=translation, dtype='>f8') @@ -1022,6 +1022,7 @@ def __init__(self, pdf_fname, config_fname='config', info['lowpass'] = lp info['acq_pars'], info['acq_stim'] = None, None info['filename'] = None + info['custom_ref_applied'] = False chs = [] ch_names = [ch['name'] for ch in bti_info['chs']] diff --git a/mne/io/constants.py b/mne/io/constants.py index caad8b7170c..224df3d9377 100644 --- a/mne/io/constants.py +++ b/mne/io/constants.py @@ -124,6 +124,7 @@ def __init__(self, **kwargs): FIFF.FIFF_DESCRIPTION = FIFF.FIFF_COMMENT # (Textual) Description of an object FIFF.FIFF_DIG_STRING = 234 # String of digitized points FIFF.FIFF_LINE_FREQ = 235 # Line frequency +FIFF.FIFF_CUSTOM_REF = 236 # Whether a custom reference was applied to the data # # HPI fitting program tags # diff --git a/mne/io/edf/edf.py b/mne/io/edf/edf.py index a3e535b691a..6862b64de2a 100644 --- a/mne/io/edf/edf.py +++ b/mne/io/edf/edf.py @@ -13,7 +13,6 @@ import re import warnings from math import ceil, floor -import warnings import numpy as np from scipy.interpolate import interp1d @@ -422,6 +421,7 @@ def _get_edf_info(fname, stim_channel, annot, annotmap, tal_channel, info['experimenter'] = None info['line_freq'] = None info['subject_info'] = None + info['custom_ref_applied'] = False edf_info = dict() edf_info['annot'] = annot diff --git a/mne/io/egi/egi.py b/mne/io/egi/egi.py index 1936b22ea82..de9bbcbe2bb 100644 --- a/mne/io/egi/egi.py +++ b/mne/io/egi/egi.py @@ -282,6 +282,7 @@ def __init__(self, input_fname, montage=None, eog=None, misc=None, info['ch_names'] = ch_names info['bads'] = [] info['comps'] = [] + info['custom_ref_applied'] = False for ii, ch_name in enumerate(ch_names): ch_info = {'cal': cal, 'logno': ii + 1, diff --git a/mne/io/fiff/raw.py b/mne/io/fiff/raw.py index 0ed334860fb..c1479628f08 100644 --- a/mne/io/fiff/raw.py +++ b/mne/io/fiff/raw.py @@ -17,7 +17,8 @@ from ..meas_info import read_meas_info from ..tree import dir_tree_find from ..tag import read_tag -from ..proj import proj_equal +from ..proj import (proj_equal, make_eeg_average_ref_proj, + _needs_eeg_average_ref_proj) from ..compensator import get_current_comp, set_current_comp, make_compensator from ..base import _BaseRaw @@ -127,7 +128,10 @@ def __init__(self, fnames, allow_maxshield=False, preload=False, self.verbose = verbose self.orig_format = raws[0].orig_format self.proj = False - self._add_eeg_ref(add_eeg_ref) + + if add_eeg_ref and _needs_eeg_average_ref_proj(self.info): + eeg_ref = make_eeg_average_ref_proj(self.info, activate=False) + self.add_proj(eeg_ref) if preload: self._preload_data(preload) @@ -305,8 +309,9 @@ def _read_raw_file(self, fname, allow_maxshield, preload, compensation, idx2 = base.rfind('-') if idx2 < 0 and next_num == 1: # this is the first file, which may not be numbered - next_fname = op.join(path, '%s-%d.%s' % (base[:idx], - next_num, base[idx + 1:])) + next_fname = op.join( + path, '%s-%d.%s' % (base[:idx], next_num, + base[idx + 1:])) continue num_str = base[idx2 + 1:idx] if not num_str.isdigit(): diff --git a/mne/io/fiff/tests/test_raw.py b/mne/io/fiff/tests/test_raw.py index c67fe043d8d..408a7a7ddc3 100644 --- a/mne/io/fiff/tests/test_raw.py +++ b/mne/io/fiff/tests/test_raw.py @@ -20,8 +20,7 @@ from mne import pick_types, pick_channels from mne.datasets import testing from mne.io.constants import FIFF -from mne.io import (Raw, concatenate_raws, - get_chpi_positions, set_eeg_reference) +from mne.io import Raw, concatenate_raws, get_chpi_positions from mne import concatenate_events, find_events, equalize_channels from mne.utils import (_TempDir, requires_nitime, requires_pandas, requires_mne, run_subprocess, run_tests_if_main) @@ -982,39 +981,6 @@ def compensate_mne(fname, grad): assert_allclose(raw_py._data, raw_c._data, rtol=1e-6, atol=1e-17) -@testing.requires_testing_data -def test_set_eeg_reference(): - """ Test rereference eeg data""" - raw = Raw(fif_fname, preload=True) - - # Rereference raw data by creating a copy of original data - reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], copy=True) - - # Separate EEG channels from other channel types - picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') - picks_other = pick_types(raw.info, meg=True, eeg=False, eog=True, - stim=True, exclude='bads') - - # Get the raw EEG data and other channel data - raw_eeg_data = raw[picks_eeg][0] - raw_other_data = raw[picks_other][0] - - # Get the rereferenced EEG data and channel other - reref_eeg_data = reref[picks_eeg][0] - unref_eeg_data = reref_eeg_data + ref_data - # Undo rereferencing of EEG channels - reref_other_data = reref[picks_other][0] - - # Check that both EEG data and other data is the same - assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) - assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) - - # Test that data is modified in place when copy=False - reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], - copy=False) - assert_true(raw is reref) - - @testing.requires_testing_data def test_drop_channels_mixin(): """Test channels-dropping functionality diff --git a/mne/io/kit/kit.py b/mne/io/kit/kit.py index 9adb75d1d99..1cf3f353747 100644 --- a/mne/io/kit/kit.py +++ b/mne/io/kit/kit.py @@ -116,6 +116,7 @@ def __init__(self, input_fname, mrk=None, elp=None, hsp=None, stim='>', self.cals = np.ones(self.info['nchan']) self.orig_format = 'double' self.rawdirs = list() + self.info['custom_ref_applied'] = False if isinstance(mrk, list): mrk = [read_mrk(marker) if isinstance(marker, string_types) diff --git a/mne/io/meas_info.py b/mne/io/meas_info.py index 15853711744..d2010d4855a 100644 --- a/mne/io/meas_info.py +++ b/mne/io/meas_info.py @@ -405,6 +405,7 @@ def read_meas_info(fid, tree, verbose=None): proj_id = None proj_name = None line_freq = None + custom_ref_applied = False p = 0 for k in range(meas_info['nent']): kind = meas_info['directory'][k].kind @@ -452,6 +453,9 @@ def read_meas_info(fid, tree, verbose=None): elif kind == FIFF.FIFF_LINE_FREQ: tag = read_tag(fid, pos) line_freq = float(tag.data) + elif kind == FIFF.FIFF_CUSTOM_REF: + tag = read_tag(fid, pos) + custom_ref_applied = bool(tag.data) # Check that we have everything we need if nchan is None: @@ -627,6 +631,7 @@ def read_meas_info(fid, tree, verbose=None): info['comps'] = comps info['acq_pars'] = acq_pars info['acq_stim'] = acq_stim + info['custom_ref_applied'] = custom_ref_applied return info, meas @@ -745,6 +750,8 @@ def write_meas_info(fid, info, data_type=None, reset_range=True): write_float(fid, FIFF.FIFF_LINE_FREQ, info['line_freq']) if data_type is not None: write_int(fid, FIFF.FIFF_DATA_PACK, data_type) + if info.get('custom_ref_applied'): + write_int(fid, FIFF.FIFF_CUSTOM_REF, info['custom_ref_applied']) # Channel information for k, c in enumerate(info['chs']): @@ -927,10 +934,11 @@ def _merge_info(infos, verbose=None): trans_name) raise ValueError(msg) other_fields = ['acq_pars', 'acq_stim', 'bads', 'buffer_size_sec', - 'comps', 'description', 'dig', 'experimenter', 'file_id', - 'filename', 'highpass', 'line_freq', 'lowpass', - 'meas_date', 'meas_id', 'orig_blocks', 'proj_id', - 'proj_name', 'projs', 'sfreq', 'subject_info', 'sfreq'] + 'comps', 'custom_ref_applied', 'description', 'dig', + 'experimenter', 'file_id', 'filename', 'highpass', + 'line_freq', 'lowpass', 'meas_date', 'meas_id', + 'orig_blocks', 'proj_id', 'proj_name', 'projs', 'sfreq', + 'subject_info', 'sfreq'] for k in other_fields: info[k] = _merge_dict_values(infos, k) @@ -997,4 +1005,5 @@ def create_info(ch_names, sfreq, ch_types=None): info['dev_head_t'] = None info['dev_ctf_t'] = None info['ctf_head_t'] = None + info['custom_ref_applied'] = False return info diff --git a/mne/io/proj.py b/mne/io/proj.py index a13f3a2418d..b66c35a093a 100644 --- a/mne/io/proj.py +++ b/mne/io/proj.py @@ -17,6 +17,7 @@ from .constants import FIFF from .pick import pick_types from ..utils import logger, verbose +from ..externals.six import string_type class Projection(dict): @@ -165,10 +166,10 @@ def del_proj(self, idx): self.info['projs'].pop(idx) return self - + def plot_projs_topomap(self, ch_type=None, layout=None): """Plot SSP vector - + Parameters ---------- ch_type : 'mag' | 'grad' | 'planar1' | 'planar2' | 'eeg' | None | List @@ -196,8 +197,8 @@ def plot_projs_topomap(self, ch_type=None, layout=None): layout = [] if ch_type is None: ch_type = [ch for ch in ['meg', 'eeg'] if ch in self] - elif isinstance(ch_type, str): - ch_type = [ch_type] + elif isinstance(ch_type, string_type): + ch_type = [ch_type] for ch in ch_type: if ch in self: layout.append(find_layout(self.info, ch, exclude=[])) @@ -207,7 +208,7 @@ def plot_projs_topomap(self, ch_type=None, layout=None): fig = plot_projs_topomap(self.info['projs'], layout) else: raise ValueError("Info is missing projs. Nothing to plot.") - + return fig @@ -584,6 +585,11 @@ def make_eeg_average_ref_proj(info, activate=True, verbose=None): eeg_proj: instance of Projection The SSP/PCA projector. """ + if info['custom_ref_applied']: + raise RuntimeError('Cannot add an average EEG reference projection ' + 'since a custom reference has been applied to the ' + 'data earlier.') + logger.info("Adding average EEG reference projection.") eeg_sel = pick_types(info, meg=False, eeg=True, ref_meg=False, exclude='bads') @@ -605,12 +611,25 @@ def make_eeg_average_ref_proj(info, activate=True, verbose=None): def _has_eeg_average_ref_proj(projs): """Determine if a list of projectors has an average EEG ref""" for proj in projs: - if proj['desc'] == 'Average EEG reference' or \ - proj['kind'] == FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF: + if (proj['desc'] == 'Average EEG reference' or + proj['kind'] == FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF): return True return False +def _needs_eeg_average_ref_proj(info): + """Determine if the EEG needs an averge EEG reference + + This returns True if no custom reference has been applied and no average + reference projection is present in the list of projections. + """ + eeg_sel = pick_types(info, meg=False, eeg=True, ref_meg=False, + exclude='bads') + return (len(eeg_sel) > 0 and + not info['custom_ref_applied'] and + not _has_eeg_average_ref_proj(info['projs'])) + + @verbose def setup_proj(info, add_eeg_ref=True, activate=True, verbose=None): @@ -636,14 +655,11 @@ def setup_proj(info, add_eeg_ref=True, activate=True, The modified measurement info (Warning: info is modified inplace). """ # Add EEG ref reference proj if necessary - eeg_sel = pick_types(info, meg=False, eeg=True, ref_meg=False, - exclude='bads') - if len(eeg_sel) > 0 and not _has_eeg_average_ref_proj(info['projs']) \ - and add_eeg_ref is True: + if _needs_eeg_average_ref_proj(info) and add_eeg_ref: eeg_proj = make_eeg_average_ref_proj(info, activate=activate) info['projs'].append(eeg_proj) - # Create the projector + # Create the projector projector, nproj = make_projector_info(info) if nproj == 0: if verbose: @@ -654,7 +670,7 @@ def setup_proj(info, add_eeg_ref=True, activate=True, logger.info('Created an SSP operator (subspace dimension = %d)' % nproj) - # The projection items have been activated + # The projection items have been activated if activate: info['projs'] = activate_proj(info['projs'], copy=False) diff --git a/mne/io/reference.py b/mne/io/reference.py new file mode 100644 index 00000000000..651aeec5255 --- /dev/null +++ b/mne/io/reference.py @@ -0,0 +1,276 @@ +# Authors: Marijn van Vliet +# Alexandre Gramfort +# +# License: BSD (3-clause) + +import numpy as np + +from .constants import FIFF +from .proj import _has_eeg_average_ref_proj, make_eeg_average_ref_proj +from .pick import pick_types +from ..utils import logger + + +def _apply_reference(raw, ref_from, ref_to=None, copy=True): + """Apply a custom EEG referencing scheme. + + Calculates a reference signal by taking the mean of a set of channels and + applies the reference to another set of channels. + + Parameters + ---------- + raw : instance of Raw + Instance of Raw with EEG channels and reference channel(s). + ref_from : list of str + The names of the channels to use to construct the reference. If an + empty list is specified, the data is assumed to already have a proper + reference and MNE will not attempt any re-referencing of the data. + ref_to : list of str | None + The names of the channels to apply the reference to. By default, + all EEG channels are chosen. + copy : bool + Specifies whether instance of Raw will be copied or modified in place. + Defaults to True. + + Returns + ------- + raw : instance of Raw + Instance of Raw with EEG channels rereferenced. + ref_data : array, shape (n_times,) + Array of reference data subtracted from EEG channels. + + Notes + ----- + 1. Do not use this function to apply an average reference. By default, an + average reference projection has already been added upon loading raw + data. + + 2. If the reference is applied to any EEG channels, this function removes + any pre-existing average reference projections. + + 3. During source localization, the EEG signal should have an average + reference. + + 4. The data must be preloaded. + + See also + -------- + set_eeg_reference : Convenience function for creating an EEG reference. + set_bipolar_reference : Convenience function for creating a bipolar + reference. + """ + # Check to see that raw data is preloaded + if not raw.preload: + raise RuntimeError('Raw data needs to be preloaded. Use ' + 'preload=True (or string) in the constructor.') + + eeg_idx = pick_types(raw.info, eeg=True, meg=False, ref_meg=False) + + ref_from = [raw.ch_names.index(ch) for ch in ref_from] + if ref_to is None: + ref_to = eeg_idx + else: + ref_to = [raw.ch_names.index(ch) for ch in ref_to] + + if copy: + raw = raw.copy() + + # Compute reference + if len(ref_from) > 0: + ref_data = raw._data[ref_from].mean(0) + raw._data[ref_to] -= ref_data + else: + ref_data = None + + # If the reference touches EEG electrodes, remove any pre-existing common + # reference and note in the info that a non-CAR has been applied. + if len(np.intersect1d(ref_to, eeg_idx)) > 0: + raw.info['custom_ref_applied'] = True + + # Remove any existing average reference projections + for i, proj in enumerate(raw.info['projs']): + if proj['desc'] == 'Average EEG reference' or \ + proj['kind'] == FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF: + logger.info('Removing existing average EEG reference ' + 'projection.') + del raw.info['projs'][i] + + return raw, ref_data + + +def set_eeg_reference(raw, ref_channels=None, copy=True): + """Rereference EEG channels to new reference channel(s). + + If multiple reference channels are specified, they will be averaged. If + no reference channels are specified, an average reference will be applied. + + Parameters + ---------- + raw : instance of Raw + Instance of Raw with EEG channels and reference channel(s). + ref_channels : list of str | None + The names of the channels to use to construct the reference. If None is + specified here, an average reference will be applied in the form of an + SSP projector. If an empty list is specified, the data is assumed to + already have a proper reference and MNE will not attempt any + re-referencing of the data. + copy : bool + Specifies whether instance of Raw will be copied (True) or modified in + place (False). Defaults to True. + + Returns + ------- + raw : instance of Raw + Instance of Raw with eeg channels rereferenced. + ref_data : array + Array of reference data subtracted from EEG channels. + + Notes + ----- + 1. If a reference is requested that is not the average reference, this + function removes any pre-existing average reference projections. + + 2. During source localization, the EEG signal should have an average + reference. + + 3. In order to apply a reference other than an average reference, the data + must be preloaded. + + See also + -------- + set_bipolar_reference : Convenience function for creating bipolar + references. + """ + if ref_channels is None: + # CAR requested + if _has_eeg_average_ref_proj(raw.info['projs']): + logger.warning('An average reference projection was already ' + 'added. The data has been left untouched.') + return raw, None + else: + raw.add_proj(make_eeg_average_ref_proj(raw.info), activate=False) + return raw, None + else: + logger.info('Applying a custom EEG reference.') + return _apply_reference(raw, ref_channels, copy=copy) + + +def set_bipolar_reference(raw, anode, cathode, ch_name=None, ch_info=None, + copy=True): + """Rereference selected channels using a bipolar referencing scheme. + + A bipolar reference takes the difference between two channels (the anode + minus the cathode) and adds it as a new virtual channel. The original + channels will be dropped. + + Multiple anodes and cathodes can be specified, in which case multiple + vitual channels will be created. The 1st anode will be substracted from the + 1st cathode, the 2nd anode from the 2nd cathode, etc. + + By default, the virtual channels will be annotated with channel info of + the anodes, their locations set to (0, 0, 0) and coil types set to + EEG_BIPOLAR. + + Parameters + ---------- + raw : instance of Raw + Instance of Raw containing the unreferenced channels. + anode : str | list of str + The name(s) of the channel(s) to use as anode in the bipolar reference. + cathode : str | list of str + The name(s) of the channel(s) to use as cathode in the bipolar + reference. + ch_name : str | list of str | None + The channel name(s) for the virtual channel(s) containing the resulting + signal. By default, bipolar channels are named after the anode and + cathode, but it is recommended to supply a more meaningful name. + ch_info : dict | list of dict | None + This parameter can be used to supply a dictionary (or a dictionary for + each bipolar channel) containing channel information to merge in, + overwriting the default values. + copy : bool + Whether to operate on a copy of the raw data (True) or modify it + in-place (False). Defaults to True. + + Returns + ------- + raw : instance of Raw + Instance of Raw with the specified channels rereferenced. + + See also + -------- + set_eeg_reference : Convenience function for creating an EEG reference. + + Notes + ----- + 1. If the anodes contain any EEG channels, this function removes + any pre-existing average reference projections. + + 2. During source localization, the EEG signal should have an average + reference. + + 3. The data must be preloaded. + """ + if not isinstance(anode, list): + anode = [anode] + + if not isinstance(cathode, list): + cathode = [cathode] + + if len(anode) != len(cathode): + raise ValueError('Number of anodes must equal the number of cathodes.') + + if ch_name is None: + ch_name = ['%s-%s' % ac for ac in zip(anode, cathode)] + elif not isinstance(ch_name, list): + ch_name = [ch_name] + if len(ch_name) != len(anode): + raise ValueError('Number of channel names must equal the number of ' + 'anodes/cathodes.') + + # Check for duplicate channel names (it is allowed to give the name of the + # anode or cathode channel, as they will be replaced). + for ch, a, c in zip(ch_name, anode, cathode): + if ch not in [a, c] and ch in raw.ch_names: + raise ValueError('There is already a channel named "%s", please ' + 'specify a different name for the bipolar ' + 'channel using the ch_name parameter.' % ch) + + if ch_info is None: + ch_info = [{} for an in anode] + elif not isinstance(ch_info, list): + ch_info = [ch_info] + if len(ch_info) != len(anode): + raise ValueError('Number of channel info dictionaries must equal the ' + 'number of anodes/cathodes.') + + # Merge specified and anode channel information dictionaries + new_ch_info = [] + for an, ci in zip(anode, ch_info): + new_info = raw.info['chs'][raw.ch_names.index(an)].copy() + + # Set channel location and coil type + if 'eeg_loc' in new_info: + new_info['eeg_loc'] = np.zeros((3, 2)) + new_info['loc'] = np.zeros(12) + new_info['coil_type'] = FIFF.FIFFV_COIL_EEG_BIPOLAR + + new_info.update(ci) + new_ch_info.append(new_info) + + if copy: + raw = raw.copy() + + # Perform bipolar referencing + for an, ca, name, info in zip(anode, cathode, ch_name, new_ch_info): + raw, _ = _apply_reference(raw, [ca], [an], copy=False) + an_idx = raw.ch_names.index(an) + raw.info['chs'][an_idx] = info + raw.info['chs'][an_idx]['ch_name'] = name + raw.info['ch_names'][an_idx] = name + logger.info('Bipolar channel added as "%s".' % name) + + # Drop cathode channels + raw.drop_channels(cathode) + + return raw diff --git a/mne/io/tests/test_reference.py b/mne/io/tests/test_reference.py new file mode 100644 index 00000000000..ad4b0db8b09 --- /dev/null +++ b/mne/io/tests/test_reference.py @@ -0,0 +1,157 @@ +# Authors: Marijn van Vliet +# Alexandre Gramfort +# +# License: BSD (3-clause) + +import os.path as op +import warnings + +from nose.tools import assert_true, assert_equal, assert_raises +from numpy.testing import assert_array_equal, assert_allclose + +from mne.io.constants import FIFF +from mne.io import set_eeg_reference, set_bipolar_reference +from mne.io.proj import _has_eeg_average_ref_proj +from mne.io.reference import _apply_reference +from mne.datasets import testing +from mne.io import Raw +from mne import pick_types + +warnings.simplefilter('always') # enable b/c these tests throw warnings + +data_dir = op.join(testing.data_path(download=False), 'MEG', 'sample') +fif_fname = op.join(data_dir, 'sample_audvis_trunc_raw.fif') + + +@testing.requires_testing_data +def test_apply_reference(): + """ Test base function for rereferencing""" + raw = Raw(fif_fname, preload=True, add_eeg_ref=True) + + # Separate EEG channels from other channel types + picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') + picks_other = pick_types(raw.info, meg=True, eeg=False, eog=True, + stim=True, exclude='bads') + + # Rereference raw data by creating a copy of original data + reref, ref_data = _apply_reference(raw, ref_from=['EEG 001', 'EEG 002'], + copy=True) + assert_true(reref.info['custom_ref_applied']) + + # The CAR reference projection should have been removed by the function + assert_true(not _has_eeg_average_ref_proj(reref.info['projs'])) + + # Check that the ref has been properly computed + ref_ch_idx = [raw.ch_names.index(ch) for ch in ['EEG 001', 'EEG 002']] + assert_array_equal(ref_data, raw[ref_ch_idx, :][0].mean(0)) + + # Get the raw EEG data and other channel data + raw_eeg_data = raw[picks_eeg][0] + raw_other_data = raw[picks_other][0] + + # Get the rereferenced EEG data and channel other + reref_eeg_data = reref[picks_eeg][0] + # Undo rereferencing of EEG channels + unref_eeg_data = reref_eeg_data + ref_data + reref_other_data = reref[picks_other][0] + + # Check that both EEG data and other data is the same + assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) + assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) + + # Test that disabling the reference does not break anything + reref, ref_data = _apply_reference(raw, []) + + # Test that data is modified in place when copy=False + reref, ref_data = _apply_reference(raw, ['EEG 001', 'EEG 002'], + copy=False) + assert_true(raw is reref) + + +@testing.requires_testing_data +def test_set_eeg_reference(): + """ Test rereference eeg data""" + raw = Raw(fif_fname, preload=True) + + # Rereference raw data by creating a copy of original data + reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], copy=True) + assert_true(reref.info['custom_ref_applied']) + + # Separate EEG channels from other channel types + picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') + picks_other = pick_types(raw.info, meg=True, eeg=False, eog=True, + stim=True, exclude='bads') + + # Get the raw EEG data and other channel data + raw_eeg_data = raw[picks_eeg][0] + raw_other_data = raw[picks_other][0] + + # Get the rereferenced EEG data and channel other + reref_eeg_data = reref[picks_eeg][0] + # Undo rereferencing of EEG channels + unref_eeg_data = reref_eeg_data + ref_data + reref_other_data = reref[picks_other][0] + + # Check that both EEG data and other data is the same + assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) + assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) + + # Test that data is modified in place when copy=False + reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], + copy=False) + assert_true(raw is reref) + + +@testing.requires_testing_data +def test_set_bipolar_reference(): + """ Test bipolar referencing""" + raw = Raw(fif_fname, preload=True) + reref = set_bipolar_reference(raw, 'EEG 001', 'EEG 002', 'bipolar', + {'kind': FIFF.FIFFV_EOG_CH, + 'extra': 'some extra value'}) + assert_true(reref.info['custom_ref_applied']) + + # Compare result to a manual calculation + a = raw.pick_channels(['EEG 001', 'EEG 002'], copy=True) + a = a._data[0, :] - a._data[1, :] + b = reref.pick_channels(['bipolar'], copy=True)._data[0, :] + assert_allclose(a, b) + + # Original channels should be replaced by a virtual one + assert_true('EEG 001' not in reref.ch_names) + assert_true('EEG 002' not in reref.ch_names) + assert_true('bipolar' in reref.ch_names) + + # Check channel information + bp_info = reref.info['chs'][reref.ch_names.index('bipolar')] + an_info = reref.info['chs'][raw.ch_names.index('EEG 001')] + for key in bp_info: + if key == 'loc' or key == 'eeg_loc': + assert_array_equal(bp_info[key], 0) + elif key == 'coil_type': + assert_equal(bp_info[key], FIFF.FIFFV_COIL_EEG_BIPOLAR) + elif key == 'kind': + assert_equal(bp_info[key], FIFF.FIFFV_EOG_CH) + else: + assert_equal(bp_info[key], an_info[key]) + assert_equal(bp_info['extra'], 'some extra value') + + # Test creating a bipolar reference that doesn't involve EEG channels: + # it should not set the custom_ref_applied flag + reref = set_bipolar_reference(raw, 'MEG 0111', 'MEG 0112', + ch_info={'kind': FIFF.FIFFV_MEG_CH}) + assert_true(not reref.info['custom_ref_applied']) + assert_true('MEG 0111-MEG 0112' in reref.ch_names) + + # Test a battery of invalid inputs + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', ['EEG 002', 'EEG 003'], 'bipolar') + assert_raises(ValueError, set_bipolar_reference, raw, + ['EEG 001', 'EEG 002'], 'EEG 003', 'bipolar') + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', ['bipolar1', 'bipolar2']) + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', 'bipolar', + ch_info=[{'foo': 'bar'}, {'foo': 'bar'}]) + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', ch_name='EEG 003') diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index e3e3b863a67..609811e1251 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -15,6 +15,7 @@ from ..io.matrix import (_read_named_matrix, _transpose_named_matrix, write_named_matrix) from ..io.proj import _read_proj, make_projector, _write_proj +from ..io.proj import _has_eeg_average_ref_proj from ..io.tree import dir_tree_find from ..io.write import (write_int, write_float_matrix, start_file, start_block, end_block, end_file, write_float, @@ -716,6 +717,15 @@ def _check_ori(pick_ori, pick_normal): return pick_ori +def _check_reference(inst): + if "eeg" in inst and not _has_eeg_average_ref_proj(inst.info['projs']): + raise ValueError('EEG average reference is mandatory for inverse ' + 'modeling.') + if inst.info['custom_ref_applied']: + raise ValueError('Custom EEG reference is not allowed for inverse ' + 'modeling.') + + def _subject_from_inverse(inverse_operator): """Get subject id from inverse operator""" return inverse_operator['src'][0].get('subject_his_id', None) @@ -752,6 +762,7 @@ def apply_inverse(evoked, inverse_operator, lambda2, method="dSPM", stc : SourceEstimate | VolSourceEstimate The source estimates """ + _check_reference(evoked) method = _check_method(method) pick_ori = _check_ori(pick_ori, pick_normal) # @@ -850,6 +861,7 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, method="dSPM", stc : SourceEstimate | VolSourceEstimate The source estimates. """ + _check_reference(raw) method = _check_method(method) pick_ori = _check_ori(pick_ori, pick_normal) @@ -1017,6 +1029,7 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, method="dSPM", stc : list of SourceEstimate or VolSourceEstimate The source estimates for all epochs. """ + _check_reference(epochs) stcs = _apply_inverse_epochs_gen(epochs, inverse_operator, lambda2, method=method, label=label, nave=nave, pick_ori=pick_ori, verbose=verbose, diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index 5cc3e9953d8..755963c7a75 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -232,6 +232,13 @@ def test_apply_inverse_operator(): assert_true(stc.data.max() < 35) assert_true(stc.data.mean() > 0.1) + # Test we get errors when using custom ref or no average proj is present + evoked.info['custom_ref_applied'] = True + assert_raises(ValueError, apply_inverse, evoked, inv_op, lambda2, "MNE") + evoked.info['custom_ref_applied'] = False + evoked.info['projs'] = [] # remove EEG proj + assert_raises(ValueError, apply_inverse, evoked, inv_op, lambda2, "MNE") + @testing.requires_testing_data def test_make_inverse_operator_fixed(): diff --git a/mne/tests/test_epochs.py b/mne/tests/test_epochs.py index 8115c65bdc8..c85ecd6aebc 100644 --- a/mne/tests/test_epochs.py +++ b/mne/tests/test_epochs.py @@ -276,6 +276,12 @@ def test_epochs_proj(): baseline=(None, 0), proj=True, add_eeg_ref=False) assert_true(not _has_eeg_average_ref_proj(epochs.info['projs'])) + # make sure we don't add avg ref when a custom ref has been applied + raw.info['custom_ref_applied'] = True + epochs = Epochs(raw, events[:4], event_id, tmin, tmax, picks=this_picks, + baseline=(None, 0), proj=True) + assert_true(not _has_eeg_average_ref_proj(epochs.info['projs'])) + def test_evoked_arithmetic(): """Test arithmetic of evoked data diff --git a/mne/tests/test_proj.py b/mne/tests/test_proj.py index b8148624d66..1363f180953 100644 --- a/mne/tests/test_proj.py +++ b/mne/tests/test_proj.py @@ -1,5 +1,5 @@ import os.path as op -from nose.tools import assert_true +from nose.tools import assert_true, assert_raises import warnings import numpy as np @@ -13,8 +13,10 @@ from mne import pick_types from mne.io import Raw from mne import compute_proj_epochs, compute_proj_evoked, compute_proj_raw -from mne.io.proj import make_projector, activate_proj -from mne.proj import read_proj, write_proj, make_eeg_average_ref_proj +from mne.io.proj import (make_projector, activate_proj, + _needs_eeg_average_ref_proj) +from mne.proj import (read_proj, write_proj, make_eeg_average_ref_proj, + _has_eeg_average_ref_proj) from mne import read_events, Epochs, sensitivity_map, read_source_estimate from mne.utils import _TempDir, run_tests_if_main, clean_warning_registry @@ -217,4 +219,51 @@ def test_compute_proj_raw(): assert_array_almost_equal(proj, np.eye(len(raw.ch_names))) +def test_make_eeg_average_ref_proj(): + """Test EEG average reference projection""" + raw = Raw(raw_fname, add_eeg_ref=False, preload=True) + eeg = mne.pick_types(raw.info, meg=False, eeg=True) + + # No average EEG reference + assert_true(not np.all(raw._data[eeg].mean(axis=0) < 1e-19)) + + # Apply average EEG reference + car = make_eeg_average_ref_proj(raw.info) + reref = raw.copy() + reref.add_proj(car) + reref.apply_proj() + assert_array_almost_equal(reref._data[eeg].mean(axis=0), 0, decimal=19) + + # Error when custom reference has already been applied + raw.info['custom_ref_applied'] = True + assert_raises(RuntimeError, make_eeg_average_ref_proj, raw.info) + + +def test_has_eeg_average_ref_proj(): + """Test checking whether an EEG average reference exists""" + assert_true(not _has_eeg_average_ref_proj([])) + + raw = Raw(raw_fname, add_eeg_ref=True, preload=False) + assert_true(_has_eeg_average_ref_proj(raw.info['projs'])) + + +def test_needs_eeg_average_ref_proj(): + """Test checking whether a recording needs an EEG average reference""" + raw = Raw(raw_fname, add_eeg_ref=False, preload=False) + assert_true(_needs_eeg_average_ref_proj(raw.info)) + + raw = Raw(raw_fname, add_eeg_ref=True, preload=False) + assert_true(not _needs_eeg_average_ref_proj(raw.info)) + + # No EEG channels + raw = Raw(raw_fname, add_eeg_ref=False, preload=True) + eeg = [raw.ch_names[c] for c in pick_types(raw.info, meg=False, eeg=True)] + raw.drop_channels(eeg) + assert_true(not _needs_eeg_average_ref_proj(raw.info)) + + # Custom ref flag set + raw = Raw(raw_fname, add_eeg_ref=False, preload=False) + raw.info['custom_ref_applied'] = True + assert_true(not _needs_eeg_average_ref_proj(raw.info)) + run_tests_if_main() diff --git a/mne/viz/raw.py b/mne/viz/raw.py index 5b1f9b5210a..b2eb31944ac 100644 --- a/mne/viz/raw.py +++ b/mne/viz/raw.py @@ -401,7 +401,7 @@ def plot_raw(raw, events=None, duration=10.0, start=0.0, n_channels=None, else: if lowpass <= highpass: raise ValueError('lowpass (%s) must be > highpass (%s)' - %(lowpass, highpass)) + % (lowpass, highpass)) ba = butter(filtorder, [highpass / nyq, lowpass / nyq], 'bandpass', analog=False)