[ENH, MRG] Add compute_csd function

[alexrockhill]

Part of #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.

[alexrockhill]

Failed tests are because this depends on #10979

[alexrockhill]

@wmvanvliet if you have a minute to review this one too that would be really great :)

[alexrockhill]

Needs to check that it gives the same results as csd_morlet etc., I think there's a scaling issue.

[alexrockhill]

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 takes n_cycles whereas csd_multitaper takes an fmin and fmax), @drammock do you know if there's a backward compatible way to match the fmin/fmax parameterization with the freqs/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 takes epochs 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 match csd_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.

[alexrockhill]

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.

[alexrockhill]

Okay to merge? This follows the definition of CSD here

(From @britta-wstnr et al. 's 2022 Neuroimage paper and pointed out by her)

[agramfort]

@drammock @larsoner @wmvanvliet ok for you?

[larsoner]

Just marking as changes requested until my comments above are addressed. CI failure is unrelated and should be fixed by #11034 eventually

[alexrockhill]

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.

[alexrockhill]

Looks good to go unless I missed anything

[larsoner]

Thanks @alexrockhill !