Skip to content
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

[MRG] ENH/DOC: demo custom spectrum creation #11493

Merged
merged 7 commits into from
Feb 19, 2023
Merged
Changes from 3 commits
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
45 changes: 45 additions & 0 deletions tutorials/simulation/10_array_objs.py
Original file line number Diff line number Diff line change
Expand Up @@ -195,3 +195,48 @@
nave=data.shape[0], comment='simulated')
print(evoked_array)
evoked_array.plot()


# %%
# In certain situations you may wish to use a custom time-frequency
# decomposition for estimation of power spectra. Or you may wish to
# to process pre-computed power spectra in MNE.
dengemann marked this conversation as resolved.
Show resolved Hide resolved
# Following the same logic, it is possible to instantiate averaged power
# spectrum using the `~mne.time_frequency.Spectrum` class. This is slightly
dengemann marked this conversation as resolved.
Show resolved Hide resolved
# experimental at the moment but works. An API for doing this may follow.

# compute power spectrum

psd, freqs = mne.time_frequency.psd_array_welch(
data, info['sfreq'], n_fft=128, n_per_seg=32)

psd_ave = psd.mean(0)

# map to `~mne.time_frequency.Spectrum` class and explore API


def spectrum_from_array(data: np.ndarray, # spectral features
freqs: np.ndarray, # frequencies
inst_info: mne.Info # the meta data of MNE instance
) -> mne.time_frequency.Spectrum: # Spectrum object
dengemann marked this conversation as resolved.
Show resolved Hide resolved
"""Create MNE averaged power spectrum object from custom data"""
state = dict(
method='my_welch',
data=data,
sfreq=inst_info['sfreq'],
dims=('channel', 'freq'),
freqs=freqs,
inst_type_str='Raw',
data_type='simulated signals',
Copy link
Member Author

Choose a reason for hiding this comment

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

@drammock I put this for illustration purposes. Are there any expectations in the Spectrum code regarding this field? It was 'Averaged EEG' in my concrete case on which this example is based.

Copy link
Member

Choose a reason for hiding this comment

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

right now it is only used in the repr of the object. I don't have any concrete plans to use it in any other way (and I would actively avoid using it in if clauses to triage different flavors of spectrum, it's probably not reliable to do so)

Copy link
Member

Choose a reason for hiding this comment

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

native values for this are Fourier Coefficients, Power Spectrum, and Averaged Power Spectrum. Averaged Fourier Coefficients is not supported, see

# averaging unaggregated spectral estimates are not supported
if hasattr(self, '_mt_weights'):
raise NotImplementedError(
'Averaging complex spectra is not supported. Consider '
'averaging the signals before computing the complex spectrum.')
elif 'segment' in self._dims:
raise NotImplementedError(
'Averaging individual Welch segments across epochs is not '
'supported. Consider averaging the signals before computing '
'the Welch spectrum estimates.')

info=inst_info,
)
defaults = dict(
method=None, fmin=None, fmax=None, tmin=None, tmax=None,
picks=None, proj=None, reject_by_annotation=None, n_jobs=None,
verbose=None
)
return mne.time_frequency.Spectrum(state, **defaults)


spectrum = spectrum_from_array(data=psd_ave, freqs=freqs, inst_info=info)
spectrum.plot(picks=[0, 1])