mne-tools
Current multitaper implementation too slow for continuous raw data
Describe the new feature or enhancement
Currently, it is basically impossible to use raw.compute_psd(method="multitaper", bandwidth=1). For long data, it calculates hundreds of DPSS windows and therefore it takes way too long. In the example below, I used bandwidth=1 which is one of the fastest settings but it takes 6 minutes for a single raw array with shape (376 x 277s * 600Hz) (taken from this tutorial). I did not have time to test how long it takes for bandwidth=4 but probably more than 20 minutes on an M1 MacBook Pro.
Describe your proposed implementation
I am not very familiar with multitapers. However, it should be really easy to speed up the implementation, by simply calculating the multitapers on epochs of the raw data and averaging afterwards.
I discussed this at the CuttingMEEG conference with @britta-wstnr and I think she agreed.
Describe possible alternatives
Multitaper too slow: Jupyter Notebook
Additional context
No response
Hi @moritz-gerster, nice to see you on Github! :)
As a little clarification, I meant that it would be possible to approximate the power in the continuous signal by cutting the signal into smaller segments e.g. using make_fixed_length_epochs, and then running multitaper on the resulting Epochs.
I have a faint memory of a discussion about this or a similar issue within the last two years or so, but did not remember it off the top of my head as I never work with continuous data myself.
I meant that it would be possible to approximate the power in the continuous signal by cutting the signal into smaller segments e.g. using
make_fixed_length_epochs, and then running multitaper on the resulting Epochs
Right @britta-wstnr, that's exactly what I do in the attached jupyter notebook. The approximation looks very similar to the continuous input but its 60 times faster. I was wondering whether mne shouldn't do this by default.
As a little clarification, I meant that it would be possible to approximate the power in the continuous signal by cutting the signal into smaller segments e.g. using make_fixed_length_epochs, and then running multitaper on the resulting Epochs.
To me this sounds like a combined Welch + multitaper PSD, right? One way to accomplish this would be to still use the n_per_seg and n_overlap option in 'multitaper' mode to say how to chop up the Raw data before computing the multitaper PSD. This is probably what most people would want anyway from mode='multitaper', and if they want the old behavior they can always set n_per_seg=len(raw.times) and live with a computation that takes forever. I think we could consider this a bugfix, and I doubt many users would be affected since the computation is so slow it's probably not widely used (or used at all) in its current form. WDYT @drammock ?
before the addition of the Spectrum class, we didn't support multitaper on Raw at all, so I agree that it's not likely to be a widespread use case yet. Using the n_per_seg and n_overlap seems like a good solution. Is this worth blocking the 1.3 release?
Hehe, sorry @moritz-gerster, I admit I did not look at your Notebook, as I am technically on vacation ;-) sorry about that! Looks like there is a good idea how to move on this - @moritz-gerster would you be up for giving this a try?
would you be up for giving this a try?
I just had a look because I thought it could be simple. However, I understand both the functions psd_array_multitaper and psd_array_welch very badly. I'm sorry, I don't feel fit for it :/
I hope you find someone more competent. Thanks again for the exciting lecture and enjoy your vacation @britta-wstnr!