MRG, BUG, ENH: Add window option

doc/changes/latest.inc

@@ -83,6 +83,8 @@ Enhancements
 
 - Add option ``on_header_missing`` to :func:`mne.channels.read_polhemus_fastscan` (:gh:`8622` by `Eric Larson`_)
 
+- Add option ``window`` to :func:`mne.time_frequency.psd_welch` and related functions (:gh:`8862` by `Eric Larson`_)
+
 - `mne.preprocessing.ICA.plot_sources` now displays an `mne.preprocessing.ICA.plot_properties` window when right-clicking on component names on the y-axis (:gh:`8381` by `Daniel McCloy`_)
 
 - :func:`mne.io.read_raw_edf`, :func:`mne.io.read_raw_bdf`, and :func:`mne.io.read_raw_gdf` now detect and handle invalid highpass/lowpass filter settings (:gh:`8584` by `Clemens Brunner`_)
@@ -173,6 +175,8 @@ Bugs
 
 - Fix bug in `mne.viz.plot_compare_evokeds` where evokeds with identical ``comment`` attributes would not plot properly if passed as a list (:gh:`8590` by `Daniel McCloy`_)
 
+- Fix bug in :func:`mne.time_frequency.psd_welch` and related functions where the window default errantly changed from ``'hamming'`` to ``('tukey', 0.25)`` (:gh:`8862` by `Eric Larson`_)
+
 - Fix bug in :func:`mne.io.read_raw_kit` where scale factors for EEG channels could be set to zero (:gh:`8542` by `Eric Larson`_)
 
 - Fix reading GDF files with excluded channels in :func:`mne.io.read_raw_gdf` (:gh:`8520` by `Clemens Brunner`_)

mne/io/base.py

@@ -1443,7 +1443,8 @@ def plot_psd(self, fmin=0, fmax=np.inf, tmin=None, tmax=None, proj=False,
                  picks=None, ax=None, color='black', xscale='linear',
                  area_mode='std', area_alpha=0.33, dB=True, estimate='auto',
                  show=True, n_jobs=1, average=False, line_alpha=None,
-                 spatial_colors=True, sphere=None, verbose=None):
+                 spatial_colors=True, sphere=None, window='hamming',
+                 verbose=None):
         return plot_raw_psd(self, fmin=fmin, fmax=fmax, tmin=tmin, tmax=tmax,
                             proj=proj, n_fft=n_fft, n_overlap=n_overlap,
                             reject_by_annotation=reject_by_annotation,
@@ -1452,7 +1453,7 @@ def plot_psd(self, fmin=0, fmax=np.inf, tmin=None, tmax=None, proj=False,
                             dB=dB, estimate=estimate, show=show, n_jobs=n_jobs,
                             average=average, line_alpha=line_alpha,
                             spatial_colors=spatial_colors, sphere=sphere,
-                            verbose=verbose)
+                            window=window, verbose=verbose)
 
     @copy_function_doc_to_method_doc(plot_raw_psd_topo)
     def plot_psd_topo(self, tmin=0., tmax=None, fmin=0, fmax=100, proj=False,

mne/time_frequency/psd.py

@@ -84,7 +84,8 @@ def _check_psd_data(inst, tmin, tmax, picks, proj, reject_by_annotation=False):
 
 @verbose
 def psd_array_welch(x, sfreq, fmin=0, fmax=np.inf, n_fft=256, n_overlap=0,
-                    n_per_seg=None, n_jobs=1, average='mean', verbose=None):
+                    n_per_seg=None, n_jobs=1, average='mean', window='hamming',
+                    verbose=None):
     """Compute power spectral density (PSD) using Welch's method.
 
     Parameters
@@ -107,13 +108,12 @@ def psd_array_welch(x, sfreq, fmin=0, fmax=np.inf, n_fft=256, n_overlap=0,
         Length of each Welch segment (windowed with a Hamming window). Defaults
         to None, which sets n_per_seg equal to n_fft.
     %(n_jobs)s
-    average : str | None
-        How to average the segments. If ``mean`` (default), calculate the
-        arithmetic mean. If ``median``, calculate the median, corrected for
-        its bias relative to the mean. If ``None``, returns the unaggregated
-        segments.
+    %(average-psd)s
 
         .. versionadded:: 0.19.0
+    %(window-psd)s
+
+        .. versionadded:: 0.22.0
     %(verbose)s
 
     Returns
@@ -154,11 +154,14 @@ def psd_array_welch(x, sfreq, fmin=0, fmax=np.inf, n_fft=256, n_overlap=0,
 
     # Parallelize across first N-1 dimensions
     x_splits = np.array_split(x, n_jobs)
+    logger.debug(
+        f'Spectogram using {n_fft}-point FFT on {n_per_seg} samples with '
+        f'{n_overlap} overlap and {window} window')
 
     from scipy.signal import spectrogram
     parallel, my_spect_func, n_jobs = parallel_func(_spect_func, n_jobs=n_jobs)
     func = partial(spectrogram, noverlap=n_overlap, nperseg=n_per_seg,
-                   nfft=n_fft, fs=sfreq)
+                   nfft=n_fft, fs=sfreq, window=window)
     f_spect = parallel(my_spect_func(d, func=func, freq_sl=freq_sl,
                                      average=average)
                        for d in x_splits)
@@ -173,7 +176,8 @@ def psd_array_welch(x, sfreq, fmin=0, fmax=np.inf, n_fft=256, n_overlap=0,
 @verbose
 def psd_welch(inst, fmin=0, fmax=np.inf, tmin=None, tmax=None, n_fft=256,
               n_overlap=0, n_per_seg=None, picks=None, proj=False, n_jobs=1,
-              reject_by_annotation=True, average='mean', verbose=None):
+              reject_by_annotation=True, average='mean', window='hamming',
+              verbose=None):
     """Compute the power spectral density (PSD) using Welch's method.
 
     Calculates periodograms for a sliding window over the time dimension, then
@@ -209,13 +213,12 @@ def psd_welch(inst, fmin=0, fmax=np.inf, tmin=None, tmax=None, n_fft=256,
     %(reject_by_annotation_raw)s
 
         .. versionadded:: 0.15.0
-    average : str | None
-        How to average the segments. If ``mean`` (default), calculate the
-        arithmetic mean. If ``median``, calculate the median, corrected for
-        its bias relative to the mean. If ``None``, returns the unaggregated
-        segments.
+    %(average-psd)s
 
         .. versionadded:: 0.19.0
+    %(window-psd)s
+
+        .. versionadded:: 0.22.0
     %(verbose)s
 
     Returns
@@ -246,7 +249,8 @@ def psd_welch(inst, fmin=0, fmax=np.inf, tmin=None, tmax=None, n_fft=256,
                                   reject_by_annotation=reject_by_annotation)
     return psd_array_welch(data, sfreq, fmin=fmin, fmax=fmax, n_fft=n_fft,
                            n_overlap=n_overlap, n_per_seg=n_per_seg,
-                           average=average, n_jobs=n_jobs, verbose=verbose)
+                           average=average, n_jobs=n_jobs, window=window,
+                           verbose=verbose)
 
 
 @verbose

mne/time_frequency/tests/test_psd.py

@@ -6,7 +6,7 @@
 
 from mne import pick_types, Epochs, read_events
 from mne.io import RawArray, read_raw_fif
-from mne.utils import run_tests_if_main
+from mne.utils import catch_logging
 from mne.time_frequency import psd_welch, psd_multitaper, psd_array_welch
 
 base_dir = op.join(op.dirname(__file__), '..', '..', 'io', 'tests', 'data')
@@ -30,6 +30,12 @@ def test_psd_nan():
         x[0], float(n_fft), n_fft=n_fft, n_overlap=n_overlap)
     assert_allclose(freqs, freqs_2)
     assert_allclose(psds[0], psds_2)
+    # defaults
+    with catch_logging() as log:
+        psd_array_welch(x, float(n_fft), verbose='debug')
+    log = log.getvalue()
+    assert 'using 256-point FFT on 256 samples with 0 overlap' in log
+    assert 'hamming window' in log
 
 
 def test_psd():
@@ -61,7 +67,15 @@ def test_psd():
     for func, kws in funcs:
         kws = kws.copy()
         kws.update(kws_psd)
-        psds, freqs = func(raw, proj=False, **kws)
+        kws.update(verbose='debug')
+        if func is psd_welch:
+            kws.update(window='hann')
+        with catch_logging() as log:
+            psds, freqs = func(raw, proj=False, **kws)
+        log = log.getvalue()
+        if func is psd_welch:
+            assert f'{n_fft}-point FFT on {n_fft} samples with 0 overl' in log
+            assert 'hann window' in log
         psds_proj, freqs_proj = func(raw, proj=True, **kws)
 
         assert psds.shape == (len(kws['picks']), len(freqs))
@@ -264,6 +278,3 @@ def test_compares_psd():
 
     assert (np.sum(psds_welch < 0) == 0)
     assert (np.sum(psds_mpl < 0) == 0)
-
-
-run_tests_if_main()

mne/utils/docs.py

@@ -578,6 +578,17 @@
     Frequency-domain window to use in resampling.
     See :func:`scipy.signal.resample`.
 """
+docdict['average-psd'] = """
+average : str | None
+    How to average the segments. If ``mean`` (default), calculate the
+    arithmetic mean. If ``median``, calculate the median, corrected for
+    its bias relative to the mean. If ``None``, returns the unaggregated
+    segments.
+"""
+docdict['window-psd'] = """
+window : str | float | tuple
+    Windowing function to use. See :func:`scipy.signal.get_window`.
+"""
 docdict['decim'] = """
 decim : int
     Factor by which to subsample the data.

mne/viz/_figure.py

@@ -2297,7 +2297,7 @@ def _line_figure(inst, axes=None, picks=None, **kwargs):
 
 def _psd_figure(inst, proj, picks, axes, area_mode, tmin, tmax, fmin, fmax,
                 n_jobs, color, area_alpha, dB, estimate, average,
-                spatial_colors, xscale, line_alpha, sphere, **kwargs):
+                spatial_colors, xscale, line_alpha, sphere, window, **kwargs):
     """Instantiate a new power spectral density figure."""
     from .. import BaseEpochs
     from ..io import BaseRaw
@@ -2309,7 +2309,7 @@ def _psd_figure(inst, proj, picks, axes, area_mode, tmin, tmax, fmin, fmax,
         if kw in kwargs:
             psd_kwargs[kw] = kwargs.pop(kw)
     if isinstance(inst, BaseRaw):
-        psd_func = psd_welch
+        psd_func = partial(psd_welch, window=window)
     elif isinstance(inst, BaseEpochs):
         psd_func = psd_multitaper
     else:

mne/viz/epochs.py

@@ -965,12 +965,15 @@ def plot_epochs_psd(epochs, fmin=0, fmax=np.inf, tmin=None, tmax=None,
     from ._figure import _psd_figure
 
     # generate figure
+    # epochs always use multitaper, not Welch, so no need to allow "window"
+    # param above
     fig = _psd_figure(
         inst=epochs, proj=proj, picks=picks, axes=ax, tmin=tmin, tmax=tmax,
         fmin=fmin, fmax=fmax, sphere=sphere, xscale=xscale, dB=dB,
         average=average, estimate=estimate, area_mode=area_mode,
         line_alpha=line_alpha, area_alpha=area_alpha, color=color,
         spatial_colors=spatial_colors, n_jobs=n_jobs, bandwidth=bandwidth,
-        adaptive=adaptive, low_bias=low_bias, normalization=normalization)
+        adaptive=adaptive, low_bias=low_bias, normalization=normalization,
+        window='hamming')
     plt_show(show)
     return fig

mne/viz/raw.py

@@ -374,7 +374,8 @@ def plot_raw_psd(raw, fmin=0, fmax=np.inf, tmin=None, tmax=None, proj=False,
                  picks=None, ax=None, color='black', xscale='linear',
                  area_mode='std', area_alpha=0.33, dB=True, estimate='auto',
                  show=True, n_jobs=1, average=False, line_alpha=None,
-                 spatial_colors=True, sphere=None, verbose=None):
+                 spatial_colors=True, sphere=None, window='hamming',
+                 verbose=None):
     """%(plot_psd_doc)s.
 
     Parameters
@@ -414,6 +415,9 @@ def plot_raw_psd(raw, fmin=0, fmax=np.inf, tmin=None, tmax=None, proj=False,
     %(plot_psd_line_alpha)s
     %(plot_psd_spatial_colors)s
     %(topomap_sphere_auto)s
+    %(window-psd)s
+
+        .. versionadded:: 0.22.0
     %(verbose)s
 
     Returns
@@ -435,7 +439,8 @@ def plot_raw_psd(raw, fmin=0, fmax=np.inf, tmin=None, tmax=None, proj=False,
         average=average, estimate=estimate, area_mode=area_mode,
         line_alpha=line_alpha, area_alpha=area_alpha, color=color,
         spatial_colors=spatial_colors, n_jobs=n_jobs, n_fft=n_fft,
-        n_overlap=n_overlap, reject_by_annotation=reject_by_annotation)
+        n_overlap=n_overlap, reject_by_annotation=reject_by_annotation,
+        window=window)
     plt_show(show)
     return fig
 

mne/viz/tests/test_figure.py

@@ -15,4 +15,4 @@ def test_browse_figure_constructor():
 def test_psd_figure_constructor():
     """Test error handling in MNELineFigure constructor."""
     with pytest.raises(TypeError, match='an instance of Raw or Epochs, got'):
-        _psd_figure('foo', *((None,) * 18))
+        _psd_figure('foo', *((None,) * 19))