Mne-python: ENH: Add classes for Spectra / PSD

[thht]

Describe the problem

mne features dedicated classes for time-series data as well as time-frequency data. however, for a current project i need to work a lot with spectra. this would preferably include good exploratory plotting solutions like we already have for time-series data, including interactive features.

Describe your solution

as i am going to need these features anyhow, i am going to start working on such a class right now in my own collection.

i would like to know whether there is any motivation to have such a class in mne python.

and if yes, i would look forward to any suggestions / input from you guys.

[larsoner]

@drammock is also fighting with this issue. One solution is to consider these TFRs with a single time point, and make our computation functions like compute_source_psd return them.

[thht]

ok, but does this let me plot spectra and then mark a time-window which results in a topoplot etc?

[larsoner]

Not sure, need to check the docs for those classes. But if not, these seem like core features that should be part of the class anyway so we should add them

[thht]

ok... so i am going to start by creating a subclass of AverageTFR in my own repo, add the methods and as soon as it works i start an mne branch and refactor it into the original class... and then we can see how to make it mne-ready....

[larsoner]

Sounds reasonable to me but let's wait to hear from at least @drammock and @agramfort that they agree before you put in too much effort

[larsoner]

... actually I'm not sure why it would need to be a subclass rather than an instance of AverageTFR directly. If it's going to have to be a subclass we'd probably just want to make a new class altogether.

I know @drammock's use case involves storing multitaper estimates, so you still end up with 3D data of shape (n_channels, n_freqs, n_tapers), so we could put this in AverageTFR easily, but it isn't really a time frequency representation, it is just a single frequency representation computed with multiple window functions.

[thht]

sorry, this seems to be a misunderstanding.

in the final mne python version this would be all implemented in AverageTFR. but as i need this now and still want to use the stable mne version in my analysis project, i will create my custom subclass of AverageTFR just adding the features i need because i do not want to change the mne code in this particular scenario.

when i am ready to integrate it into mne, i will refactor the code so that my custom subclass is not needed anymore but instead all new logic/functionality would be implemented directly in AverageTFR.

i hope this clears up the confusion?

[larsoner]

Indeed, makes sense

[drammock]

Maybe I'm still not understanding your use case, @thht, but I don't see how a subclass approach will work in the long run. EpochsTFR, for example, does not even have a .plot() method. And the default .plot() for AverageTFR is an imageplot, not a line spectrum.

If you can make AverageTFR work for your use case, definitely go for it... I'm just skeptical about how it will integrate into the main codebase later, mostly because the kinds of plots that are useful / make sense are different between spectral and spectrotemporal data.

[thht]

good point... but the good thing here is, if i implement my personal version like this, we would still have both options for the final mne implementation. either we create a specialized class for spectrum data so we just need to copy the missing pieces from AverageTFR and then make it an independent class. or we copy the code i write to AverageTFR. just boils down to what code is going to be copied where.

alright, i'm going to do some coding now.

thanks for the input and let's see where this discussion is going!

[larsoner]

Maybe a Specta class that can be a single or multiple spectral estimates is the way to go?

• stores data as (n_channels, n_frequencies, n_estimates)
• has a plot method that produces amplitude + phase as a function of freq
• if you do my_avg_tfr.as_spectra(time) this returns an instance of Spectra

Something like this? There would be remaining open questions that we could then solve during development/actual use, like:

• in .plot, how do we combine the spectra (last dim)? If we start out by assuming that the length is either 1 (means a standard FFT) or multiple (means a multitaper estimate) then it's trivial to do so (hopefully)

Other questions would be, do we need an [epochs,] dimension? Hopefully not to start...

Would such a Spectra class work for your use cas @thht ?

[thht]

yes, that would be exactly what i am looking for.

[thht]

one issue might be that when doing multitaper, you might end up with a non-uniform amount of tapers for the frequencies. no idea how this is implemented in mne, but fieldtrip offers to calculate the amount of tapers individually per frequency...

[larsoner]

but fieldtrip offers to calculate the amount of tapers individually per frequency...

I'm guessing this is like the adaptive=True in MNE. If so, this is done after computing all tapers for all frequencies (window with Slepian -> FFT -> combine across tapers to get the PSD) so is really a question of how to combine the complex multitaper values rather than what complex values to compute in the first place. And I'm not too worried about storage unless we start talking about allowing (tons of) epochs and also being in source space, etc.

[dengemann]

Continuation of #1434. @thht I agree this is a good idea.

[agramfort]

make it work for your needs, use it, iterate until you're happy with API
and we'll make it nice to be integrated into MNE

my 2c

>

[NeuroLaunch]

one issue might be that when doing multitaper, you might end up with a non-uniform amount of tapers for the frequencies. no idea how this is implemented in mne, but fieldtrip offers to calculate the amount of tapers individually per frequency...

If as @larsoner suggests for data storage as (n_channels, n_frequencies, n_estimates), then perhaps the class could hold a function (standard or user-customized) that acts on all stored estimates/tapers to produce a desired final spectrum for each channel. A 1-D spectrum plot method, for instance, would invoke that aggregation function.

[dominikwelke]

hi there,

I am also following this discussion, since I'm interested in analysing frequency tagging data using mne.
I met @kalenkovich via the mne mailing list and we started discussing working on a public facing package for that.

Anyway, on the way he found ssvepy that does exactly this. It consists of a class carrying frequency spectra (for epoched data - simply the 3d numpy array spit out by the various psd methods), and features some plot methods, plus specific frequency tagging stuff, and tfr.

I don't think it is optimal for the purpose discussed here - it lacks flexibility, and the specific frequency-tagging functionality might rather go to a separate class - but its nice to get an impression how something like this could look and feel. and might help to avoid some mistakes :)

Cheers, Dominik

[mmagnuski]

If anyone's interested I have a simple mne-compatible psd class implementation here:
https://github.com/mmagnuski/borsar/blob/master/borsar/freq.py#L244

It allows for things like:

# metadata selection and epochs averaging
psd_fast = psd['RT < 0.5'].average()
# window of interest selection just like in Evokeds
psd_fast.pick_channels('Oz').crop(fmin=8, fmax=12).data.mean()

# familiar mne plots
psd.plot(fmax=15, dB=True)
psd.plot_topomap(freqs=[6, 8, 10, 12])
psd.plot_joint(freqs=[12, 23])

I even have a hack to use it in mne.viz.plot_compare_evokeds:

dct = {'cond1': psd[cond1].to_evoked(), 'cond2': psd[cond2].to_evoked()}
mne.viz.plot_compare_evokeds(dct, ...)

[mmagnuski]

my implementation does not store n_estimates but it could be easily extended for that.