nb-4-2-wavelets · STFT and Morlet wavelets (L4.2)¶
Lesson L4.2 · Level 4 · Status draft — for expert review; uncertain points carry TODO(confirm).
A Morlet wavelet is a windowed sinusoid, and the window's width fixes everything: how sharply the estimate can localise an event in time, how sharply it can separate two frequencies, and how far from the edges of an epoch the estimate is still made of data rather than of padding. Those three facts are one fact, and this notebook measures it rather than asserting it.
The last cell prints the L4.2 answer key: the temporal and spectral full width at half maximum of a 10 Hz wavelet with 7 cycles, with the convention named, because there are two conventions in use and they differ by √2.
Data. ds-eegbci — EEGMMIDB, Schalk et al. (2004), IEEE Trans Biomed Eng 51(6), DOI
10.1109/TBME.2004.827072; dataset DOI
10.13026/C28G6P. From data/directory.yaml: 64-channel 10-10 cap, 160 Hz,
no online filters, 60 Hz mains, access: open, licence ODC-By 1.0. Runs R04, R08, R12 (left-versus-right-fist
imagery), subject S001, channel C3.
Three EDF files, about 2.4 MB each, are downloaded and deleted in a finally at the end; free disk is printed
before and after.
In [1]:
# Setup: dependencies, the shared helpers, non-interactive plotting.
import importlib.util
import subprocess
import sys
import warnings
from pathlib import Path
# 1. Dependencies are pinned in notebooks/requirements.txt. Nothing is installed when the
# pinned stack is already present (local runs, CI); a fresh Colab or Binder kernel installs
# it once. On Colab, run from a clone of the repository so that notebooks/_shared/ is
# available (repository URL: TODO(confirm), spec section 13 item 3).
_needed = ('mne', 'scipy', 'matplotlib', 'pooch')
_missing = [p for p in _needed if importlib.util.find_spec(p) is None]
if _missing:
_req = next((d / "requirements.txt" for d in (Path.cwd(), *Path.cwd().parents)
if (d / "requirements.txt").exists()), None)
_cmd = [sys.executable, "-m", "pip", "install", "-q"]
_cmd += ["-r", str(_req)] if _req else ["mne==1.10.2", "pooch>=1.8"]
subprocess.check_call(_cmd)
# 2. Shared helpers, located relative to the working directory -- notebooks/<level>/ or
# notebooks/ -- never through an absolute path.
_shared = next((d / "_shared" for d in (Path.cwd(), *Path.cwd().parents)
if (d / "_shared" / "helpers_l4.py").exists()), None)
if _shared is None:
raise FileNotFoundError("start the kernel in notebooks/L4/ (or notebooks/) so that "
"_shared/helpers_l4.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l1
import helpers_l4 as L4
# 3. Plotting: Jupyter's default inline backend renders static PNGs through Agg (no windows,
# nothing blocks); outside Jupyter the helpers select Agg. Every MNE figure is requested
# with show=False and each figure cell ends with plt.show().
import matplotlib.pyplot as plt
import numpy as np
import mne
mne.set_log_level("WARNING")
plt.rcParams["figure.dpi"] = 72
# 4. Quiet the downloader. pooch, which MNE uses to fetch datasets, logs
# "Downloading file '...' from '...' to '<cache directory>'" at INFO, and that last field is an
# ABSOLUTE PATH from whichever machine executed the notebook. Absolute paths are not allowed in a
# stored notebook (scripts/scrub-notebooks.py is a CI gate) and re-executing would put them straight
# back, so the message is suppressed at the source rather than cleaned up afterwards. Nothing is
# hidden by this: every cell below prints the file NAMES it fetched and helpers_l4.Downloads prints
# free disk before and after. Please do not delete this as noise.
try:
import pooch
pooch.get_logger().setLevel("WARNING")
except Exception: # pooch absent or its API moved: the scrub script is the backstop
pass
print(f"MNE {mne.__version__}; helpers_l4 imported from notebooks/_shared")
print(f"downloads go to {helpers_l1.download_dir().name}/ (resolved relative to the working directory, "
"or $EEG_COURSE_DOWNLOADS) and are deleted at the end of this notebook")
1 · The wavelet, measured rather than described¶
MNE builds the wavelet from two numbers: the frequency f and the number of cycles n. The Gaussian envelope has
sigma_t = n / (2 pi f) seconds
sigma_f = 1 / (2 pi sigma_t) = f / n hertzand the full width at half maximum of a Gaussian of standard deviation σ is 2 sqrt(2 ln 2) σ.
Which half maximum, though? The half maximum of the amplitude envelope, or of the power (the squared
envelope)? The second is narrower by √2, and both appear in the literature. This notebook uses the amplitude
convention by default, for two reasons: it is what mne.time_frequency.fwhm computes, and it is the width the
kernel MNE actually convolves with has. The other is printed beside it, because a pair of numbers without the
convention has two right answers.
The cell below does not trust the algebra. It builds the kernel with mne.time_frequency.morlet, finds the half
maximum of |W(t)| by interpolation, does the same on |FFT(W)|, and compares all three routes.
In [2]:
FREQ, NCYC = 10.0, 7.0
an = L4.wavelet_fwhm(FREQ, NCYC, "amplitude")
pw = L4.wavelet_fwhm(FREQ, NCYC, "power")
lib = float(mne.time_frequency.fwhm(FREQ, NCYC))
print(f"A {FREQ:g} Hz Morlet wavelet with {NCYC:g} cycles")
print(f" sigma_t = n / (2 pi f) = {an['sigma_t_s'] * 1000:.6f} ms")
print(f" sigma_f = f / n = {an['sigma_f_hz']:.6f} Hz")
print()
print(f" {'convention':10s} {'FWHM in time':>16s} {'FWHM in frequency':>19s} {'product':>12s} identity")
for d in (an, pw):
print(f" {d['convention']:10s} {d['fwhm_t_ms']:13.5f} ms {d['fwhm_f_hz']:16.6f} Hz "
f"{d['product_s_hz']:12.6f} {d['identity']} = {d['identity_value']:.6f} s.Hz")
print()
print(f" mne.time_frequency.fwhm({FREQ:g}, {NCYC:g}) = {lib * 1000:.6f} ms "
f"-- {'the amplitude convention, to the last printed digit' if abs(lib - an['fwhm_t_s']) < 1e-12 else 'DIFFERENT from both conventions above'}")
for sf in (500.0, 1000.0, 2000.0):
m = L4.measured_wavelet_fwhm(FREQ, NCYC, sfreq=sf)
print(f" measured off MNE's own kernel at {sf:6.0f} Hz: {m['fwhm_t_ms']:10.4f} ms in time, "
f"{m['fwhm_f_hz']:.5f} Hz in frequency ({m['kernel_samples']} samples, "
f"{m['kernel_duration_s']:.3f} s long)")
print()
print("The measured widths converge on the analytic amplitude-convention values as the kernel is sampled more "
"finely, which is the check that the algebra describes the object the library convolves with.")
In [3]:
fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.0))
L4.plot_wavelet(FREQ, NCYC, sfreq=1000.0, axes=axes)
fig.suptitle(f"The kernel behind every map in this notebook: {FREQ:g} Hz, {NCYC:g} cycles "
f"(amplitude convention; FWHM {an['fwhm_t_ms']:.3f} ms and {an['fwhm_f_hz']:.5f} Hz)", y=1.03)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
print("The uncertainty relation, stated as an identity rather than as an inequality:")
print(f" FWHM_t * FWHM_f = (2 sqrt(2 ln2) n / (2 pi f)) * (2 sqrt(2 ln2) f / n) = 8 ln2 / (2 pi) "
f"= 4 ln2 / pi = {4 * np.log(2) / np.pi:.6f} s.Hz, for every f and every n.")
print(" On the power convention both widths shrink by sqrt(2), so the product is half that: "
f"2 ln2 / pi = {2 * np.log(2) / np.pi:.6f} s.Hz.")
print()
print(" TODO(confirm) for the data track: site/public/data/widgets/w-tf-baseline-explorer/tfpower.json "
"states this identity as '2 ln 2 / pi' while its own fwhm_t_s and fwhm_f_hz arrays multiply to "
f"{4 * np.log(2) / np.pi:.6f} = 4 ln2 / pi. The arrays are right and the sentence beside them is "
"off by a factor of two; nothing computes from the sentence, so no number moves, but it should be "
"corrected.")
print(" Checked here on that file's own numbers for the 4 Hz row: "
f"fwhm_t 0.2811 s x fwhm_f 3.1398 Hz = {0.2811 * 3.1398:.4f} s.Hz.")
2 · Cycles on real epochs: the trade-off is visible, not theoretical¶
Four maps of the same 22 trials, differing only in the cycle scheme:
- 3 cycles everywhere — sharp in time, blurred in frequency.
- 7 cycles everywhere — the reverse, and at 4 Hz the kernel is 1.4 s long, so it eats the epoch's edges.
- 12 cycles everywhere — frequency-sharp to the point of being useless for anything transient.
max(3, f/2)— the variable scheme this course uses: roughly constant time resolution above 6 Hz, with a floor below it so that the low-frequency kernels do not grow without limit.
The shaded band at each end is where the estimate is partly made of the zero padding MNE adds. Its width is
5 σ_t, so it is a function of the cycle scheme and of nothing else — which is why it changes between panels.
In [4]:
SUBJECT, CH, COND = "S001", "C3", "T2"
dl = L4.Downloads("nb-4-2").start()
epochs, info = L4.load_imagery_epochs(SUBJECT, L4.MI_RUNS, dl)
x = epochs[COND].get_data(picks=[CH])[:, 0, :] * 1e6
times_raw, sf = epochs.times, float(epochs.info["sfreq"])
print(f"{SUBJECT} {COND} ({L4.MI_CONDITIONS[COND]}): {x.shape[0]} trials, channel {CH}, "
f"{times_raw[0]:g}..{times_raw[-1]:g} s at {sf:g} Hz")
print("downloaded for this notebook: " + ", ".join(p.name for p in dl.paths))
In [5]:
SCHEMES = {
"3 cycles (fixed)": np.full(len(L4.TF_FREQS), 3.0),
"7 cycles (fixed)": np.full(len(L4.TF_FREQS), 7.0),
"12 cycles (fixed)": np.full(len(L4.TF_FREQS), 12.0),
"max(3, f/2) (this course)": L4.TF_CYCLES,
}
fig, axes = plt.subplots(1, 4, figsize=(21, 4.0))
scheme_stats = []
for ax, (name, ncyc) in zip(axes, SCHEMES.items()):
r = L4.morlet_power(x[:, None, :], sfreq=sf, freqs=L4.TF_FREQS, n_cycles=ncyc, decim=L4.TF_DECIM,
times=times_raw)
P = L4.baseline_normalise(r["power"].mean(0)[0], r["times"], L4.TF_BASELINE, mode="db")
e = L4.edge_seconds(L4.TF_FREQS, ncyc)
L4.plot_tfr(P, L4.TF_FREQS, r["times"], ax=ax, vlim=(-6, 6), edge_s=e, baseline=L4.TF_BASELINE,
cbar_label="Power change from baseline (dB)",
title=f"{name}\n{SUBJECT} {CH} {COND}, dB re {L4.TF_BASELINE[0]:g}..{L4.TF_BASELINE[1]:g} s")
at10 = L4.wavelet_fwhm(10.0, float(np.interp(10.0, L4.TF_FREQS, ncyc)))
scheme_stats.append({"scheme": name,
"n_cycles at 10 Hz": float(np.interp(10.0, L4.TF_FREQS, ncyc)),
"FWHM_t at 10 Hz (ms)": at10["fwhm_t_ms"],
"FWHM_f at 10 Hz (Hz)": at10["fwhm_f_hz"],
"edge at 4 Hz (s)": float(e[0]),
"edge at 40 Hz (s)": float(e[-1])})
fig.suptitle("The same 22 trials under four cycle schemes (dB re baseline; shaded = partly zero padding)", y=1.04)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
print(L4.fmt_table(scheme_stats, floatfmt="{:.3f}"))
print()
print("Reading the table: the time and frequency columns move in opposite directions and their product is "
f"{4 * np.log(2) / np.pi:.6f} s.Hz in every row. Choosing cycles is choosing where on that curve to sit; "
"there is no setting that is sharp in both.")
3 · Edge effects, and why padding does not fix them¶
MNE zero-pads before convolving, so an estimate near the edge of an epoch is an average of data and of zeros. It does not error, it does not warn, and it looks like a real result — usually like a drop in power, because the zeros carry none.
The demonstration below is a controlled one. A synthetic burst of constant amplitude runs the whole length of an
epoch, so the true map is flat. Anything the estimate does at the edges is an artefact of the estimator, and
the size of the region it happens in is exactly 5 σ_t.
In [6]:
SF_SYN, DUR = 160.0, 6.0
t_syn = np.arange(int(DUR * SF_SYN)) / SF_SYN - DUR / 2
sine = 10.0 * np.cos(2 * np.pi * 10.0 * t_syn) # 10 uV, 10 Hz, constant for the whole epoch
fig, axes = plt.subplots(1, 3, figsize=(17, 3.9))
for ax, ncyc_val in zip(axes, (3.0, 7.0, 12.0)):
r = L4.morlet_power(sine[None, None, :], sfreq=SF_SYN, freqs=np.arange(4.0, 20.5, 1.0),
n_cycles=ncyc_val, decim=1, times=t_syn)
row = r["power"][0, 0, np.argmin(np.abs(r["freqs"] - 10.0))]
ed = 5 * ncyc_val / (2 * np.pi * 10.0)
inner = np.abs(t_syn) <= (DUR / 2 - ed)
ax.plot(t_syn, 100 * row / np.median(row[inner]), lw=1.4, color="k")
ax.axvspan(t_syn[0], t_syn[0] + ed, color="0.4", alpha=0.28, lw=0)
ax.axvspan(t_syn[-1] - ed, t_syn[-1], color="0.4", alpha=0.28, lw=0)
ax.axhline(100, color="tab:blue", lw=1.0, ls="--")
ax.set_ylim(0, 130)
ax.set(xlabel="Time (s)", ylabel="Estimated 10 Hz power (% of the interior median)",
title=f"SYNTHETIC constant 10 Hz, 10 uV: {ncyc_val:g} cycles\n5 sigma_t = {ed:.3f} s shaded; "
f"power at the first sample = {100 * row[0] / np.median(row[inner]):.0f} %")
ax.grid(alpha=0.3)
fig.suptitle("A flat signal does not give a flat map at the edges (synthetic; the true answer is 100 % "
"everywhere)", y=1.04)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
print("Two ways of dealing with it, and only one of them is honest:")
print(" 1. Cut the epoch longer than you need and analyse only the interior. The padding then falls outside "
"the window you report. This is what MI_PIPELINE does: the epoch runs -2.5..4.5 s and the baseline "
f"({L4.TF_BASELINE[0]:g}..{L4.TF_BASELINE[1]:g} s) starts "
f"{L4.TF_BASELINE[0] - (-2.5 + L4.edge_seconds()[0]):.3f} s after the edge region ends.")
print(" 2. Pad with more zeros. This moves nothing: the estimate near the edge is still an average of data "
"and of not-data. Padding changes the interpolation of the spectrum, not the resolution, and it cannot "
"create samples that were never recorded (pf-tf-edge-effects).")
4 · The STFT, and what "variable resolution" buys¶
The short-time Fourier transform uses one window length for every frequency; a wavelet transform uses a window of a fixed number of cycles, so the window shortens as the frequency rises. Below, the same epoch through two STFTs (a long window and a short one) and through the course's wavelet scheme.
A 0.5 s STFT window gives a 2 Hz frequency resolution at every frequency — generous at 40 Hz, too coarse at 5 Hz. A 2 s window gives 0.5 Hz everywhere — good at 5 Hz, and hopeless for anything that happens faster than 2 s.
In [7]:
from scipy import signal as sps
avg_trials = x # the same 22 trials
fig, axes = plt.subplots(1, 3, figsize=(17.5, 4.1))
for ax, win_s in zip(axes[:2], (0.5, 2.0)):
nper = int(round(win_s * sf))
f_st, t_st, Z = sps.stft(avg_trials, fs=sf, window="hann", nperseg=nper,
noverlap=nper - max(1, nper // 16), boundary=None, padded=False)
P = (np.abs(Z) ** 2).mean(0) # average of single-trial power
t_abs = t_st + times_raw[0]
keep = (f_st >= 4) & (f_st <= 40)
Pn = L4.baseline_normalise(P[keep], t_abs, L4.TF_BASELINE, mode="db")
L4.plot_tfr(Pn, f_st[keep], t_abs, ax=ax, vlim=(-6, 6), baseline=L4.TF_BASELINE,
edge_s=win_s / 2, cbar_label="Power change from baseline (dB)",
title=f"STFT, {win_s:g} s Hann window\nresolution {1 / win_s:.2f} Hz at EVERY frequency; "
f"{len(t_st)} time bins")
r = L4.morlet_power(x[:, None, :], sfreq=sf, times=times_raw)
Pm = L4.baseline_normalise(r["power"].mean(0)[0], r["times"], L4.TF_BASELINE, mode="db")
L4.plot_tfr(Pm, r["freqs"], r["times"], ax=axes[2], vlim=(-6, 6), baseline=L4.TF_BASELINE,
edge_s=L4.edge_seconds(), cbar_label="Power change from baseline (dB)",
title="Morlet, max(3, f/2) cycles\nresolution varies with frequency")
fig.suptitle(f"Same trials, three estimators ({SUBJECT} {CH} {COND}, dB re "
f"{L4.TF_BASELINE[0]:g}..{L4.TF_BASELINE[1]:g} s)", y=1.04)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
print(f"{'estimator':28s} {'resolution at 5 Hz':>20s} {'at 10 Hz':>12s} {'at 40 Hz':>12s}")
for win_s in (0.5, 2.0):
print(f"{'STFT, ' + format(win_s, 'g') + ' s Hann':28s} "
+ "".join(f"{1 / win_s:12.3f} Hz" for _ in (5, 10, 40)))
row = "".join(f"{L4.wavelet_fwhm(f0, max(3.0, f0 / 2))['fwhm_f_hz'] / (2 * np.sqrt(2 * np.log(2))):12.3f} Hz"
for f0 in (5.0, 10.0, 40.0))
print(f"{'Morlet, max(3, f/2)':28s} {row}")
print("(quoted as sigma_f so that the three are on one definition; multiply by "
f"{2 * np.sqrt(2 * np.log(2)):.4f} for the FWHM)")
5 · The numbers¶
In [8]:
try:
WIDGET_FWHM_T_MS, WIDGET_FWHM_F_HZ = 262.347, 3.36403 # site/notes/integration-phase3.md, widgets-J
print("nb-4-2-wavelets -- L4.2 exercise numbers (draft; TODO(confirm) at author review)")
print(f"Data: ds-eegbci {SUBJECT}, runs {'+'.join(info['runs'])} "
f"(PhysioNet DOI {L4.DATASETS_L4['ds-eegbci']['dataset_doi']}; ODC-By 1.0); "
f"{info['n_eeg']} EEG channels at {info['sfreq']:g} Hz, average reference, no filtering. "
f"Condition {COND}, {x.shape[0]} trials, channel {CH}.")
print(f"Wavelets: mne.time_frequency.morlet / tfr_array_morlet, MNE {mne.__version__}.")
print()
print(f"ANSWER KEY -- ex-4-2 (numeric pair): temporal and spectral FWHM of a {FREQ:g} Hz wavelet with "
f"{NCYC:g} cycles")
print(f" {an['fwhm_t_ms']:.3f} ms and {an['fwhm_f_hz']:.5f} Hz "
f"-- AMPLITUDE-ENVELOPE convention (sigma_t = n / (2 pi f), FWHM = 2 sqrt(2 ln2) sigma)")
print(f" THE PROMPT MUST NAME THE CONVENTION. On the POWER convention the same wavelet gives "
f"{pw['fwhm_t_ms']:.3f} ms and {pw['fwhm_f_hz']:.5f} Hz, narrower by sqrt(2).")
print(f" The amplitude convention is the one mne.time_frequency.fwhm returns "
f"({lib * 1000:.5f} ms) and the one the kernel measurably has.")
print(f" Neighbouring cycle counts, so a tolerance can be checked against them:")
for nc in (6.0, 7.0, 8.0):
d = L4.wavelet_fwhm(FREQ, nc)
print(f" {nc:g} cycles: {d['fwhm_t_ms']:8.3f} ms, {d['fwhm_f_hz']:.5f} Hz")
print(f" A tolerance of +/- 2 ms and +/- 0.03 Hz excludes the power convention and both neighbours.")
print(f" Cross-check: site/notes/integration-phase3.md records {WIDGET_FWHM_T_MS:g} ms and "
f"{WIDGET_FWHM_F_HZ:g} Hz from widgets-J. This notebook "
f"{'AGREES to the digits quoted' if (abs(an['fwhm_t_ms'] - WIDGET_FWHM_T_MS) < 5e-3 and abs(an['fwhm_f_hz'] - WIDGET_FWHM_F_HZ) < 5e-6) else 'DISAGREES -- report both'}.")
print()
print("Supporting -- the uncertainty product, which is the reason the pair is a pair:")
print(f" FWHM_t * FWHM_f = {an['product_s_hz']:.6f} s.Hz = 4 ln2 / pi, for every frequency and every "
f"cycle count (half that on the power convention).")
print(f" Measured off MNE's kernel at 1000 Hz: {L4.measured_wavelet_fwhm(FREQ, NCYC, 1000.0)['fwhm_t_ms']:.4f} ms "
f"and {L4.measured_wavelet_fwhm(FREQ, NCYC, 1000.0)['fwhm_f_hz']:.5f} Hz.")
print()
print("Supporting -- the edge region, which the cycle scheme fixes and nothing else does:")
for nm, nc in (("3 cycles", 3.0), ("7 cycles", 7.0), ("max(3, f/2)", None)):
e = L4.edge_seconds(L4.TF_FREQS, L4.TF_CYCLES if nc is None else np.full(len(L4.TF_FREQS), nc))
print(f" {nm:12s}: {e[0]:.4f} s at 4 Hz, {e[np.argmin(abs(L4.TF_FREQS - 10))]:.4f} s at 10 Hz, "
f"{e[-1]:.4f} s at 40 Hz")
print()
print("Pitfall: pf-tf-edge-effects. Widget: w-wavelet-explorer (mode morlet).")
print("TODO(confirm) for the data track: w-tf-baseline-explorer/tfpower.json states the uncertainty "
"identity as '2 ln 2 / pi' where its own arrays give 4 ln 2 / pi. The arrays are right.")
finally:
dl.finish()