Multitaper keyword "bandwidth" does not behave as expected

[moritz-gerster]

Proposed documentation enhancement

I try to calculate power spectral densities using multitaper at different frequency resolutions but I don't get it to work. The documentation states

bandwidth[float]
The bandwidth of the multi taper windowing function in Hz. The default value is a window half-bandwidth of 4.

I understand this doc such that mne.time_frequency.psd_multitaper(raw, bandwidth=None)

should yield the same as mne.time_frequency.psd_multitaper(raw, bandwidth=4)

but it does not. Looking at the code, I see how half_nbw is obtained from bandwidth: https://github.com/mne-tools/mne-python/blob/2d90ecaee807d3d7ceab31bdac0a734bf5418c87/mne/time_frequency/multitaper.py#L341

So I try mne.time_frequency.psd_multitaper(raw, bandwidth=(4 * raw.last_samp) / (2 * raw.info["sfreq"])) which again yields something different.

My question: Which float do I need to pass to bandwidth to reproduce bandwith=None?

Here is my code:

import mne
import numpy as np
import matplotlib.pyplot as plt

info = mne.create_info(["eeg1"], sfreq=256, ch_types="eeg")
data = data = np.random.randn(1, 256*10)
raw = mne.io.RawArray(data, info)

bandwidths = [None, 
              4, 
              (4 * raw.last_samp) / (2 * raw.info["sfreq"])]

_, axes = plt.subplots(1, 3, figsize=(10, 5), sharey=True)
for ax, bandwidth in zip(axes, bandwidths):
    raw.compute_psd(method="multitaper", bandwidth=bandwidth).plot(axes=ax)
    ax.set_xlabel("Frequency (Hz)")
    ax.set_title(f"Bandwidth={bandwidth}")
# plt.savefig("bandwidths.png", dpi=150)
plt.show()

This is the result:

[mscheltienne]

Hello, have also a look at the updated documentation in 1.3 which is a bit more detailed: https://mne.tools/dev/generated/mne.time_frequency.psd_array_multitaper.html

But there seems to be a bug here. The value that yields the same as bandwidth=None in your example is 0.8. I don't get the multiplication with n_times / (2. * sfreq) here https://github.com/mne-tools/mne-python/blob/815bee6b38c160e6f83d86bf929e79040536c83d/mne/time_frequency/multitaper.py#L320-L323

EDIT: It looks like if you provide a float, it gets converted to the normalized bandwidth; while when you provide None it is not normalized. I guess the same multiplication should be present in the case bandwidth=None?

EDIT2: The normalization looks also weird.. bandwidth=4 yields half_nbw=20.

[moritz-gerster]

But there seems to be a bug here. The value that yields the same as bandwidth=None in your example is 0.8.

I guess the label of the issue should be switched from "DOC" to "BUG" but I cannot do that.

[mmagnuski]

@moritz-gerster @mscheltienne I've recently noticed this, the bandwidth that is used in psd is the frequency bandwidth, not the time-bandwidth product, as in time-frequency. In general I still think it would be simpler to expose something like frequency smoothing (frequency bandwidth / 2), like fieldtrip does.

[mmagnuski]

I didn't check in detail but I think the default bandwidth of None chooses bandwidth and then time-bandwidth product that gives 7 tapers. In your case the signal is 10 seconds long, so a bandwidth of 0.8 gives time-bandwidth product of 8 (0.8 * 10), and we subtract 1 when we calculate the number of tapers which gives the default 7. A default of 7 tapers is not great, I admit, we may think about changing that in the future to a constant bandwidth - but that has the downside of raising errors on short signals for example, when default bandwidth would lead to < 1 taper.

[larsoner]

that has the downside of raising errors on short signals for example, when default bandwidth would lead to < 1 taper

We could just say that the default is a constant bandwidth X, unless that yields no tapers in which case the smallest bandwidth is chosen that still yields a single taper.

[mmagnuski]

(I changed the label back to DOC, as don't think this is a bug, but I agree that the inconsistency in psd and time-frequency multitaper is source of many confusions)

[moritz-gerster]

I didn't check in detail but I think the default bandwidth of None chooses bandwidth and then time-bandwidth product that gives 7 tapers. In your case the signal is 10 seconds long, so a bandwidth of 0.8 gives time-bandwidth product of 8 (0.8 * 10), and we subtract 1 when we calculate the number of tapers which gives the default 7. A default of 7 tapers is not great, I admit, we may think about changing that in the future to a constant bandwidth - but that has the downside of raising errors on short signals for example, when default bandwidth would lead to < 1 taper.

Yes, 0.8 is the solution.

[larsoner]

@moritz-gerster @mmagnuski @mscheltienne @drammock you all are working/thinking about this the most it seems, so happy to go with whatever you think is the most reasonable behavior here. It seems like at the very least, a default change (which could be considered a bugfix I think) as suggested by @mmagnuski above would probably make sense, anything else? Do we also have some doc problem and/or bug with the actual computations to fix?

[mmagnuski]

One issue I think is worth fixing is the discrepancy between psd and time-frequency multitaper usage of bandwidth. psd uses frequency bandwidth, while time-frequency multitaper uses time-bandwidth product (full or half, I don't remember), which is less intuitive. It would be best if both used frequency bandwidth, but I'm not sure this can classify as bug fix.

[tomdstone]

https://github.com/mne-tools/mne-python/blob/8dfb23be27536d654838377622a898b08f86cd31/mne/time_frequency/multitaper.py#L291-L292

Doesn't this show that the bandwidth variable is assumed to be the full frequency bandwidth? The documentation for the psd function says it is the half bandwidth, so could part of the issue be that when we plug back in the half bandwidth it then halves it again?

https://github.com/mne-tools/mne-python/blob/8dfb23be27536d654838377622a898b08f86cd31/mne/time_frequency/multitaper.py#L333-L336

[mmagnuski]

@tomdstone Yes, that should be the full bandwidth. The docs changed in #11265, but I didn't noticed then that it went from bandwidth to half-bandwidth. Feel free to open a PR to fix the docs.

The next step would be to change the default bandwidth in psd as @larsoner mentioned. Then it would be a good idea to make tfr multitaper also use bandwidth instead of time-bandwidth product. But this would have to be a deprecation, so it requires some careful planning.

[moritz-gerster]

🥳