Possible bug in computation of multitaper TFR power

[tsbinns]

Description of the problem

When working through some bugs with plotting TFR multitaper data in #12910 I was looking over the _time_frequency_loop code for converting complex multitaper coefficients into real-valued power.

Here the timeseries data, X, and windows, Ws, for each taper, W, are convolved to give the time-frequency transformed data, coefs:

mne-python/mne/time_frequency/tfr.py

Lines 709 to 713 in b329515

	# Loops across tapers.
	for taper_idx, W in enumerate(Ws):
	# No need to check here, it's done earlier (outside parallel part)
	nfft = _get_nfft(W, X, use_fft, check=False)
	coefs = _cwt_gen(X, W, fsize=nfft, mode=mode, decim=decim, use_fft=use_fft)

When power is requested, the complex coefficients for each taper are converted to power:

mne-python/mne/time_frequency/tfr.py

Lines 719 to 723 in b329515

	# Loop across epochs
	for epoch_idx, tfr in enumerate(coefs):
	# Transform complex values
	if output in ["power", "avg_power"]:
	tfr = (tfr * tfr.conj()).real # power

... and summed together:

mne-python/mne/time_frequency/tfr.py

Line 740 in b329515

tfrs[epoch_idx] += tfr

But I don't see at any point a weighting of the coefficients according to the taper weights. There is at the end a normalisation according to the number of tapers:

mne-python/mne/time_frequency/tfr.py

Lines 752 to 754 in b329515

	# Normalization by number of taper
	if n_tapers > 1 and output not in ["complex", "phase"]:
	tfrs /= n_tapers

... but that just means each taper is treated as having an equal weighting.

I also didn't notice anywhere upstream of the TFR computation where the taper weights were already applied.

This could be a misunderstanding on my part, but at least based on how PSD computation is handled, taper weights are applied when converting complex coefficients to power:

mne-python/mne/time_frequency/multitaper.py

Lines 187 to 206 in b329515

	def _psd_from_mt(x_mt, weights):
	"""Compute PSD from tapered spectra.
	Parameters
	----------
	x_mt : array, shape=(..., n_tapers, n_freqs)
	Tapered spectra
	weights : array, shape=(n_tapers,)
	Weights used to combine the tapered spectra
	Returns
	-------
	psd : array, shape=(..., n_freqs)
	The computed PSD
	"""
	psd = weights * x_mt
	psd *= psd.conj()
	psd = psd.real.sum(axis=-2)
	psd *= 2 / (weights * weights.conj()).real.sum(axis=-2)
	return psd

@larsoner do you have any ideas?

Steps to reproduce

N/A

Link to data

No response

Expected results

N/A

Actual results

N/A

Additional information

N/A

[larsoner]

Yeah just the fact that we call _make_dpss without the return_weights=True suggests that they aren't/can't be used. It seems like they should be applied to the Ws and then normalized properly at the end.

[tsbinns]

Okay, so #12910 makes the weights available in the class which could be used. Is it worth just fixing this there, or open a new PR? Don't want to overburden with a big PR review.

[larsoner]

We're cutting a release now-ish. After that let's merge #12910 then a separate PR (which you could start now based on #12910 if you don't mind some rebase/merge mess) let's fix the computation.

[tsbinns]

Sure, sounds good!