Skip to content

align FFT windows to good data spans in psd_array_welch#12536

Merged
larsoner merged 3 commits intomne-tools:mainfrom
drammock:fix-welch-nan
Apr 9, 2024
Merged

align FFT windows to good data spans in psd_array_welch#12536
larsoner merged 3 commits intomne-tools:mainfrom
drammock:fix-welch-nan

Conversation

@drammock
Copy link
Copy Markdown
Member

@drammock drammock commented Apr 9, 2024

closes #11413

@drammock drammock requested a review from cbrnr April 9, 2024 19:54
@larsoner larsoner merged commit a92d798 into mne-tools:main Apr 9, 2024
@larsoner
Copy link
Copy Markdown
Member

larsoner commented Apr 9, 2024

Thanks @drammock !

drammock
drammock commented Apr 9, 2024
View reviewed changes
mne/time_frequency/psd.py
Comment on lines +235 to +237
weights = [
split.shape[-1] for split in x_splits if split.shape[-1] >= n_per_seg
]
Copy link
Copy Markdown
Member Author

@drammock drammock Apr 9, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

was actually just wondering if this is in fact the most accurate weighting. Perhaps this is better?

        # weights reflect the number of welch windows in each span
        stepsize = n_per_seg - n_overlap
        weights = [(span.shape[-1] - n_overlap) // stepsize for span in x_splits]
        weights = np.array(weights)[np.array(weights) > 0]

Old code weights based on total samples in the span (as long as span > n_per_seg). New code weights based on number of welch windows that fit into the span, i.e., should reflect how much data in the span is actually used in the computation. If that seems correct, I can open a follow-up PR. cc @agramfort @larsoner

Copy link
Copy Markdown
Member

@larsoner larsoner Apr 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new way seems more correct to me

@cbrnr
Copy link
Copy Markdown
Contributor

cbrnr commented Apr 10, 2024

I just tested this with the MWE in #11413, and at least for me this example is not fixed. I now get a different error:

>>> raw.plot_psd()
NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Setting 1201 of 6007 (19.99%) samples to NaN, retaining 4806 (80.01%) samples.
Effective window size : 3.410 (s)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<decorator-gen-150>", line 5, in plot_psd
  File "<decorator-gen-149>", line 12, in plot_psd
  File "/Users/clemens/Projects/mne-python/mne/time_frequency/spectrum.py", line 139, in plot_psd
    return self.compute_psd(**init_kw).plot(**plot_kw)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<decorator-gen-334>", line 12, in compute_psd
  File "/Users/clemens/Projects/mne-python/mne/io/base.py", line 2230, in compute_psd
    return Spectrum(
           ^^^^^^^^^
  File "/Users/clemens/Projects/mne-python/mne/time_frequency/spectrum.py", line 1146, in __init__
    self._compute_spectra(data, fmin, fmax, n_jobs, method_kw, verbose)
  File "/Users/clemens/Projects/mne-python/mne/time_frequency/spectrum.py", line 445, in _compute_spectra
    result = self._psd_func(
             ^^^^^^^^^^^^^^^
  File "<decorator-gen-148>", line 12, in psd_array_welch
  File "/Users/clemens/Projects/mne-python/mne/time_frequency/psd.py", line 251, in psd_array_welch
    psds = agg_func(f_spect, axis=0)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/clemens/.pyenv/versions/mne/lib/python3.12/site-packages/numpy/lib/function_base.py", line 541, in average
    raise ValueError(
ValueError: Length of weights not compatible with specified axis.

This is the second plot (after adding an annotation). The first plot works (as previously).

@drammock drammock mentioned this pull request Apr 10, 2024
Merged
@drammock drammock deleted the fix-welch-nan branch April 10, 2024 14:28
@drammock
Copy link
Copy Markdown
Member Author

drammock commented Apr 10, 2024

I just tested this with the MWE in #11413, and at least for me this example is not fixed. I now get a different error

argh sorry @cbrnr I forgot to go back and verify with your MWE. Should be fixed now by #12538

larsoner added a commit to larsoner/mne-python that referenced this pull request Apr 15, 2024
@larsoner
Merge remote-tracking branch 'upstream/main' into rank
455a3dc 
* upstream/main: (50 commits)
  ENH: Improve OPM auditory dataset and example (mne-tools#12539)
  MAINT: Bump to latest pydata-sphinx-theme (mne-tools#12228)
  MRG: Simplify manual installation instructions a little by dropping explicit mention of (lib)mamba (mne-tools#12362)
  fix PSD weights handling when bad annotations present (mne-tools#12538)
  Fix phase loading (mne-tools#12537)
  align FFT windows to good data spans in psd_array_welch (mne-tools#12536)
  explicitly disallow multitaper in presence of bad annotations (mne-tools#12535)
  MAINT: Clean up PyVista contexts (mne-tools#12533)
  MAINT: Complete API change of ordered (mne-tools#12534)
  MAINT: Reinstall statsmodels and improve logging (mne-tools#12532)
  MAINT: Remove scipy.signal.morlet2 (mne-tools#12531)
  Update README badge links (mne-tools#12529)
  BUG: Fix bug with reading his_id from snirf (mne-tools#12526)
  [pre-commit.ci] pre-commit autoupdate (mne-tools#12524)
  Fix file format check in _check_eeglab_fname function (mne-tools#12523)
  MAINT: Reenable picard in pre testing (mne-tools#12525)
  MAINT: Bump to large resource class (mne-tools#12522)
  MAINT: Restore 2 jobs on Windows (mne-tools#12520)
  Add exclude_after_unique option to mne.io.read_raw_edf/read_raw_edf to search for exclude channels after making channel names unique (mne-tools#12518)
  Improve consistency of sensor types in code and documentation (mne-tools#12509)
  ...
@moritz-gerster moritz-gerster mentioned this pull request Dec 21, 2024
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@agramfort agramfort agramfort approved these changes

@cbrnr cbrnr Awaiting requested review from cbrnr

@mscheltienne mscheltienne Awaiting requested review from mscheltienne mscheltienne is a code owner

@larsoner larsoner Awaiting requested review from larsoner larsoner is a code owner

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

All NaNs in raw.plot_psd()

4 participants

@drammock @larsoner @cbrnr @agramfort