ENH: Add polyphase resampling (#12268)

Co-authored-by: Daniel McCloy <[email protected]>

Author: larsoner

doc/Makefile

@@ -76,6 +76,6 @@ doctest:
 	      "results in _build/doctest/output.txt."
 
 view:
-	@python -c "import webbrowser; webbrowser.open_new_tab('file://$(PWD)/_build/html/index.html')"
+	@python -c "import webbrowser; webbrowser.open_new_tab('file://$(PWD)/_build/html/sg_execution_times.html')"
 
 show: view

doc/changes/devel.rst

@@ -36,6 +36,7 @@ Enhancements
 ~~~~~~~~~~~~
 - Speed up export to .edf in :func:`mne.export.export_raw` by using ``edfio`` instead of ``EDFlib-Python`` (:gh:`12218` by :newcontrib:`Florian Hofer`)
 - We added type hints for the return values of :func:`mne.read_evokeds` and :func:`mne.io.read_raw`. Development environments like VS Code or PyCharm will now provide more help when using these functions in your code. (:gh:`12250` by `Richard Höchenberger`_ and `Eric Larson`_)
+- Add ``method="polyphase"`` to :meth:`mne.io.Raw.resample` and related functions to allow resampling using :func:`scipy.signal.upfirdn` (:gh:`12268` by `Eric Larson`_)
 
 Bugs
 ~~~~

examples/datasets/spm_faces_dataset_sgskip.py renamed to examples/datasets/spm_faces_dataset.py

@@ -5,28 +5,15 @@
 From raw data to dSPM on SPM Faces dataset
 ==========================================
 
-Runs a full pipeline using MNE-Python:
-
-    - artifact removal
-    - averaging Epochs
-    - forward model computation
-    - source reconstruction using dSPM on the contrast : "faces - scrambled"
-
-.. note:: This example does quite a bit of processing, so even on a
-          fast machine it can take several minutes to complete.
+Runs a full pipeline using MNE-Python. This example does quite a bit of processing, so
+even on a fast machine it can take several minutes to complete.
 """
 # Authors: Alexandre Gramfort <[email protected]>
 #          Denis Engemann <[email protected]>
 #
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
-# %%
-
-# sphinx_gallery_thumbnail_number = 10
-
-import matplotlib.pyplot as plt
-
 import mne
 from mne import combine_evoked, io
 from mne.datasets import spm_face
@@ -40,114 +27,76 @@
 spm_path = data_path / "MEG" / "spm"
 
 # %%
-# Load and filter data, set up epochs
+# Load data, filter it, and fit ICA.
 
 raw_fname = spm_path / "SPM_CTF_MEG_example_faces1_3D.ds"
-
 raw = io.read_raw_ctf(raw_fname, preload=True)  # Take first run
 # Here to save memory and time we'll downsample heavily -- this is not
 # advised for real data as it can effectively jitter events!
-raw.resample(120.0, npad="auto")
-
-picks = mne.pick_types(raw.info, meg=True, exclude="bads")
-raw.filter(1, 30, method="fir", fir_design="firwin")
+raw.resample(100)
+raw.filter(1.0, None)  # high-pass
+reject = dict(mag=5e-12)
+ica = ICA(n_components=0.95, max_iter="auto", random_state=0)
+ica.fit(raw, reject=reject)
+# compute correlation scores, get bad indices sorted by score
+eog_epochs = create_eog_epochs(raw, ch_name="MRT31-2908", reject=reject)
+eog_inds, eog_scores = ica.find_bads_eog(eog_epochs, ch_name="MRT31-2908")
+ica.plot_scores(eog_scores, eog_inds)  # see scores the selection is based on
+ica.plot_components(eog_inds)  # view topographic sensitivity of components
+ica.exclude += eog_inds[:1]  # we saw the 2nd ECG component looked too dipolar
+ica.plot_overlay(eog_epochs.average())  # inspect artifact removal
 
+# %%
+# Epoch data and apply ICA.
 events = mne.find_events(raw, stim_channel="UPPT001")
-
-# plot the events to get an idea of the paradigm
-mne.viz.plot_events(events, raw.info["sfreq"])
-
 event_ids = {"faces": 1, "scrambled": 2}
-
 tmin, tmax = -0.2, 0.6
-baseline = None  # no baseline as high-pass is applied
-reject = dict(mag=5e-12)
-
 epochs = mne.Epochs(
     raw,
     events,
     event_ids,
     tmin,
     tmax,
-    picks=picks,
-    baseline=baseline,
+    picks="meg",
+    baseline=None,
     preload=True,
     reject=reject,
 )
-
-# Fit ICA, find and remove major artifacts
-ica = ICA(n_components=0.95, max_iter="auto", random_state=0)
-ica.fit(raw, decim=1, reject=reject)
-
-# compute correlation scores, get bad indices sorted by score
-eog_epochs = create_eog_epochs(raw, ch_name="MRT31-2908", reject=reject)
-eog_inds, eog_scores = ica.find_bads_eog(eog_epochs, ch_name="MRT31-2908")
-ica.plot_scores(eog_scores, eog_inds)  # see scores the selection is based on
-ica.plot_components(eog_inds)  # view topographic sensitivity of components
-ica.exclude += eog_inds[:1]  # we saw the 2nd ECG component looked too dipolar
-ica.plot_overlay(eog_epochs.average())  # inspect artifact removal
+del raw
 ica.apply(epochs)  # clean data, default in place
-
 evoked = [epochs[k].average() for k in event_ids]
-
 contrast = combine_evoked(evoked, weights=[-1, 1])  # Faces - scrambled
-
 evoked.append(contrast)
-
 for e in evoked:
     e.plot(ylim=dict(mag=[-400, 400]))
 
-plt.show()
-
-# estimate noise covarariance
-noise_cov = mne.compute_covariance(epochs, tmax=0, method="shrunk", rank=None)
-
 # %%
-# Visualize fields on MEG helmet
-
-# The transformation here was aligned using the dig-montage. It's included in
-# the spm_faces dataset and is named SPM_dig_montage.fif.
-trans_fname = spm_path / "SPM_CTF_MEG_example_faces1_3D_raw-trans.fif"
-
-maps = mne.make_field_map(
-    evoked[0], trans_fname, subject="spm", subjects_dir=subjects_dir, n_jobs=None
-)
-
-evoked[0].plot_field(maps, time=0.170, time_viewer=False)
-
-# %%
-# Look at the whitened evoked daat
+# Estimate noise covariance and look at the whitened evoked data
 
+noise_cov = mne.compute_covariance(epochs, tmax=0, method="shrunk", rank=None)
 evoked[0].plot_white(noise_cov)
 
 # %%
 # Compute forward model
 
+trans_fname = spm_path / "SPM_CTF_MEG_example_faces1_3D_raw-trans.fif"
 src = subjects_dir / "spm" / "bem" / "spm-oct-6-src.fif"
 bem = subjects_dir / "spm" / "bem" / "spm-5120-5120-5120-bem-sol.fif"
 forward = mne.make_forward_solution(contrast.info, trans_fname, src, bem)
 
 # %%
-# Compute inverse solution
+# Compute inverse solution and plot
+
+# sphinx_gallery_thumbnail_number = 8
 
 snr = 3.0
 lambda2 = 1.0 / snr**2
-method = "dSPM"
-
-inverse_operator = make_inverse_operator(
-    contrast.info, forward, noise_cov, loose=0.2, depth=0.8
-)
-
-# Compute inverse solution on contrast
-stc = apply_inverse(contrast, inverse_operator, lambda2, method, pick_ori=None)
-# stc.save('spm_%s_dSPM_inverse' % contrast.comment)
-
-# Plot contrast in 3D with mne.viz.Brain if available
+inverse_operator = make_inverse_operator(contrast.info, forward, noise_cov)
+stc = apply_inverse(contrast, inverse_operator, lambda2, method="dSPM", pick_ori=None)
 brain = stc.plot(
     hemi="both",
     subjects_dir=subjects_dir,
     initial_time=0.170,
     views=["ven"],
     clim={"kind": "value", "lims": [3.0, 6.0, 9.0]},
 )
-# brain.save_image('dSPM_map.png')

examples/decoding/receptive_field_mtrf.py

@@ -17,7 +17,7 @@
 .. _figure 1: https://www.frontiersin.org/articles/10.3389/fnhum.2016.00604/full#F1
 .. _figure 2: https://www.frontiersin.org/articles/10.3389/fnhum.2016.00604/full#F2
 .. _figure 5: https://www.frontiersin.org/articles/10.3389/fnhum.2016.00604/full#F5
-"""  # noqa: E501
+"""
 
 # Authors: Chris Holdgraf <[email protected]>
 #          Eric Larson <[email protected]>
@@ -26,9 +26,6 @@
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
-# %%
-# sphinx_gallery_thumbnail_number = 3
-
 from os.path import join
 
 import matplotlib.pyplot as plt
@@ -58,8 +55,8 @@
 speech = data["envelope"].T
 sfreq = float(data["Fs"].item())
 sfreq /= decim
-speech = mne.filter.resample(speech, down=decim, npad="auto")
-raw = mne.filter.resample(raw, down=decim, npad="auto")
+speech = mne.filter.resample(speech, down=decim, method="polyphase")
+raw = mne.filter.resample(raw, down=decim, method="polyphase")
 
 # Read in channel positions and create our MNE objects from the raw data
 montage = mne.channels.make_standard_montage("biosemi128")
@@ -131,6 +128,8 @@
 # across the scalp. We will recreate `figure 1`_ and `figure 2`_ from
 # :footcite:`CrosseEtAl2016`.
 
+# sphinx_gallery_thumbnail_number = 3
+
 # Print mean coefficients across all time delays / channels (see Fig 1)
 time_plot = 0.180  # For highlighting a specific time.
 fig, ax = plt.subplots(figsize=(4, 8), layout="constrained")

examples/time_frequency/source_power_spectrum_opm.py

@@ -58,16 +58,16 @@
 raw_erms = dict()
 new_sfreq = 60.0  # Nyquist frequency (30 Hz) < line noise freq (50 Hz)
 raws["vv"] = mne.io.read_raw_fif(vv_fname, verbose="error")  # ignore naming
-raws["vv"].load_data().resample(new_sfreq)
+raws["vv"].load_data().resample(new_sfreq, method="polyphase")
 raws["vv"].info["bads"] = ["MEG2233", "MEG1842"]
 raw_erms["vv"] = mne.io.read_raw_fif(vv_erm_fname, verbose="error")
-raw_erms["vv"].load_data().resample(new_sfreq)
+raw_erms["vv"].load_data().resample(new_sfreq, method="polyphase")
 raw_erms["vv"].info["bads"] = ["MEG2233", "MEG1842"]
 
 raws["opm"] = mne.io.read_raw_fif(opm_fname)
-raws["opm"].load_data().resample(new_sfreq)
+raws["opm"].load_data().resample(new_sfreq, method="polyphase")
 raw_erms["opm"] = mne.io.read_raw_fif(opm_erm_fname)
-raw_erms["opm"].load_data().resample(new_sfreq)
+raw_erms["opm"].load_data().resample(new_sfreq, method="polyphase")
 # Make sure our assumptions later hold
 assert raws["opm"].info["sfreq"] == raws["vv"].info["sfreq"]
 

mne/cuda.py

@@ -330,7 +330,7 @@ def _fft_resample(x, new_len, npads, to_removes, cuda_dict=None, pad="reflect_li
         Number of samples to remove after resampling.
     cuda_dict : dict
         Dictionary constructed using setup_cuda_multiply_repeated().
-    %(pad)s
+    %(pad_resample)s
         The default is ``'reflect_limited'``.
 
         .. versionadded:: 0.15