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:
|
# 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:
|
# 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:
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:
|
# 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:
|
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
Link to data
No response
Expected results
N/A
Actual results
N/A
Additional information
N/A
Description of the problem
When working through some bugs with plotting TFR multitaper data in #12910 I was looking over the
_time_frequency_loopcode 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
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
... and summed together:
mne-python/mne/time_frequency/tfr.py
Line 740 in b329515
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
... 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
@larsoner do you have any ideas?
Steps to reproduce
Link to data
No response
Expected results
N/A
Actual results
N/A
Additional information
N/A