mne-python
mne-python copied to clipboard
[ENH, MRG] Add compute_csd function
Part of https://github.com/mne-tools/mne-python/pull/10920.
This uses power to compute the cross-spectral density. It would be nice if this could be done on complex numbers but that's causing issues with positive semi-definite matrices. It might be nice to just go ahead with this version and then improve it later since it matches how other functions like csd_morlet
do CSD and then modify all the functions to be consistent in a followup PR.
Failed tests are because this depends on https://github.com/mne-tools/mne-python/pull/10979
@wmvanvliet if you have a minute to review this one too that would be really great :)
Needs to check that it gives the same results as csd_morlet
etc., I think there's a scaling issue.
Okay, so this passes tests locally that calling tfr_morlet(epochs)
and then compute_csd(epochs_tfr)
is the same as csd_morlet(epochs)
by changing it to match. But, a couple important questions remain:
- Does this code work for
csd_multitaper
(csd_fourier
I think is not applicable)? I wasn't able to check because the parameterization is currently different (tfr_multitaper
takesn_cycles
whereascsd_multitaper
takes anfmin
andfmax
), @drammock do you know if there's a backward compatible way to match thefmin
/fmax
parameterization with thefreqs
/n_cycles
parameterization? - Would this theoretically work with Hilbert transformed data? Maybe this is an empirical question to compare with the filtered raw data approach, tagging @britta-wstnr as I think she might be interested
- There are many different methods in
mne.compute_covariance
, can we allow those as well? It takesepochs
with dimension (n_epochs
,n_channels
,n_times
) so it would need to take 4D input or somehow handle frequencies. The way that I had it in the previous commits returned something sensible but didn't matchcsd_morlet
even when normalized so I think there's some scaling going on that may make things different but it would be interesting to know if those give better results.
I think question 1 should be addressed before merging but questions 2 and 3 are more long-term trying to figure out how best to do this kind of thing.
I took a bit more of a look at it appears that tfr_mutlitaper
doesn't even accept an output
argument to get complex data. This was easy to add (you have to average across the tapers which I'm not sure works really well for phase though) but did not give quantitatively similar results using the same frequencies (you still have the freqs
+ n_cycles
parameterization vs the fmin
/fmax
parameterization) so I think it does not work. Up to people who maintain csd
functions more whether to go ahead and merge knowing that it only works for Morlet or add a warning that it only works for EpochsTFR
derived from tfr_morlet
or dig into the multitaper more or what. This is a nice efficiency improvement for applying DICS to time-frequency epochs though so it would be nice to get merged at some point.
Okay to merge? This follows the definition of CSD here
(From @britta-wstnr et al. 's 2022 Neuroimage paper and pointed out by her)
Thanks for the review @larsoner, I added those changes. Now, https://github.com/mne-tools/mne-python/blob/main/mne/time_frequency/csd.py#L1344-L1345 is done the old none-einsum way (thanks for figuring that out by the way, I was pretty close actually but just didn't get the right combination of the first index being the only different one, I had the first and second and both the same but now that you wrote the correct one, it seems so easy) but it's done per channel so it's a bit more refactoring than should probably go in this PR.
Looks good to go unless I missed anything
Thanks @alexrockhill !