Skip to content

ENH: Add polyphase resampling #12268

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Dec 6, 2023
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion doc/Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -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
1 change: 1 addition & 0 deletions doc/changes/devel.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
~~~~
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -40,114 +27,77 @@
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')
4 changes: 2 additions & 2 deletions examples/decoding/receptive_field_mtrf.py
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sphinx_gallery_thumbnail_number=3 is visible in the rendered tutorial, can you fix that here please?

Original file line number Diff line number Diff line change
Expand Up @@ -58,8 +58,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")
Expand Down
8 changes: 4 additions & 4 deletions examples/time_frequency/source_power_spectrum_opm.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"]

Expand Down
2 changes: 1 addition & 1 deletion mne/cuda.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading