add spectrum class

[drammock]

Progress so far:

• .compute_psd() method for Raw, Epochs, Evoked yields Spectrum object
• supports method='welch' and method='multitaper'
• supports unaggregated welch via average=False
• supports unaggregated multitaper via output='complex'
• support passing a callable to aggregate welch/multitaper estimates support passing a callable to aggregate welch/multitaper estimates in Spectrum constructor #11093
• properties ch_names, freqs, method (welch/multitaper)
• private property ._dims, e.g., ('epoch', 'channel', 'freq', 'taper') for an unaggregated multitaper derived from Epochs
• save to / load from HDF5 (spectrum.save(fname) / mne.time_frequency.read_spectrum(fname)) for Spectrum object
• save/load for EpochsSpectrum object
• __eq__ method for comparisons / testing
• __getitem__ works for raw/epochs/evoked (raw/evoked-derived spectra work like np.array; epochs-derived spectra work like epochs: select by condition name, pandas query, etc)
• .pick() method
• .units() method
• .to_data_frame() method (works with raw/epochs/evoked, agged/unagged estimates)
• .plot() method with equivalent API to old inst.plot_psd() method
• .plot_topomap() method
• .plot_topo() method
• tests
• tutorial
• deprecations
• changelog

import os
import mne

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                                    'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file, verbose=False, preload=False)

spectrum = raw.compute_psd()
spectrum.plot()

yields this:

spectrum has properties _data, freqs, ch_names, method (welch or multitaper), an info dict, and spec._dims (description of the data dims). For example

In [8]: unagg_spec = raw.to_spectrum(average=None)
In [9]: unagg_spec._dims
Out[9]: ('channel', 'freq', 'segment')

Here is the corresponding repr <Spectrum (from Raw) | 364 channels × 129 freqs × 651 segments, 0.0-300.3 Hz>

Eventually (in another PR) I'd like to refactor the multitaper functions to make it possible to return the unaggregated tapers EDIT I think this was done already while I was on leave, in which case you'd expect something like raw.to_spectrum(multitaper, aggregate=False)._dims to be ('channel', 'freq', 'taper').

epochs.compute_psd() and evoked.compute_psd() also work:

events = mne.find_events(raw, stim_channel='STI 014')
event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3}
epochs = mne.Epochs(raw, events, event_id=event_dict, preload=False,
                    proj=False)
evk = epochs.average()

epo_spec = epochs.compute_psd()
evk_spec = evk.compute_psd()

[drammock]

this is still a draft but now would be a good time for some early feedback about code organization and API (ping @larsoner and @agramfort in particular, but input is welcome from anyone).

[agramfort]

I am wondering about the use of spectrum. We use psd everywhere else.

my reflex would have been

psd = raw.compute_psd()
psd.plot()

would be the same as what we have now

raw.plot_psd()

[drammock]

I am wondering about the use of spectrum. We use psd everywhere else.

I chose "spectrum" because it's more general: e.g., depending on what parameters people use, sometimes we plot the amplitude spectrum, or we plot 10*log10 of the power spectrum.

[drammock]

Other reasons to pick a new name:

• people are used to things like psds, freqs = psd_welch(raw) so maybe (?) raw.compute_psd() might be similar enough to generate wrong expectations about the return value?
• having a more distinct name makes it less surprising that there are API changes. For example the previous psd plotting functions have area_mode parameter which I'm proposing to change to ci in spectrum.plot() (ci is used in seaborn, and also in our own plot_compare_evokeds).

[agramfort]

ok makes sense

…

Message ID: ***@***.***>

[dengemann]

@dengemann tagging myself to follow. Great this is alive!

[dengemann]

Really happy to see this @drammock. Once basics are settled it would be cool to also add plot_topomap, based on Epochs.plot_psd_topomap. Not to forget adjusted .to_data_frame. Not suggesting all this should happen in this PR though :)

[drammock]

@larsoner this test fails for evoked because of a limitation in pytest. The failure is:

_____________________________________________ test_spectrum_io[evoked] ______________________________________________
The requested fixture has no parameter defined for test:
    mne/time_frequency/tests/test_spectrum.py::test_spectrum_io[evoked]

Requested fixture '_evoked' defined in:
mne/conftest.py:306

Requested here:
/opt/miniforge3/envs/mnedev/lib/python3.10/site-packages/_pytest/fixtures.py:659

I can make it go away by removing params=[testing._pytest_param()] from the signature here:

mne-python/mne/conftest.py

Lines 305 to 312 in 2ba483a

	@pytest.fixture(scope='session', params=[testing._pytest_param()])
	def _evoked():
	# This one is session scoped, so be sure not to modify it (use evoked
	# instead)
	evoked = mne.read_evokeds(fname_evoked, condition='Left Auditory',
	baseline=(None, 0))
	evoked.crop(0, 0.2)
	return evoked

But presumably that will break your magic auto-skipping behavior. The limitation is that request.getfixturevalue() cannot handle fixtures that have params, and this is apparently by design. A workaround is to use the pytest-lazy-fixture package. Another alternative is to copy-paste a separate test for Evoked IO, like this:

@requires_h5py
def test_evoked_spectrum_io(evoked, tmp_path):
    """Test save/load of spectrum objects."""
    fname = tmp_path / 'spectrum.h5'
    orig = evoked.compute_psd()
    orig.save(fname)
    loaded = read_spectrum(fname)
    assert orig == loaded

... which works because now the evoked fixture is listed directly in the test signature (instead of being accessed via getfixturevalue()), but is not very DRY.

What do you recommend?

[larsoner]

I would work around it with:

@requires_h5py
@pytest.mark.parametrize('inst', ('raw', 'epochs', 'evoked'))
def test_spectrum_io(inst, tmp_path, request, evoked):
    """Test save/load of spectrum objects."""
    # work around parametrized fixtures not being accessible
    # via request.getfixturevalue
    # https://github.com/pytest-dev/pytest/issues/4666#issuecomment-456593913  # noqa: E501
    if inst == 'evoked':
        inst = evoked
    else:
        inst = request.getfixturevalue(inst)
    fname = tmp_path / 'spectrum.h5'
    orig = inst.compute_psd()
    orig.save(fname)
    loaded = read_spectrum(fname)
    assert orig == loaded

[drammock]

Ah, OK. Inelegant but practical. Will do.

[drammock]

• new tutorial
• changed tutorial sleep stage classification
• changed tutorial raw plotting methods
• changed tutorial visualizing epochs

[wmvanvliet]

It's a work of art! 🤩

Do we at some point want to support users creating their own Spectral instances like EpochsArray and EvokedArray?

[drammock]

Do we at some point want to support users creating their own Spectral instances like EpochsArray and EvokedArray?

this is something I struggled with, and ultimately decided "YAGNI". But, I suspect it wouldn't be too difficult given that it's already basically implemented in the _from_file method for loading from disk.

[wmvanvliet]

ultimately decided "YAGNI"

Fair enough.

[larsoner]

LGTM +1 for merge after the _ci fix!

[larsoner]

https://output.circle-artifacts.com/output/job/79b51cb7-962f-4095-87ce-5aa72b49dcc1/artifacts/0/dev/auto_tutorials/time-freq/10_spectrum_class.html#sphx-glr-auto-tutorials-time-freq-10-spectrum-class-py

[larsoner]

Congrats @drammock !

[wmvanvliet]

Nice one!

[agramfort]

🎉 🍻 👏 @drammock