add spectrum class

[drammock]

Progress so far:

• [x] .compute_psd() method for Raw, Epochs, Evoked yields Spectrum object
• [x] supports method='welch' and method='multitaper'
• [x] supports unaggregated welch via average=False
• [x] supports unaggregated multitaper via output='complex'
• [x] properties ch_names, freqs, method (welch/multitaper)
• [x] private property ._dims, e.g., ('epoch', 'channel', 'freq', 'taper') for an unaggregated multitaper derived from Epochs
• [x] save to / load from HDF5 (spectrum.save(fname) / mne.time_frequency.read_spectrum(fname)) for Spectrum object
• [x] save/load for EpochsSpectrum object
• [x] __eq__ method for comparisons / testing
• [x] __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)
• [x] .pick() method
• [x] .units() method
• [x] .to_data_frame() method (works with raw/epochs/evoked, agged/unagged estimates)
• [x] .plot() method with equivalent API to old inst.plot_psd() method
• [ ] .plot_topomap() method (partially implemnted)
• [x] tests
• [ ] tutorial

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 :)

[larsoner]

@drammock the bullet point is still unchecked:

support passing a callable to aggregate welch/multitaper estimates

Should we wait for a follow-up PR, or did you want to do it here?

Either way I can look at the code early tomorrow

[drammock]

That was added to the list this week after conversations with @ajquinn. I think a follow up PR is probably best but if you disagree I can add it here. I don't think it will be too hard.

-------- Original Message -------- On Aug 25, 2022, 00:13, Eric Larson wrote:

@.***(https://github.com/drammock) the bullet point is still unchecked:

support passing a callable to aggregate welch/multitaper estimates

Should we wait for a follow-up PR, or did you want to do it here?

Either way I can look at the code early tomorrow

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>

[drammock]

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

[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]

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