MRG : EEG Reference reloaded !

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

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
 ~~~
 

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 <w.m.vanvliet@gmail.com>
+#
+# 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()

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, '

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, '

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
 
 

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)

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
 

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

mne/io/base.py

@@ -3,6 +3,7 @@
 #          Martin Luessi <mluessi@nmr.mgh.harvard.edu>
 #          Denis Engemann <denis.engemann@gmail.com>
 #          Teon Brooks <teon.brooks@gmail.com>
+#          Marijn van Vliet <w.m.vanvliet@gmail.com>
 #
 # 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)

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 = {}