[DOC,BUG] Add note about FWHM to Morlet, fix norm (mne-tools#11353)

Co-authored-by: Eric Larson <[email protected]>
Co-authored-by: Daniel McCloy <[email protected]>

doc/changes/latest.inc

@@ -30,6 +30,7 @@ Enhancements
 - Add size information to the ``repr`` of :class:`mne.Report` (:gh:`11357` by `Eric Larson`_)
 - Add support for ``image_format='webp'`` to :class:`mne.Report` when using Matplotlib 3.6+, which can reduce file sizes by up to 50% compared to ``'png'``. The new default ``image_format='auto'`` will automatically use this format if it's available on the system (:gh:`11359` by `Eric Larson`_)
 - Add :func:`mne.beamformer.apply_dics_tfr_epochs` to apply a DICS beamformer to time-frequency resolved epochs (:gh:`11096` by `Alex Rockhill`_)
+- Add :func:`mne.time_frequency.fwhm` to determine the full-width half maximum for :func:`mne.time_frequency.morlet` (:gh:`11353` by `Britta Westner`_, `Daniel McCloy`_, and `Eric Larson`_)
 - Check whether head radius (estimated from channel positions) is correct when reading EEGLAB data with :func:`~mne.io.read_raw_eeglab` and :func:`~mne.read_epochs_eeglab`. If head radius is not within likely values, warn informing about possible units mismatch and the new ``montage_units`` argument (:gh:`11283` by `Mikołaj Magnuski`_).
 - Add support for a callable passed in ``combine`` for `mne.time_frequency.AverageTFR.plot` and `mne.time_frequency.AverageTFR.plot_joint` (:gh:`11329` by `Mathieu Scheltienne`_)
 

doc/conf.py

@@ -92,6 +92,7 @@
     'mne_substitutions',
     'newcontrib_substitutions',
     'gen_names',
+    'matplotlib.sphinxext.plot_directive',
     'sphinxcontrib.bibtex',
     'sphinx_copybutton',
     'sphinx_design',
@@ -916,6 +917,31 @@ def append_attr_meth_examples(app, what, name, obj, options, lines):
 htmlhelp_basename = 'mne-doc'
 
 
+# -- Options for plot_directive ----------------------------------------------
+
+# Adapted from SciPy
+plot_include_source = True
+plot_formats = [('png', 96)]
+plot_html_show_formats = False
+plot_html_show_source_link = False
+font_size = 13 * 72 / 96.0  # 13 px
+plot_rcparams = {
+    'font.size': font_size,
+    'axes.titlesize': font_size,
+    'axes.labelsize': font_size,
+    'xtick.labelsize': font_size,
+    'ytick.labelsize': font_size,
+    'legend.fontsize': font_size,
+    'figure.figsize': (6, 5),
+    'figure.subplot.bottom': 0.2,
+    'figure.subplot.left': 0.2,
+    'figure.subplot.right': 0.9,
+    'figure.subplot.top': 0.85,
+    'figure.subplot.wspace': 0.4,
+    'text.usetex': False,
+}
+
+
 # -- Options for LaTeX output ------------------------------------------------
 
 # Grouping the document tree into LaTeX files. List of tuples
@@ -978,6 +1004,8 @@ def reset_warnings(gallery_conf, fname):
         'Jupyter is migrating its paths to use standard',
         # PyQt6
         'Enum value .* is marked as deprecated',
+        # matplotlib PDF output
+        'The py23 module has been deprecated',
     ):
         warnings.filterwarnings(  # deal with other modules having bad imports
             'ignore', message=".*%s.*" % key, category=DeprecationWarning)

doc/references.bib

@@ -272,6 +272,16 @@ @book{Cohen2014
   year={2014}
 }
 
+@article{Cohen2019,
+  author={Cohen, Michael X},
+  doi = {10.1016/j.neuroimage.2019.05.048},
+  journal = {NeuroImage},
+  pages = {81-86},
+  title = {A better way to define and describe {Morlet} wavelets for time-frequency analysis},
+  volume = {199},
+  year = {2019}
+}
+
 @article{CohenHosaka1976,
   author = {Cohen, David and Hosaka, Hidehiro},
   doi = {10.1016/S0022-0736(76)80041-6},
@@ -1752,6 +1762,15 @@ @article{TadelEtAl2011
   year = {2011}
 }
 
+@article{Tallon-BaudryEtAl1997,
+	author = {Tallon-Baudry, Catherine and Bertrand, Olivier and Delpuech, Claude and Pernier, Jacques},
+  doi = {10.1523/JNEUROSCI.17-02-00722.1997},
+	journal = {Journal of Neuroscience},
+	pages = {722-734},
+	title = {Oscillatory {Gamma}-{Band} (30–70 {Hz}) {Activity} {Induced} by a {Visual} {Search} {Task} in {Humans}},
+	year = {1997},
+}
+
 @article{TannerEtAl2015,
   author = {Tanner, Darren and {Morgan-Short}, Kara and Luck, Steven J.},
   doi = {10.1111/psyp.12437},

doc/time_frequency.rst

@@ -47,6 +47,7 @@ Functions that operate on ``np.ndarray`` objects:
    csd_array_multitaper
    csd_array_morlet
    dpss_windows
+   fwhm
    morlet
    stft
    istft

mne/conftest.py

@@ -127,6 +127,7 @@ def pytest_configure(config):
     ignore:Passing a schema to Validator\.iter_errors is deprecated.*:
     ignore:Unclosed context <zmq.asyncio.Context.*:ResourceWarning
     ignore:Jupyter is migrating its paths.*:DeprecationWarning
+    ignore:Widget\..* is deprecated\.:DeprecationWarning
     # hopefully temporary https://github.com/matplotlib/matplotlib/pull/24455#issuecomment-1319318629
     ignore:The circles attribute was deprecated in Matplotlib.*:
     # PySide6

mne/time_frequency/__init__.py

@@ -1,7 +1,7 @@
 """Time frequency analysis tools."""
 
 from .tfr import (morlet, tfr_morlet, AverageTFR, tfr_multitaper, _BaseTFR,
-                  read_tfrs, write_tfrs, EpochsTFR, tfr_array_morlet)
+                  read_tfrs, write_tfrs, EpochsTFR, tfr_array_morlet, fwhm)
 from .psd import psd_array_welch
 from .csd import (CrossSpectralDensity, csd_fourier, csd_multitaper,
                   csd_morlet, csd_array_fourier, csd_array_multitaper,

mne/time_frequency/tests/test_tfr.py

@@ -7,6 +7,7 @@
 from numpy.testing import (assert_array_equal, assert_equal, assert_allclose)
 import pytest
 import matplotlib.pyplot as plt
+from scipy.signal import morlet2
 
 import mne
 from mne import (Epochs, read_events, pick_types, create_info, EpochsArray,
@@ -17,7 +18,7 @@
 from mne.time_frequency.tfr import (morlet, tfr_morlet, _make_dpss,
                                     tfr_multitaper, AverageTFR, read_tfrs,
                                     write_tfrs, combine_tfr, cwt, _compute_tfr,
-                                    EpochsTFR)
+                                    EpochsTFR, fwhm)
 from mne.time_frequency import tfr_array_multitaper, tfr_array_morlet
 from mne.viz.utils import _fake_click, _fake_keypress, _fake_scroll
 from mne.tests.test_epochs import assert_metadata_equal
@@ -38,13 +39,35 @@ def test_tfr_ctf():
         method(epochs, [10], 1)  # smoke test
 
 
-def test_morlet():
+@pytest.mark.parametrize('sfreq', [1000., 100 + np.pi])
+@pytest.mark.parametrize('freq', [10., np.pi])
+@pytest.mark.parametrize('n_cycles', [7, 2])
+def test_morlet(sfreq, freq, n_cycles):
     """Test morlet with and without zero mean."""
-    Wz = morlet(1000, [10], 2., zero_mean=True)
-    W = morlet(1000, [10], 2., zero_mean=False)
+    Wz = morlet(sfreq, freq, n_cycles, zero_mean=True)
+    W = morlet(sfreq, freq, n_cycles, zero_mean=False)
 
-    assert (np.abs(np.mean(np.real(Wz[0]))) < 1e-5)
-    assert (np.abs(np.mean(np.real(W[0]))) > 1e-3)
+    assert np.abs(np.mean(np.real(Wz))) < 1e-5
+    if n_cycles == 2:
+        assert np.abs(np.mean(np.real(W))) > 1e-3
+    else:
+        assert np.abs(np.mean(np.real(W))) < 1e-5
+
+    assert_allclose(np.linalg.norm(W), np.sqrt(2), atol=1e-6)
+
+    # Convert to SciPy nomenclature and compare
+    M = len(W)
+    w = n_cycles
+    s = w * sfreq / (2 * freq * np.pi)  # from SciPy docs
+    Ws = morlet2(M, s, w) * np.sqrt(2)
+    assert_allclose(W, Ws)
+
+    # Check FWHM
+    fwhm_formula = fwhm(freq, n_cycles)
+    half_max = np.abs(W).max() / 2.
+    fwhm_empirical = (np.abs(W) >= half_max).sum() / sfreq
+    # Could be off by a few samples
+    assert_allclose(fwhm_formula, fwhm_empirical, atol=3 / sfreq)
 
 
 def test_time_frequency():
@@ -101,7 +124,8 @@ def test_time_frequency():
                    freqs=freqs, n_cycles=n_cycles, use_fft=True,
                    return_itc=False, picks=picks, average=False)
     assert_allclose(
-        epochs_power_picks.data[0, 0, 0, 0], 9.130315e-23, rtol=1e-4)
+        epochs_power_picks.data[0, 0, 0, 0], 9.130315e-23,
+        rtol=1e-4)
     power_picks_avg = epochs_power_picks.average()
     # the actual data arrays here are equivalent, too...
     assert_allclose(power.data, power_picks.data)

mne/time_frequency/tfr.py

@@ -39,17 +39,19 @@
                          _setup_vmin_vmax, _set_title_multiple_electrodes)
 
 
+@fill_doc
 def morlet(sfreq, freqs, n_cycles=7.0, sigma=None, zero_mean=False):
     """Compute Morlet wavelets for the given frequency range.
 
     Parameters
     ----------
     sfreq : float
         The sampling Frequency.
-    freqs : array
-        Frequency range of interest (1 x Frequencies).
-    n_cycles : float | array of float, default 7.0
-        Number of cycles. Fixed number or one per frequency.
+    freqs : float | array-like, shape (n_freqs,)
+        Frequencies to compute Morlet wavelets for.
+    n_cycles : float | array-like, shape (n_freqs,)
+        Number of cycles. Can be a fixed number (float) or one per frequency
+        (array-like).
     sigma : float, default None
         It controls the width of the wavelet ie its temporal
         resolution. If sigma is None the temporal resolution
@@ -63,45 +65,143 @@ def morlet(sfreq, freqs, n_cycles=7.0, sigma=None, zero_mean=False):
 
     Returns
     -------
-    Ws : list of array
-        The wavelets time series.
-    """
+    Ws : list of ndarray | ndarray
+        The wavelets time series. If ``freqs`` was a float, a single
+        ndarray is returned instead of a list of ndarray.
+
+    See Also
+    --------
+    mne.time_frequency.fwhm
+
+    Notes
+    -----
+    %(morlet_notes)s
+    %(fwhm_morlet_notes)s
+
+    References
+    ----------
+    .. footbibliography::
+
+    Examples
+    --------
+    Let's show a simple example of the relationship between ``n_cycles`` and
+    the FWHM using :func:`mne.time_frequency.fwhm`, as well as the equivalent
+    call using :func:`scipy.signal.morlet2`:
+
+    .. plot::
+
+        import numpy as np
+        from scipy.signal import morlet2 as sp_morlet
+        import matplotlib.pyplot as plt
+        from mne.time_frequency import morlet, fwhm
+
+        sfreq, freq, n_cycles = 1000., 10, 7  # i.e., 700 ms
+        this_fwhm = fwhm(freq, n_cycles)
+        wavelet = morlet(sfreq=sfreq, freqs=freq, n_cycles=n_cycles)
+        M, w = len(wavelet), n_cycles # convert to SciPy convention
+        s = w * sfreq / (2 * freq * np.pi)  # from SciPy docs
+        wavelet_sp = sp_morlet(M, s, w) * np.sqrt(2)  # match our normalization
+
+        _, ax = plt.subplots(constrained_layout=True)
+        colors = {
+            ('MNE', 'real'): '#66CCEE',
+            ('SciPy', 'real'): '#4477AA',
+            ('MNE', 'imag'): '#EE6677',
+            ('SciPy', 'imag'): '#AA3377',
+        }
+        lw = dict(MNE=2, SciPy=4)
+        zorder = dict(MNE=5, SciPy=4)
+        t = np.arange(-M // 2 + 1, M // 2 + 1) / sfreq
+        for name, w in (('MNE', wavelet), ('SciPy', wavelet_sp)):
+            for kind in ('real', 'imag'):
+                ax.plot(t, getattr(w, kind), label=f'{name} {kind}',
+                        lw=lw[name], color=colors[(name, kind)],
+                        zorder=zorder[name])
+        ax.plot(t, np.abs(wavelet), label=f'MNE abs', color='k', lw=1., zorder=6)
+        half_max = np.max(np.abs(wavelet)) / 2.
+        ax.plot([-this_fwhm / 2., this_fwhm / 2.], [half_max, half_max],
+                color='k', linestyle='-', label='FWHM', zorder=6)
+        ax.legend(loc='upper right')
+        ax.set(xlabel='Time (s)', ylabel='Amplitude')
+    """  # noqa: E501
     Ws = list()
-    n_cycles = np.atleast_1d(n_cycles)
+    n_cycles = np.array(n_cycles, float).ravel()
 
-    freqs = np.array(freqs)
+    freqs = np.array(freqs, float)
     if np.any(freqs <= 0):
         raise ValueError("all frequencies in 'freqs' must be "
                          "greater than 0.")
 
     if (n_cycles.size != 1) and (n_cycles.size != len(freqs)):
         raise ValueError("n_cycles should be fixed or defined for "
                          "each frequency.")
+    _check_option('freqs.ndim', freqs.ndim, [0, 1])
+    singleton = freqs.ndim == 0
+    if singleton:
+        freqs = freqs[np.newaxis]
     for k, f in enumerate(freqs):
         if len(n_cycles) != 1:
             this_n_cycles = n_cycles[k]
         else:
             this_n_cycles = n_cycles[0]
-        # fixed or scale-dependent window
+        # sigma_t is the stddev of gaussian window in the time domain; can be
+        # scale-dependent or fixed across freqs
         if sigma is None:
             sigma_t = this_n_cycles / (2.0 * np.pi * f)
         else:
             sigma_t = this_n_cycles / (2.0 * np.pi * sigma)
-        # this scaling factor is proportional to (Tallon-Baudry 98):
-        # (sigma_t*sqrt(pi))^(-1/2);
+        # time vector. We go 5 standard deviations out to make sure we're
+        # *very* close to zero at the ends. We also make sure that there's a
+        # sample at exactly t=0
         t = np.arange(0., 5. * sigma_t, 1.0 / sfreq)
         t = np.r_[-t[::-1], t[1:]]
         oscillation = np.exp(2.0 * 1j * np.pi * f * t)
-        gaussian_enveloppe = np.exp(-t ** 2 / (2.0 * sigma_t ** 2))
-        if zero_mean:  # to make it zero mean
+        if zero_mean:
+            # this offset is equivalent to the κ_σ term in Wikipedia's
+            # equations, and satisfies the "admissibility criterion" for CWTs
             real_offset = np.exp(- 2 * (np.pi * f * sigma_t) ** 2)
             oscillation -= real_offset
-        W = oscillation * gaussian_enveloppe
+        gaussian_envelope = np.exp(-t ** 2 / (2.0 * sigma_t ** 2))
+        W = oscillation * gaussian_envelope
+        # the scaling factor here is proportional to what is used in
+        # Tallon-Baudry 1997: (sigma_t*sqrt(pi))^(-1/2).  It yields a wavelet
+        # with norm sqrt(2) for the full wavelet / norm 1 for the real part
         W /= np.sqrt(0.5) * np.linalg.norm(W.ravel())
         Ws.append(W)
+    if singleton:
+        Ws = Ws[0]
     return Ws
 
 
+def fwhm(freq, n_cycles):
+    """Compute the full-width half maximum of a Morlet wavelet.
+
+    Uses the formula from :footcite:t:`Cohen2019`.
+
+    Parameters
+    ----------
+    freq : float
+        The oscillation frequency of the wavelet.
+    n_cycles : float
+        The duration of the wavelet, expressed as the number of oscillation
+        cycles.
+
+    Returns
+    -------
+    fwhm : float
+        The full-width half maximum of the wavelet.
+
+    Notes
+    -----
+     .. versionadded:: 1.3
+
+    References
+    ----------
+    .. footbibliography::
+    """
+    return n_cycles * np.sqrt(2 * np.log(2)) / (np.pi * freq)
+
+
 def _make_dpss(sfreq, freqs, n_cycles=7., time_bandwidth=4.0, zero_mean=False):
     """Compute DPSS tapers for the given frequency range.
 
@@ -744,7 +844,16 @@ def tfr_morlet(inst, freqs, n_cycles, use_fft=False, return_itc=True, decim=1,
 
     Notes
     -----
+    %(morlet_notes)s
     %(temporal-window_tfr_notes)s
+    %(fwhm_morlet_notes)s
+
+    See :func:`mne.time_frequency.morlet` for more information about the
+    Morlet wavelet.
+
+    References
+    ----------
+    .. footbibliography::
     """
     tfr_params = dict(n_cycles=n_cycles, n_jobs=n_jobs, use_fft=use_fft,
                       zero_mean=zero_mean, output=output)
@@ -811,9 +920,14 @@ def tfr_array_morlet(epoch_data, sfreq, freqs, n_cycles=7.0,
 
     Notes
     -----
+    %(morlet_notes)s
     %(temporal-window_tfr_notes)s
 
     .. versionadded:: 0.14.0
+
+    References
+    ----------
+    .. footbibliography::
     """
     return _compute_tfr(epoch_data=epoch_data, freqs=freqs,
                         sfreq=sfreq, method='morlet', n_cycles=n_cycles,