is this a bug in the new `compute_psd` method?

💬 Support & Discussions
1

In some cases, I’m getting an assertion error with the new compute_psd method. Here are some examples of what does and does not work for me.

import mne


misc_path = mne.datasets.misc.data_path()
raw = mne.io.read_raw(misc_path / 'seeg' / 'sample_seeg_ieeg.fif')
n_per_seg = int(raw.info['sfreq'])
fmax = raw.info['sfreq']//2
overlap = int((0.25)*raw.info['sfreq'] + 0.5)
print(n_per_seg, fmax, overlap)

spectrum = raw.compute_psd(method='welch')
psd = spectrum.get_data()
print(psd.shape)

spectrum2 = raw.compute_psd(method='welch', average=None)
psd2 = spectrum.get_data()
print(psd2.shape)

spectrum3 = raw.compute_psd(method='welch', fmin=0, fmax=fmax,
                            n_fft=n_per_seg,
                            n_per_seg=n_per_seg,
                            n_overlap=overlap)
psd3 = spectrum.get_data()
print(psd3.shape)


spectrum4 = raw.compute_psd(method='welch', fmin=0, fmax=fmax,
                            n_fft=n_per_seg,
                            n_per_seg=n_per_seg,
                            n_overlap=overlap, average=None)
psd4 = spectrum4.get_data()
print(psd4.shape)

the spectrum4 calculation give the following error

Traceback (most recent call last):
  File "c:\Users\leisenman\Box\Box_Documents\Python\test_mne_1_2.py", line 27, in <module>
    spectrum4 = raw.compute_psd(method='welch', fmin=0, fmax=fmax,
  File "<decorator-gen-238>", line 12, in compute_psd
  File "C:\Users\leisenman\Miniconda3\envs\seeg\lib\site-packages\mne\io\base.py", line 1841, in compute_psd
    return Spectrum(
  File "C:\Users\leisenman\Miniconda3\envs\seeg\lib\site-packages\mne\time_frequency\spectrum.py", line 944, in __init__
    self._check_values()
  File "C:\Users\leisenman\Miniconda3\envs\seeg\lib\site-packages\mne\time_frequency\spectrum.py", line 362, in _check_values
    assert self._data.shape == self._shape
AssertionError

Am I doing something wrong or is this a bug? Thanks.

Larry

Platform:         Windows-10-10.0.19042-SP0
Python:           3.10.6 | packaged by conda-forge | (main, Aug 22 2022, 20:30:19) [MSC v.1929 64 bit (AMD64)]
Executable:       C:\Users\leisenman\Miniconda3\envs\seeg\python.exe
CPU:              Intel64 Family 6 Model 140 Stepping 1, GenuineIntel: 8 cores
Memory:           15.7 GB

mne:              1.2.0
numpy:            1.23.3 {MKL 2022.1-Product with 4 threads}
scipy:            1.9.1
matplotlib:       3.6.1 {backend=QtAgg}

sklearn:          1.1.2
numba:            Not found
nibabel:          4.0.2
nilearn:          0.9.2
dipy:             Not found
openmeeg:         Not found
cupy:             Not found
pandas:           1.5.0
pyvista:          0.36.1 {OpenGL 4.5.0 - Build 30.0.100.9955 via Intel(R) Iris(R) Xe Graphics}
pyvistaqt:        0.9.0
ipyvtklink:       0.2.2
vtk:              9.1.0
qtpy:             2.2.1 {PyQt5=5.15.4}
ipympl:           Not found
pyqtgraph:        Not found
pooch:            v1.6.0

mne_bids:         0.11
mne_nirs:         Not found
mne_features:     Not found
mne_qt_browser:   Not found
mne_connectivity: Not found
mne_icalabel:     Not found
2

that is a bug in MNE. As a safety check, we try to predict how many Welch segments will result and it seems we are not computing it correctly, it is off-by-one. I’m working on a fix and will post back here when I’ve opened a Pull Request that fixes it.

1 Like
3

fix is in: Fix welch windows by drammock · Pull Request #11248 · mne-tools/mne-python · GitHub

1 Like
4

the fix has been merged into current main branch and has also been backported to maint/1.2, so it will be included in the patch release 1.2.1 (probably later this week).

1 Like
5

That did the trick. Thanks for the rapid reply and fix.

Larry

1 Like