Morlet wavelets

The Morlet (Gabor) wavelet transform is the central, load-bearing feature of PTSA. Most lab analyses use PTSA primarily to obtain time-resolved power and phase estimates of an EEG signal in a set of frequency bands. This page is the canonical reference for the formula PTSA implements, the meaning of the user-facing parameters, and the differences between PTSA’s parameterization and the ones used by SciPy and MNE.

Formula

PTSA’s Morlet wavelet is the standard complex Morlet (Gabor) wavelet parameterized by the number of cycles under the Gaussian envelope. For a target frequency \(f\) and a width \(w\) (number of cycles):

\[\sigma_t = \frac{w}{2 \pi f}, \qquad \omega = 2 \pi f,\]

\[\psi(t) = \left(\sigma_t \sqrt{\pi}\right)^{-1/2} \, e^{-t^2 / (2 \sigma_t^2)} \, e^{i \omega t}.\]

The prefactor \((\sigma_t \sqrt{\pi})^{-1/2}\) gives the wavelet unit \(L^2\) energy. The Gaussian envelope \(e^{-t^2/(2 \sigma_t^2)}\) is evaluated out to \(\pm 3.5 \sigma_t\); values beyond that window are truncated to zero.

Time-frequency decomposition is carried out by convolving the input signal with \(\psi\). PTSA performs the convolution in the frequency domain using FFTW:

\[(\Psi * \psi)(t) = \mathcal{F}^{-1}\!\left[ \widehat{\Psi}(\xi) \cdot \widehat{\psi}(\xi) \right](t),\]

where \(\widehat{\,\cdot\,}\) denotes the discrete Fourier transform. The wavelet FFT is precomputed once per frequency and reused across signals and trials, which is the source of PTSA’s speed for long recording sessions.

The complete correction (default)

When complete=True (the default since PTSA 2.0.6), the wavelet is zero-mean corrected so that the integral \(\int \psi(t) \, dt\) vanishes exactly. The corrected wavelet is:

\[\tilde{\psi}(t) = \left(\sigma_t \sqrt{\pi}\right)^{-1/2} \, e^{-t^2 / (2 \sigma_t^2)} \, \left( e^{i \omega t} - e^{-w^2 / 2} \right).\]

The \(e^{-w^2/2}\) subtraction would slightly shift the peak frequency away from \(f\) if applied naively, so PTSA also rescales the time axis by

\[\alpha = \frac{2}{\pi} \arccos\!\left( e^{-w^2 / 2} \right),\]

so that the corrected wavelet’s peak in the frequency domain stays at \(f\).

When to flip complete to False:

• The correction matters most at small widths (e.g.width < 6), where the uncorrected wavelet has a non-trivial DC component that biases low-frequency power estimates.

• For larger widths the correction is numerically negligible (\(e^{-w^2/2}\) shrinks very fast: w=5 already gives \(\sim 4 \times 10^{-6}\)).

• Setting complete=False is mainly useful for reproducing analyses that explicitly used the uncorrected wavelet, or for cross-checking against another tool that does not apply a zero-mean correction.

Parameters

freqsarray of float

Frequencies (Hz) at which to evaluate the transform. The output has a frequency dimension of the same length. There is no requirement that the frequencies be log-spaced; PTSA simply builds one cached wavelet FFT per entry.

widthint, default 5

Number of cycles of the carrier sinusoid under one standard deviation of the Gaussian envelope (the same convention as MNE’s n_cycles). This is not a SciPy-style scale factor. The trade-off is the usual time-frequency uncertainty:

• Higher width → narrower bandwidth (tighter frequency resolution), wider time extent (looser time resolution).

• Lower width → broader bandwidth, sharper temporal localization.

Common values in the EEG literature are width=5 (the PTSA default) through width=7; width=4 is at the low end of what is typically reported, and the complete=True correction becomes important there.

outputstr or sequence of str, default ('power', 'phase')

Any combination of 'power' (squared magnitude of the complex coefficient), 'phase' (angle of the coefficient, in radians, in \((-\pi, \pi]\)), and 'complex' (the complex coefficient itself). 'complex' may not be combined with the other two on the same call.

completebool, default True

Apply the zero-mean correction described above.

cpusint, default 1

Number of worker threads used by the C++ kernel.

verbosebool, default True

Print kernel timing information.

Worked example

The example below runs end-to-end on synthetic data (no lab files required). It computes the time-resolved power of a 10 Hz sinusoid embedded in white noise and checks that the dominant frequency estimate is at 10 Hz.

import numpy as np
from ptsa.data.timeseries import TimeSeries
from ptsa.data.filters import MorletWaveletFilter

# 2 s of data at 500 Hz: 10 Hz sine plus 0.1-amplitude white noise.
rng = np.random.default_rng(0)
samplerate = 500.0
t = np.arange(0.0, 2.0, 1.0 / samplerate)
signal = np.sin(2 * np.pi * 10.0 * t) + 0.1 * rng.standard_normal(t.size)

# Wrap as a TimeSeries (samplerate coord is required).
ts = TimeSeries.create(
    signal[np.newaxis, :],
    samplerate,
    dims=('channel', 'time'),
    coords={'channel': [0], 'time': t},
)

# Decompose at 5, 10, 20, 40 Hz.
freqs = np.array([5.0, 10.0, 20.0, 40.0])
wf = MorletWaveletFilter(freqs=freqs, output='power',
                         width=5, complete=True,
                         cpus=1, verbose=False)
power = wf.filter(ts)

# power has shape (freq, channel, time).
assert power.shape == (4, 1, t.size)

# Mean power over the central 1 s (drop edge artifacts).
middle = (t > 0.5) & (t < 1.5)
mean_power = power.values[:, 0, middle].mean(axis=1)
dominant_freq = freqs[mean_power.argmax()]
assert dominant_freq == 10.0

The returned power is a TimeSeries with dims ('frequency', 'channel', 'time'). The leading boundary samples (roughly \(3.5 \sigma_t\) on each side, i.e.a few wavelet widths) carry edge artifacts; users typically discard them or add a mirror buffer with add_mirror_buffer() before filtering.

Parameterization comparison

Different libraries parameterize the Morlet wavelet differently. Mixing them is a common source of analysis bugs. The table below summarizes the three conventions you are most likely to encounter.

Library	Width parameter	Notes
PTSA (this package)	width = number of cycles under the Gaussian envelope (one \(\sigma_t\)).	complete=True by default applies the zero-mean correction.
mne.time_frequency.morlet	n_cycles = number of cycles, same convention as PTSA’s width.	Compatible parameterization. MNE also offers a zero-mean option (zero_mean=True) analogous to PTSA’s complete.
scipy.signal.morlet (removed in SciPy 1.15)	w = “omega-zero”, related to but not equal to the cycle count. SciPy used \(\sigma = 1\) in the wavelet and scaled the time axis via an additional s parameter.	The naming conflict (both w for SciPy and width for PTSA can be called “the wavelet width”) routinely confuses new users. SciPy values are not drop-in replacements for PTSA values.

Warning

Do not pass a SciPy-style w value to PTSA’s width (or vice versa). The two parameterizations differ by both definition and scale, and using the wrong one silently produces a well-formed but incorrect time-frequency map.