Frequency-Resolved Kernels¶

This notebook asks which frequencies contribute at different lags of a fitted kernel. A simulation contains an early transient and a time-locked 10 Hz burst, giving each representation a known feature to recover.

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

Simulate a response with a time-locked alpha burst¶

In [2]:

rng = np.random.default_rng(11)
fs = 128.0
tmin, tmax = 0.0, 0.32
times = np.arange(round(tmax * fs)) / fs
true_kernel = (
    0.4 * np.exp(-0.5 * ((times - 0.04) / 0.012) ** 2)
    + 0.55 * np.exp(-0.5 * ((times - 0.16) / 0.055) ** 2)
    * np.cos(2 * np.pi * 10 * (times - 0.16))
)

stimulus, response = [], []
for _ in range(7):
    x = rng.standard_normal((3072, 1))
    y = fftconvolve(x[:, 0], true_kernel, mode="full")[:3072]
    y += 0.12 * rng.standard_normal(3072)
    stimulus.append(x)
    response.append(y[:, None])

train_x, test_x = stimulus[:-1], stimulus[-1]
train_y, test_y = response[:-1], response[-1]

rng = np.random.default_rng(11) fs = 128.0 tmin, tmax = 0.0, 0.32 times = np.arange(round(tmax * fs)) / fs true_kernel = ( 0.4 * np.exp(-0.5 * ((times - 0.04) / 0.012) ** 2) + 0.55 * np.exp(-0.5 * ((times - 0.16) / 0.055) ** 2) * np.cos(2 * np.pi * 10 * (times - 0.16)) ) stimulus, response = [], [] for _ in range(7): x = rng.standard_normal((3072, 1)) y = fftconvolve(x[:, 0], true_kernel, mode="full")[:3072] y += 0.12 * rng.standard_normal(3072) stimulus.append(x) response.append(y[:, None]) train_x, test_x = stimulus[:-1], stimulus[-1] train_y, test_y = response[:-1], response[-1]

Fit and inspect the ordinary kernel first¶

The lag-domain kernel establishes the timing and polarity that the frequency-resolved views will decompose.

In [3]:

model = TRF(direction=1)
model.train(
    train_x, train_y,
    fs=fs, tmin=tmin, tmax=tmax,
    regularization=1e-2,
    segment_duration=2.0,
    overlap=0.5,
    window="hann",
)
_, heldout_r = model.predict(test_x, test_y)
print(f"held-out r: {float(heldout_r):.3f}")

fig, ax = plt.subplots(figsize=(8, 3.5))
ax.plot(times * 1e3, true_kernel, "--", color="black", label="True")
ax.plot(model.times * 1e3, model.weights[0, :, 0], label="Recovered")
ax.set(xlabel="Lag (ms)", ylabel="Weight", title="Ordinary kernel")
ax.legend()
plt.show()

model = TRF(direction=1) model.train( train_x, train_y, fs=fs, tmin=tmin, tmax=tmax, regularization=1e-2, segment_duration=2.0, overlap=0.5, window="hann", ) _, heldout_r = model.predict(test_x, test_y) print(f"held-out r: {float(heldout_r):.3f}") fig, ax = plt.subplots(figsize=(8, 3.5)) ax.plot(times * 1e3, true_kernel, "--", color="black", label="True") ax.plot(model.times * 1e3, model.weights[0, :, 0], label="Recovered") ax.set(xlabel="Lag (ms)", ylabel="Weight", title="Ordinary kernel") ax.legend() plt.show()

held-out r: 0.997

Signed frequency-resolved weights¶

value_mode="real" preserves polarity and cancellation. Use it when the sign of the response matters.

In [4]:

signed = model.frequency_resolved_weights(
    n_bands=20,
    fmax=30.0,
    value_mode="real",
)
model.plot_frequency_resolved_weights(
    resolved=signed,
    title="Signed band-limited kernels",
)
plt.show()

signed = model.frequency_resolved_weights( n_bands=20, fmax=30.0, value_mode="real", ) model.plot_frequency_resolved_weights( resolved=signed, title="Signed band-limited kernels", ) plt.show()

Magnitude of frequency-resolved weights¶

Magnitude discards sign and emphasizes the strength of band-limited structure. It is easier to scan, but positive and negative contributions are no longer distinguishable.

In [5]:

magnitude = model.frequency_resolved_weights(
    n_bands=20,
    fmax=30.0,
    value_mode="magnitude",
)
model.plot_frequency_resolved_weights(
    resolved=magnitude,
    title="Magnitude of band-limited kernels",
)
plt.show()

magnitude = model.frequency_resolved_weights( n_bands=20, fmax=30.0, value_mode="magnitude", ) model.plot_frequency_resolved_weights( resolved=magnitude, title="Magnitude of band-limited kernels", ) plt.show()

Time-frequency power¶

Hilbert power smooths the oscillatory signed weights into a positive envelope. Here it should localize energy near 10 Hz and 160 ms. It describes the fitted kernel—not induced EEG power.

In [6]:

power = model.time_frequency_power(n_bands=20, fmax=30.0)
model.plot_time_frequency_power(
    power=power,
    title="Hilbert power of the fitted kernel",
)
plt.show()

power = model.time_frequency_power(n_bands=20, fmax=30.0) model.plot_time_frequency_power( power=power, title="Hilbert power of the fitted kernel", ) plt.show()

Band count, bandwidth, and frequency range are analysis choices. Prespecify them when confirmatory interpretation is intended. See the Frequency-Resolved guide.