MRG: Revise filtering (mne-tools#3277)

* ENH: Add generic filter function

* DOC: Minor fix [ci skip]

* DOC: Add some tutorial materials [circle full]

* FIX: Fix tests

* STY: Fix plots [ci skip]

* FIX: Minor fixes

* FIX: Spelling

* FIX: Minor fixes

* FIX: Typos

* FIX: Fix name

Makefile

@@ -5,8 +5,7 @@
 PYTHON ?= python
 NOSETESTS ?= nosetests
 CTAGS ?= ctags
-# The *.fif had to be there twice to be properly ignored (!)
-CODESPELL_SKIPS ?= "*.fif,*.fif,*.eve,*.gz,*.tgz,*.zip,*.mat,*.stc,*.label,*.w,*.bz2,*.annot,*.sulc,*.log,*.local-copy,*.orig_avg,*.inflated_avg,*.gii,*.pyc,*.doctree,*.pickle,*.inv,*.png,*.edf,*.touch,*.thickness,*.nofix,*.volume,*.defect_borders,*.mgh,lh.*,rh.*,COR-*,FreeSurferColorLUT.txt,*.examples,.xdebug_mris_calc,bad.segments,BadChannels,*.hist,empty_file,*.orig"
+CODESPELL_SKIPS ?= "*.fif,*.eve,*.gz,*.tgz,*.zip,*.mat,*.stc,*.label,*.w,*.bz2,*.annot,*.sulc,*.log,*.local-copy,*.orig_avg,*.inflated_avg,*.gii,*.pyc,*.doctree,*.pickle,*.inv,*.png,*.edf,*.touch,*.thickness,*.nofix,*.volume,*.defect_borders,*.mgh,lh.*,rh.*,COR-*,FreeSurferColorLUT.txt,*.examples,.xdebug_mris_calc,bad.segments,BadChannels,*.hist,empty_file,*.orig,*.js,*.map,*.ipynb"
 CODESPELL_DIRS ?= mne/ doc/ tutorials/ examples/
 all: clean inplace test test-doc
 

circle.yml

@@ -44,7 +44,7 @@ dependencies:
 
   override:
     - cd /home/ubuntu/mne-python && python setup.py develop;
-    - if [ "$CIRCLE_BRANCH" == "master" ] || [ "$CIRCLE_BRANCH" == "maint/0.12" ]; then
+    - if [ "$CIRCLE_BRANCH" == "master" ] || [ "$CIRCLE_BRANCH" == "maint/0.12" ] || [[ `git log -1 --pretty=%B` == *"[circle full]"* ]]; then
         mkdir -p ~/mne_data;
         python -c "import mne; mne.datasets._download_all_example_data()";
       fi
@@ -57,12 +57,12 @@ dependencies:
 
 test:
   override:
-    - if [ "$CIRCLE_BRANCH" == "master" ] || [ "$CIRCLE_BRANCH" == "maint/0.12" ]; then
+    - if [ "$CIRCLE_BRANCH" == "master" ] || [ "$CIRCLE_BRANCH" == "maint/0.12" ] || [[ `git log -1 --pretty=%B` == *"[circle full]"* ]]; then
         make test-doc;
       else
         cd doc && make html_dev-noplot;
       fi
-    - if [ "$CIRCLE_BRANCH" == "master" ]; then cd doc && make html_dev; fi:
+    - if [ "$CIRCLE_BRANCH" == "master" ] || [[ `git log -1 --pretty=%B` == *"[circle full]"* ]]; then cd doc && make html_dev; fi:
           timeout: 1500
     - if [ "$CIRCLE_BRANCH" == "maint/0.12" ]; then cd doc && make html_stable; fi:
           timeout: 1500

doc/conf.py

@@ -264,7 +264,11 @@
 trim_doctests_flags = True
 
 # Example configuration for intersphinx: refer to the Python standard library.
-intersphinx_mapping = {'http://docs.python.org/': None}
+intersphinx_mapping = {
+    'python': ('http://docs.python.org/', None),
+    'numpy': ('http://docs.scipy.org/doc/numpy-dev/', None),
+    'scipy': ('http://scipy.github.io/devdocs/', None),
+}
 
 examples_dirs = ['../examples', '../tutorials']
 gallery_dirs = ['auto_examples', 'auto_tutorials']
@@ -290,4 +294,5 @@
     'gallery_dirs': gallery_dirs,
     'find_mayavi_figures': find_mayavi_figures,
     'default_thumb_file': os.path.join('_static', 'mne_helmet.png'),
+    'mod_example_dir': 'generated',
     }

doc/faq.rst

@@ -84,8 +84,7 @@ If you want to write your own data to disk (e.g., subject behavioral
 scores), we strongly recommend using `h5io <https://github.com/h5io/h5io>`_,
 which is based on the
 `HDF5 format <https://en.wikipedia.org/wiki/Hierarchical_Data_Format>`_ and
-`h5py <http://www.h5py.org/>`_,
-to save data in a fast, future-compatible, standard format.
+h5py_, to save data in a fast, future-compatible, standard format.
 
 
 Resampling and decimating data

doc/manual/cookbook.rst

@@ -384,8 +384,6 @@ ways:
   stationary with respect to background brain activity. This can also
   use :func:`mne.compute_raw_covariance`.
 
-See :ref:`covariance` for more information.
-
 .. _CIHCFJEI:
 
 Calculating the inverse operator

doc/manual/io.rst

@@ -373,7 +373,7 @@ Arbitrary (e.g., simulated or manually read in) raw data can be constructed
 from memory by making use of :class:`mne.io.RawArray`, :class:`mne.EpochsArray`
 or :class:`mne.EvokedArray` in combination with :func:`mne.create_info`.
 
-This functionality is illustrated in :ref:`example_io_plot_objects_from_arrays.py` .
+This functionality is illustrated in :ref:`sphx_glr_auto_examples_io_plot_objects_from_arrays.py`.
 Using 3rd party libraries such as NEO (https://pythonhosted.org/neo/) in combination
 with these functions abundant electrophysiological file formats can be easily loaded
 into MNE.

doc/manual/source_localization/inverse.rst

@@ -184,7 +184,7 @@ Using the UNIX tools :ref:`mne_inverse_operator`, the values
 :math:`\varepsilon_k` can be adjusted with the regularization options
 ``--magreg`` , ``--gradreg`` , and ``--eegreg`` specified at the time of the
 inverse operator decomposition, see :ref:`inverse_operator`. The convenience script
-:ref:`mne_do_inverse_solution` has the ``--magreg`` and ``--gradreg`` combined to
+:ref:`mne_do_inverse_operator` has the ``--magreg`` and ``--gradreg`` combined to
 a single option, ``--megreg`` , see :ref:`CIHCFJEI`.
 Suggested range of values for :math:`\varepsilon_k` is :math:`0.05 \dotso 0.2`.
 

doc/python_reference.rst

@@ -819,7 +819,6 @@ Functions that operate on ``np.ndarray`` objects:
    cwt_morlet
    dpss_windows
    morlet
-   multitaper_psd
    single_trial_power
    stft
    istft

doc/tutorials.rst

@@ -50,6 +50,17 @@ For further reading:
     auto_tutorials/plot_python_intro.rst
     tutorials/seven_stories_about_mne.rst
 
+.. container:: span box
+
+  .. raw:: html
+
+    <h2>Background information</h2>
+
+  .. toctree::
+    :maxdepth: 1
+
+    auto_tutorials/plot_background_filtering.rst
+
 .. container:: span box
 
   .. raw:: html

doc/whats_new.rst

@@ -17,6 +17,8 @@ Changelog
 
     - Added the ability to decimate :class:`mne.Evoked` objects with :func:`mne.Evoked.decimate` by `Eric Larson`_
 
+    - Add generic array-filtering function :func:`mne.filter.filter_` by `Eric Larson`_
+
 BUG
 ~~~
 

mne/epochs.py

@@ -1538,29 +1538,30 @@ def __getitem__(self, item):
         -----
         Epochs can be accessed as ``epochs[...]`` in several ways:
 
-        1. ``epochs[idx]``: Return ``Epochs`` object with a subset of
-           epochs (supports single index and python-style slicing).
-
-        2. ``epochs['name']``: Return ``Epochs`` object with a copy of the
-           subset of epochs corresponding to an experimental condition as
-           specified by 'name'.
-
-            If conditions are tagged by names separated by '/' (e.g.
-            'audio/left', 'audio/right'), and 'name' is not in itself an
-            event key, this selects every event whose condition contains
-            the 'name' tag (e.g., 'left' matches 'audio/left' and
-            'visual/left'; but not 'audio_left'). Note that tags like
-            'auditory/left' and 'left/auditory' will be treated the
-            same way when accessed using tags.
-
-        3. ``epochs[['name_1', 'name_2', ... ]]``: Return ``Epochs`` object
-           with a copy of the subset of epochs corresponding to multiple
-            experimental conditions as specified by
-            ``'name_1', 'name_2', ...`` .
-
-            If conditions are separated by '/', selects every item containing
-            every list tag (e.g. ['audio', 'left'] selects 'audio/left' and
-            'audio/center/left', but not 'audio/right').
+            1. ``epochs[idx]``: Return ``Epochs`` object with a subset of
+               epochs (supports single index and python-style slicing).
+
+            2. ``epochs['name']``: Return ``Epochs`` object with a copy of the
+               subset of epochs corresponding to an experimental condition as
+               specified by 'name'.
+
+               If conditions are tagged by names separated by '/' (e.g.
+               'audio/left', 'audio/right'), and 'name' is not in itself an
+               event key, this selects every event whose condition contains
+               the 'name' tag (e.g., 'left' matches 'audio/left' and
+               'visual/left'; but not 'audio_left'). Note that tags like
+               'auditory/left' and 'left/auditory' will be treated the
+               same way when accessed using tags.
+
+            3. ``epochs[['name_1', 'name_2', ... ]]``: Return ``Epochs`` object
+               with a copy of the subset of epochs corresponding to multiple
+               experimental conditions as specified by
+               ``'name_1', 'name_2', ...`` .
+
+               If conditions are separated by '/', selects every item
+               containing every list tag (e.g. ['audio', 'left'] selects
+               'audio/left' and 'audio/center/left', but not 'audio/right').
+
         """
         data = self._data
         del self._data

mne/filter.py

@@ -278,11 +278,10 @@ def _filter(x, Fs, freq, gain, filter_length='10s', picks=None, n_jobs=1,
     gain : 1d array
         Filter gain at frequency sampling points.
     filter_length : str (Default: '10s') | int | None
-        Length of the filter to use. If None or "len(x) < filter_length",
-        the filter length used is len(x). Otherwise, if int, overlap-add
-        filtering with a filter of the specified length in samples) is
-        used (faster for long signals). If str, a human-readable time in
-        units of "s" or "ms" (e.g., "10s" or "5500ms") will be converted
+        Length of the filter to use. If None or ``len(x) < filter_length``,
+        the filter length used is ``len(x)``. If int, a filter of the
+        specified length in samples is used. If str, a human-readable time
+        in units of "s" or "ms" (e.g., "10s" or "5500ms") will be converted
         to the shortest power-of-two length at least that duration.
     picks : array-like of int | None
         Indices to filter. If None all indices will be filtered.
@@ -409,6 +408,7 @@ def _filtfilt(x, b, a, padlen, picks, n_jobs, copy):
 
 def _estimate_ringing_samples(b, a):
     """Helper function for determining IIR padding"""
+    # XXX Need to extend this to more than 1000 samples for long IIR filters!
     from scipy.signal import lfilter
     x = np.zeros(1000)
     x[0] = 1
@@ -572,6 +572,127 @@ def _check_method(method, iir_params, extra_types):
     return iir_params
 
 
+@verbose
+def filter_data(data, sfreq, l_freq, h_freq, picks=None, filter_length='10s',
+                l_trans_bandwidth=0.5, h_trans_bandwidth=0.5, n_jobs=1,
+                method='fft', iir_params=None, copy=True, verbose=None):
+    """Filter a subset of channels.
+
+    Applies a zero-phase low-pass, high-pass, band-pass, or band-stop
+    filter to the channels selected by ``picks``.
+
+    ``l_freq`` and ``h_freq`` are the frequencies below which and above
+    which, respectively, to filter out of the data. Thus the uses are:
+
+        * ``l_freq < h_freq``: band-pass filter
+        * ``l_freq > h_freq``: band-stop filter
+        * ``l_freq is not None and h_freq is None``: high-pass filter
+        * ``l_freq is None and h_freq is not None``: low-pass filter
+
+    .. note:: If n_jobs > 1, more memory is required as
+              ``len(picks) * n_times`` additional time points need to
+              be temporaily stored in memory.
+
+    Parameters
+    ----------
+    data : ndarray, shape (n_channels, n_times)
+        The data to filter.
+    sfreq : float
+        The sample frequency in Hz.
+    l_freq : float | None
+        Low cut-off frequency in Hz. If None the data are only low-passed.
+    h_freq : float | None
+        High cut-off frequency in Hz. If None the data are only
+        high-passed.
+    picks : array-like of int | None
+        Indices of channels to filter. If None all channels will be
+        filtered.
+    filter_length : str (Default: '10s') | int | None
+        Length of the filter to use. If None or ``len(x) < filter_length``,
+        the filter length used is ``len(x)``. If int, a filter of the
+        specified length in samples is used. If str, a human-readable time
+        in units of "s" or "ms" (e.g., "10s" or "5500ms") will be converted
+        to the shortest power-of-two length at least that duration.
+        Not used for 'iir' filters.
+    l_trans_bandwidth : float
+        Width of the transition band at the low cut-off frequency in Hz
+        (high pass or cutoff 1 in bandpass). Not used for 'iir' filters.
+    h_trans_bandwidth : float
+        Width of the transition band at the high cut-off frequency in Hz
+        (low pass or cutoff 2 in bandpass). Not used for 'iir' filters.
+    n_jobs : int | str
+        Number of jobs to run in parallel. Can be 'cuda' if scikits.cuda
+        is installed properly, CUDA is initialized, and method='fft'.
+    method : str
+        'fft' will use overlap-add FIR filtering, 'iir' will use IIR
+        forward-backward filtering (via filtfilt).
+    iir_params : dict | None
+        Dictionary of parameters to use for IIR filtering.
+        See mne.filter.construct_iir_filter for details. If iir_params
+        is None and method="iir", 4th order Butterworth will be used.
+    copy : bool
+        If True, a copy of x, filtered, is returned. Otherwise, it operates
+        on x in place.
+    verbose : bool, str, int, or None
+        If not None, override default verbose level (see mne.verbose).
+        Defaults to self.verbose.
+
+    Returns
+    -------
+    data : ndarray, shape (n_channels, n_times)
+        The filtered data.
+    """
+    if not isinstance(data, np.ndarray) or data.ndim != 2:
+        raise ValueError('data must be an array with two dimensions')
+    sfreq = float(sfreq)
+    if sfreq < 0:
+        raise ValueError('sfreq must be positive')
+    if h_freq is not None:
+        h_freq = float(h_freq)
+        if h_freq > (sfreq / 2.):
+            raise ValueError('h_freq (%s) must be less than the Nyquist '
+                             'frequency %s' % (h_freq, sfreq / 2.))
+    if l_freq is not None:
+        l_freq = float(l_freq)
+        if l_freq == 0:
+            l_freq = None
+    if l_freq is None and h_freq is not None:
+        logger.info('Low-pass filtering at %0.2g Hz' % h_freq)
+        low_pass_filter(data, sfreq, h_freq,
+                        filter_length=filter_length,
+                        trans_bandwidth=h_trans_bandwidth, method=method,
+                        iir_params=iir_params, picks=picks, n_jobs=n_jobs,
+                        copy=copy)
+    if l_freq is not None and h_freq is None:
+        logger.info('High-pass filtering at %0.2g Hz' % l_freq)
+        high_pass_filter(
+            data, sfreq, l_freq, filter_length=filter_length,
+            trans_bandwidth=l_trans_bandwidth, method=method,
+            iir_params=iir_params, picks=picks, n_jobs=n_jobs, copy=copy)
+    if l_freq is not None and h_freq is not None:
+        if l_freq < h_freq:
+            logger.info('Band-pass filtering from %0.2g - %0.2g Hz'
+                        % (l_freq, h_freq))
+            data = band_pass_filter(
+                data, sfreq, l_freq, h_freq,
+                filter_length=filter_length,
+                l_trans_bandwidth=l_trans_bandwidth,
+                h_trans_bandwidth=h_trans_bandwidth,
+                method=method, iir_params=iir_params, picks=picks,
+                n_jobs=n_jobs, copy=copy)
+        else:
+            logger.info('Band-stop filtering from %0.2g - %0.2g Hz'
+                        % (h_freq, l_freq))
+            data = band_stop_filter(
+                data, sfreq, h_freq, l_freq,
+                filter_length=filter_length,
+                l_trans_bandwidth=h_trans_bandwidth,
+                h_trans_bandwidth=l_trans_bandwidth, method=method,
+                iir_params=iir_params, picks=picks, n_jobs=n_jobs,
+                copy=copy)
+    return data
+
+
 @verbose
 def band_pass_filter(x, Fs, Fp1, Fp2, filter_length='10s',
                      l_trans_bandwidth=0.5, h_trans_bandwidth=0.5,