MNE-Python's implementation of Morlet

đź’¬ Support & Discussions
1

I’m trying to wrap my head around MNE-Python’s implementation of Morlet wavelets, especially as compared to SciPy’s implementation. Here’s what I think I’ve figured out so far:

scipy’s morlet2 function parameterizes in terms of M (total number of samples in the wavelet), s (“width”) and w (“omega0”). MNE parameterizes in terms of freq (frequency in Hz), n_cycles (float), and sigma (float or None). For now I’m only considering the sigma=None condition (in which n_cycles is fixed across frequencies and temporal resolution varies with frequency). I’m pretty sure that SciPy’s w AKA “omega0” is the same as n_cycles, so I think that means that:

w = n_cycles
s = sfreq * w / (2 * pi * freq)
M = 2 * ceil(s) - 1

However, if I define the SciPy inputs in terms of the MNE-Python inputs in this way, The results don’t match: I have to add in a factor of 5 into the definition of M, and then the SciPy output matches the MNE output (up to a scaling factor):

M = 2 * ceil(5 * s) - 1  # note the 5 here

This multiplier 5 shows up in the MNE-Python code here, in the definition of the time vector: mne-python/tfr.py at f26528d78764c83f754873c40f17e40d5eb08d2d · mne-tools/mne-python · GitHub, but I can’t figure out why it is there. I know it’s been over a decade since you wrote this code @agramfort, but any chance you can explain how that 5 got there? Based on the comments you seem to be following the math in https://doi.org/10.1523/JNEUROSCI.18-11-04244.1998 but I can’t find a 5 anywhere in there either.

MWE:

import matplotlib.pyplot as plt
import numpy as np
from mne.time_frequency import morlet
from scipy.signal import morlet2

sfreq = 1000
freqs = [17]
n_cycles = 3

w = n_cycles
s = sfreq * w / (2 * np.pi * freqs[0])
M_1 = int(2 * np.ceil(s) - 1)
M_5 = int(2 * np.ceil(5 * s) - 1)  # note the 5 here

mne_wav = morlet(sfreq=sfreq, freqs=freqs, n_cycles=n_cycles, zero_mean=False)[0]
scipy_wav_1 = morlet2(M=M_1, s=s, w=w)
scipy_wav_5 = morlet2(M=M_5, s=s, w=w)

with plt.ion():
    fig, axs = plt.subplots(2, 1, sharex=True)
    axs[0].plot(mne_wav.real)
    axs[1].plot(scipy_wav_5.real, label='with five')
    axs[1].plot(scipy_wav_1.real, label='without five')
    axs[0].set_title('MNE')
    axs[1].set_title('SciPy')
    axs[1].legend()

Figure_4

Another question I have is why the energy normalization is different between SciPy and MNE; SciPy multiplies by pi ** (-1/4) * sqrt(1 / s) whereas MNE divides by sqrt(0.5) * np.linalg.norm(W.ravel()). Wikipedia seems to agree more closely with SciPy on this one (at least, there is a pi ** -1/4 term in their equation: Morlet wavelet - Wikipedia).

1 Like
2

discussed with @agramfort during office hours. the factor of 5 is related to the width at which the gaussian envelope gets sufficiently close to zero. The second question about extra normalization in MNE is probably related to the fact that we allow users to pass a fixed sigma for all frequencies if they want, which might lead to the wavelets not having total energy of 1. I’ll do some experiments and report back if that all turns out to be accurate.

2 Likes
3

I found this thread because I am trying to port an analysis workflow from matlab to MNE-python and was getting discrepancies for my wavelet analysis. It seems that I get nearly equivalent results with matlab and scipy, but the MNE results differ, presumably due to the difference you identified here. Have you figured out the reasoning for the MNE implementation? I’m curious as to which method is preferable…

Thanks!

4

yes, I’ve continued looking into this. Here’s a figure that illustrates part of the issue:

Figure_1
Figure_11533Ă—1067 250 KB

As you can see the gaussian envelope is narrower for MNE than for SciPy. The MATLAB envelope is intermediate between the two, at least when using the default fb=1 (caveat: I implemented the MATLAB version in Python, based on the equations published in their docs. I didn’t actually use the MATLAB functions directly, because I only have access to Octave which doesn’t have the wavelets toolbox.)

A second issue relates to the y-axis scales on each subplot. There seems to be quite some variety in how these things are normalized. MNE guarantees unit L2 norm for each wavelet (i.e., total energy of 1). The two SciPy functions differ only in their normalization. We didn’t have time to discuss that aspect of it in today’s dev meeting, so I don’t really have an update on that front, other than to say that MNE will probably keep normalizing the same way, but based on today’s discussion it’s likely that we’ll change our envelope. Once I have a PR open I’ll link to it here, so you can keep abreast of the discussion there.

Here’s the code to generate that figure:

import numpy as np
import matplotlib.pyplot as plt


def morlet_matlab(t, f, bandwidth=1):
    # https://www.mathworks.com/help/wavelet/ug/wavelet-families-additional-discussion.html#f8-40178
    # https://www.mathworks.com/help/wavelet/ref/cmorwavf.html
    envelope = np.exp(-t ** 2 / bandwidth)
    wave = np.exp(2 * np.pi * 1j * f * t)
    normalizer = 1 / np.sqrt(np.pi * bandwidth)
    return normalizer * envelope * wave


def morlet_wikipedia(t, f):
    # https://en.wikipedia.org/wiki/Morlet_wavelet#Definition
    # cites doi:10.12743/quanta.v1i1.5 (eqns 14 & 26 are the relevant ones)
    f *= 2 * np.pi  # this step isn't mentioned in the article; I inferred its necessity
    envelope = np.exp(-0.5 * t ** 2)
    wave = np.exp(1j * f * t) - np.exp(-0.5 * f ** 2)
    normalizer = (
        (np.pi ** -0.25) *
        (1 + np.exp(-1 * f ** 2) - 2 * np.exp(-0.75 * f ** 2)) ** -0.5
    )
    return normalizer * envelope * wave


def morlet_scipy(t, f, sfreq):
    from scipy.signal import morlet
    M = t.size
    s = t[-1] / (2 * np.pi)
    w = f * M / (2 * s * sfreq)
    return morlet(M=M, w=w, s=s, complete=True)


def morlet2_scipy(t, f, sfreq):
    from scipy.signal import morlet2
    s = sfreq
    w = 2 * np.pi * f * s / sfreq
    return morlet2(M=len(t), s=s, w=w)


def morlet_mne(t, f, sfreq):
    from mne.time_frequency import morlet
    sigma_t = t[-1] / 5
    n_cycles = sigma_t * 2 * np.pi * f
    return morlet(sfreq=sfreq, freqs=[f], n_cycles=n_cycles)[0]


sfreq = 1000
duration = 6
tmin, tmax = np.array([-0.5, 0.5]) * duration
n_samp = 1 + np.rint(duration * sfreq).astype(int)
time = np.linspace(tmin, tmax, num=n_samp, endpoint=True)
freq = 11

methods = dict(MATLAB=morlet_matlab,
               MNE=morlet_mne,
               SciPy1=morlet_scipy,
               SciPy2=morlet2_scipy,
               Wikipedia=morlet_wikipedia)

fig, axs = plt.subplots(len(methods), 1, sharex=True,
                        constrained_layout=True)

for ax, (method, func) in zip(axs, methods.items()):
    kwargs = dict()
    if method.startswith('SciPy') or method.startswith('MNE'):
        kwargs = dict(sfreq=sfreq)
    ax.plot(time, func(time, freq, **kwargs).real, label=method)
    ax.set_title(method)
2 Likes
5

cross-ref to GitHub PR: https://github.com/mne-tools/mne-python/pull/11353

1 Like