Multitaper Spectral Estimation¶

Multitaper estimation averages spectra across orthogonal DPSS tapers. It can stabilize spectral estimates when trials are noisy or short, at the cost of additional smoothing and computation. It is an estimator choice, not a universally better fit.

In [1]:

import matplotlib
matplotlib.use("module://matplotlib_inline.backend_inline")
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import fftconvolve

from fftrf import TRF, pearsonr, r2_score

import matplotlib matplotlib.use("module://matplotlib_inline.backend_inline") import matplotlib.pyplot as plt import numpy as np from scipy.signal import fftconvolve from fftrf import TRF, pearsonr, r2_score

In [2]:

def simulate_trials(*, seed=7, n_trials=6, n_samples=2048, noise_scale=0.12):
    rng = np.random.default_rng(seed)
    fs = 128.0
    tmin, tmax = 0.0, 0.25
    times = np.arange(round(tmin * fs), round(tmax * fs)) / fs
    true_kernel = (
        0.9 * np.exp(-0.5 * ((times - 0.055) / 0.014) ** 2)
        - 0.5 * np.exp(-0.5 * ((times - 0.125) / 0.024) ** 2)
    )

    stimulus, response = [], []
    for _ in range(n_trials):
        x = rng.standard_normal((n_samples, 1))
        y = fftconvolve(x[:, 0], true_kernel, mode="full")[:n_samples]
        y += noise_scale * rng.standard_normal(n_samples)
        stimulus.append(x)
        response.append(y[:, np.newaxis])
    return fs, tmin, tmax, times, true_kernel, stimulus, response

def simulate_trials(*, seed=7, n_trials=6, n_samples=2048, noise_scale=0.12): rng = np.random.default_rng(seed) fs = 128.0 tmin, tmax = 0.0, 0.25 times = np.arange(round(tmin * fs), round(tmax * fs)) / fs true_kernel = ( 0.9 * np.exp(-0.5 * ((times - 0.055) / 0.014) ** 2) - 0.5 * np.exp(-0.5 * ((times - 0.125) / 0.024) ** 2) ) stimulus, response = [], [] for _ in range(n_trials): x = rng.standard_normal((n_samples, 1)) y = fftconvolve(x[:, 0], true_kernel, mode="full")[:n_samples] y += noise_scale * rng.standard_normal(n_samples) stimulus.append(x) response.append(y[:, np.newaxis]) return fs, tmin, tmax, times, true_kernel, stimulus, response

In [3]:

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials(
    seed=44,
    noise_scale=0.3,
)
train_x, test_x = stimulus[:-1], stimulus[-1]
train_y, test_y = response[:-1], response[-1]

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials( seed=44, noise_scale=0.3, ) train_x, test_x = stimulus[:-1], stimulus[-1] train_y, test_y = response[:-1], response[-1]

Fit a segmented Hann estimate and a DPSS multitaper estimate¶

The same ridge value and segment duration isolate the spectral estimator choice.

In [4]:

hann_model = TRF(direction=1)
hann_model.train(
    train_x, train_y,
    fs=fs, tmin=tmin, tmax=tmax,
    regularization=1e-2,
    segment_duration=1.0,
    overlap=0.5,
    window="hann",
)

multitaper_model = TRF(direction=1)
multitaper_model.train_multitaper(
    train_x, train_y,
    fs=fs, tmin=tmin, tmax=tmax,
    regularization=1e-2,
    segment_duration=1.0,
    overlap=0.5,
    time_bandwidth=3.5,
    n_tapers=4,
)

_, hann_r = hann_model.predict(test_x, test_y)
_, multitaper_r = multitaper_model.predict(test_x, test_y)
print(f"Hann held-out r: {float(hann_r):.3f}")
print(f"multitaper held-out r: {float(multitaper_r):.3f}")

hann_model = TRF(direction=1) hann_model.train( train_x, train_y, fs=fs, tmin=tmin, tmax=tmax, regularization=1e-2, segment_duration=1.0, overlap=0.5, window="hann", ) multitaper_model = TRF(direction=1) multitaper_model.train_multitaper( train_x, train_y, fs=fs, tmin=tmin, tmax=tmax, regularization=1e-2, segment_duration=1.0, overlap=0.5, time_bandwidth=3.5, n_tapers=4, ) _, hann_r = hann_model.predict(test_x, test_y) _, multitaper_r = multitaper_model.predict(test_x, test_y) print(f"Hann held-out r: {float(hann_r):.3f}") print(f"multitaper held-out r: {float(multitaper_r):.3f}")

Hann held-out r: 0.987
multitaper held-out r: 0.944

Compare the estimated kernels¶

Judge the choice using held-out performance and stability across independent observations. A smoother-looking kernel is not by itself evidence of a better model.

In [5]:

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(times * 1e3, true_kernel, "--", color="black", label="True")
ax.plot(hann_model.times * 1e3, hann_model.weights[0, :, 0], label="Hann")
ax.plot(multitaper_model.times * 1e3, multitaper_model.weights[0, :, 0],
        label="Multitaper")
ax.set(xlabel="Lag (ms)", ylabel="Weight", title="Estimator comparison")
ax.legend()
plt.show()

fig, ax = plt.subplots(figsize=(8, 4)) ax.plot(times * 1e3, true_kernel, "--", color="black", label="True") ax.plot(hann_model.times * 1e3, hann_model.weights[0, :, 0], label="Hann") ax.plot(multitaper_model.times * 1e3, multitaper_model.weights[0, :, 0], label="Multitaper") ax.set(xlabel="Lag (ms)", ylabel="Weight", title="Estimator comparison") ax.legend() plt.show()

Tune time_bandwidth and n_tapers together. Larger values increase smoothing. See the Multitaper guide.