Multitaper and filter-Hilbert

Why this matters

Morlet is not the only way to get a time-frequency map, and the alternatives are not interchangeable. Multitaper buys a lower-variance estimate over a bandwidth you state explicitly, which is what a broad high-frequency band needs. Filter-Hilbert hands you the instantaneous amplitude and phase of one band, which is what phase analyses need — and it is also the easiest way in this whole level to manufacture a rhythm that is not there. Run all three on the same epoch and they disagree, and reading the disagreement correctly is the point of the lesson.

Concepts

The problem tapering solves

Cutting a segment out of a recording multiplies it by a rectangular window, and a rectangle has a badly behaved spectrum: a narrow main lobe with side lobes that decay slowly, so power at one frequency leaks across the whole axis (L1.4). Tapering the segment — Hann, Hamming, Gaussian — suppresses the side lobes at the cost of a wider main lobe. Every window is a different point on that trade-off, and every one of them throws away the samples at the edges of the segment, which is a real loss of data.

Multitaper takes a different route. Instead of one taper, use several that are mutually orthogonal, transform the segment with each, and average the resulting power estimates. Because the tapers are orthogonal, the estimates are approximately independent, so averaging K of them cuts the variance by roughly a factor of K; and because different tapers weight different parts of the segment, the data near the edges is used rather than discarded.

DPSS tapers and the time-bandwidth product

The tapers that make this work are the discrete prolate spheroidal sequences (DPSS, or Slepian sequences): for a segment of length T and a chosen half-bandwidth W, they are the sequences that concentrate the largest possible fraction of their energy inside ±W of the analysis frequency. The first taper concentrates the most, the second the next most, and so on; their concentration ratios (eigenvalues) fall toward zero, and the useful ones are those still close to 1.

Two numbers control everything:

• the time-bandwidth product, conventionally written NW or, in MNE’s parameterisation, time_bandwidth = the product of the segment length and the full bandwidth;
• the number of tapers K, which cannot exceed the number of well-concentrated sequences available: the usual rule is K = floor(time_bandwidth − 1).

So the choice is one choice, made twice over. A larger time-bandwidth product allows more tapers and therefore a lower-variance estimate, and it does so by smoothing over a wider band. For an analysis window of length T = n_cycles / f, the full frequency smoothing is time_bandwidth / T = time_bandwidth · f / n_cycles Hz. At 10 Hz with 7 cycles, time_bandwidth = 4 passes about 5.7 Hz and gives 3 tapers; time_bandwidth = 8 passes about 11.4 Hz and gives 7. TODO(confirm): that expression for the full smoothing is this site’s derivation from the window length and the half-bandwidth, checked numerically by the widget’s tests; whether MNE documents the same expression in its own text has not been verified.

Set against Morlet, the difference is where the bandwidth comes from. A Morlet wavelet’s spectral width is fixed by the cycle count, 2 f √(2 ln 2) / n_cycles, and cannot be changed without also changing the temporal width. Multitaper lets you hold the window length and set the bandwidth separately — and then spends the extra bandwidth on variance reduction rather than throwing it away.

Filter-Hilbert

The third route stays in the time domain. Band-pass filter the signal to the band of interest, form the analytic signal with the Hilbert transform, and read off its magnitude and angle sample by sample:

z(t) = x_band(t) + i · H{x_band}(t) = a(t) · e^{iφ(t)}

a(t) is the band-limited amplitude envelope and φ(t) the instantaneous phase. Both are exactly what the burst analysis of L4.6 and the phase analyses of L4.5 need, and neither requires a time-frequency plane to be computed first.

Three properties decide when this is the right tool:

• The filter is explicit and yours to design. Type, order, transition band, causal or zero-phase — all the choices of L1.5, all reportable. That is an advantage over a wavelet, whose filter properties are implied by a cycle count.
• The result is one band. You get a time course, not a map. Asking for a plane means repeating the whole procedure per band, and a bank of narrow filters is slower and harder to report than one wavelet family.
• The narrower the band, the more the output looks like an oscillation whether or not one exists. A 2 Hz-wide filter applied to white noise returns something that rises and falls smoothly at the band centre, because that is what the filter passes. This is pf-narrowband-filter-oscillation, and it is not a subtle effect: a phase-randomised surrogate of a real recording — identical power spectrum, every waveform feature destroyed — comes out of a narrow filter looking exactly as rhythmic as the real thing.

There is a further practical constraint. A zero-phase FIR needs a filter long enough to realise the transition band you asked for, and then padding of a few times its own length; on a short epoch a narrow band cannot be built at all, and a library that silently shortens the filter gives you a wider passband than the one you reported. The widget says when it has had to do that, and by how much.

What is the same, and what is not

At the level of the arithmetic, the three are closer than the vocabulary suggests. Convolving with a complex Morlet wavelet is band-pass filtering with a Gaussian-shaped, analytic filter and then taking magnitude and angle — so “Morlet power” and “the squared envelope from filter-Hilbert with a Gaussian band” are two names for one operation, and a multitaper estimate with a single taper is a third. (Cohen, 2014) makes this equivalence explicitly. Where they differ is in what is easy: Morlet gives the whole plane from one parameter, multitaper gives a stated bandwidth with reduced variance, filter-Hilbert gives one band with an explicit filter and a phase you can use.

What they do not do is return the same numbers. On the shipped epoch, the client’s Morlet map and the precomputed multitaper map of the same 3 seconds of data, on the same frequency grid and with the same window length at every frequency, correlate at r = 0.7099 on log₁₀ power over the region more than one edge margin from either end — and the multitaper estimate is a median factor of 2.40 larger. Two things follow, and both are easy to get wrong:

• The factor is not a scaling mistake. On a pure tone the two estimators agree to within about 13 % (a ratio of 1.1255). The factor on the real epoch is the estimators differing on a broadband input: the multitaper window is several times longer in time and its taper bank sums over a wider band, so a 1/f background contributes more to every one of its estimates.
• The correlation depends on the units, so the units have to be named. Both maps are the analysis library’s own |W ∗ x|², which is proportional to power but carries a frequency-dependent constant. Dividing each frequency row by its own constant to get mean-square µV² shifts that row in log space, and a correlation pooled across rows moves when the rows move relative to each other: the same comparison is r = 0.7099 in library units and r = 0.7001 calibrated, with the median ratio moving from 2.398 to 2.131. Neither is “the” agreement. “The estimators correlate at r = 0.71” is a true sentence only with the units stated.

The moderate correlation is the two estimators disagreeing about where in time power sits, which is what different temporal smoothing means. That is the subject of this lesson, not an error to be tuned away.

Choosing

• A narrow, well-defined low-frequency rhythm — alpha, mu, theta — where you want the time course and, later, the phase: Morlet at a stated cycle count, or filter-Hilbert with a filter you designed. Both are defensible; say which and give the parameters.
• A broad high-frequency band — 60–90 Hz, or any “gamma” range: multitaper. The activity there is broadband rather than a line, a single estimator over a band that wide is noisy, and multitaper lets you state the bandwidth you intend and average several tapers over it. A fixed-cycle Morlet also smooths broadly at high frequency, but its bandwidth is a side effect of the cycle count rather than a stated parameter. And be aware that the scalp carries cranial EMG in that range (pf-muscle-as-gamma), which no estimator fixes.
• Phase — of one band, for coherence, phase-amplitude coupling or an ITC time course at a single frequency: filter-Hilbert, or the angle of the Morlet coefficient, which is the same quantity. Whichever you use, the band is the thing to report, because phase is only defined relative to one.
• A first look at an unfamiliar epoch: Morlet, because one parameter gives the whole plane and the widths are predictable from the formula.

The data behind this lesson

• The same single epoch as L4.2: ds-eegbci S001 R08, channel C3, a right-fist motor-imagery trial, 3 s at the dataset’s native 160 Hz, average-referenced offline, unfiltered. ODC-By 1.0, open access. The trial is the one with the largest mu drop among that subject’s 22 — a selected teaching example, with the full ranking in the sidecar.
• The precomputed multitaper map ships beside it: 37 frequency rows (4–40 Hz), 480 time points, time_bandwidth = 4.0, 3 tapers (read back from MNE’s own taper weights rather than assumed), the same n_cycles = max(3, f/2) scheme as the Morlet map so the window length at each frequency is identical, and the exact mne.time_frequency.tfr_array_multitaper call recorded in the sidecar.
• Both maps carry unit_cosine_power_uv2, one value per frequency row: what the identical call returns for a unit-amplitude cosine at that frequency, measured rather than derived. That is the constant the calibration above divides out.
• The notebook runs all three estimators on the same epochs and reports the agreement, so the numbers in this lesson can be checked rather than taken.

Explore

Set the time-bandwidth product to 4 and look at the three tapers it implies and their concentration ratios; each taper alone is a noisy estimate and the average is not. Raise the product and watch the map smooth across frequency while the bandwidth read-out grows. Switch to filter-Hilbert and confirm that with the passband matched to the wavelet’s own spectral FWHM the two maps agree closely — then narrow the band and watch the widget report that the filter had to be truncated. Finally put the client’s Morlet map beside the stored multitaper map and read the agreement panel, in both unit systems.

Practice

The notebook computes Morlet, multitaper and filter-Hilbert estimates of the same ds-eegbci epochs with the parameters stated in one cell, and prints the correlation and the ratio between them in both library units and calibrated µV².

Exercises