使用 scipy 创建脑电波的时频表示

Question

我想创建脑电波的时频表示（通过 MEG 神经成像获得的神经活动数据）。

我的问题由我的previous question跟进

我需要的是脑电波不同频段的时频表示（幅度值或功率值）。

更具体地说，脑电波的频带定义如下：

增量：0.5-4 赫兹 θ：4-8赫兹 阿尔法：8-12 赫兹 测试版：12-30 赫兹 low_gamma 30-80 赫兹 high_gamma = 80 - 120 赫兹

previous thread 中我的（固定）代码：

import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal


N = 5000
rnd = np.random.RandomState(12345)

brain_signal = np.sin(np.linspace(0, 1000, N)) + rnd.uniform(0, 1, N)
widths = np.arange(1, N//8)
cwtmatr = signal.cwt(brain_signal, signal.ricker, widths)

fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(10, 6))
axes = ax.flatten()

sns.lineplot(np.linspace(0, 1000, N), brain_signal, ax=axes[0], lw=2)
sns.heatmap(cwtmatr, cmap='Spectral', ax=axes[1]);

axes[0].set_title('Brain signal')
axes[1].set_title('CWT of brain signal')

plt.tight_layout()

我的问题是，y 轴上的值是否正确显示频率信息？如果没有，我应该如何解决？

更新： 上述代码的结果：

提前致谢

Answer 1

就像我们在 cmets 上向您的 older question 讨论的那样，CWT 的 Y 轴只是基于某个小波的输入波的缩放版本，其中 X 轴是原始波上的时间点。小波变换只是小波和原始波的卷积。正如 Tim 在 cmets 中提到的，您可能正在寻找 FFT。这是一个例子：

import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq

# Number of data points    
N = 5000
rnd = np.random.RandomState(12345)
# Sampling interval
T = 1/1000

我没有你的数据，所以我只是将你上面的频率用于不同的大脑信号来制作一个虚构的大脑信号。我还添加了一些噪点以使其随机且更逼真。

## Frequencies
delta = 2
theta = 5
alpha = 10
beta = 20
low_gamma = 75
high_gamma = 100


x = np.linspace(0, N*T, N)
brain_signal = np.sin(delta * 2.0*np.pi*x) + \
               np.sin(theta * 2.0*np.pi*x) + \
               4 *np.sin(alpha * 2.0*np.pi*x) + \
               np.sin(beta * 2.0*np.pi*x) + \
               2 * np.sin(low_gamma * 2.0*np.pi*x) + \
               10 * np.sin(high_gamma * 2.0*np.pi*x) + rnd.uniform(-10, 10, N)

请注意，在此brain_signal 中，high_gamma 具有最高幅度 (10)。现在让我们看看 FFT 是否可以从一个信号中恢复任何这些不同类型的波。

yf = fft(brain_signal)
xf = fftfreq(N, T)[:N//2]

fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(10, 6))
axes = ax.flatten()

sns.lineplot(x, brain_signal, ax=axes[0], lw=2)
sns.lineplot(xf[:1000], 2.0/N * np.abs(yf[0:N//2])[:1000], ax=axes[1]);


axes[0].set_title('Brain signal \nTime domain')
axes[0].set_xlabel('Time')
axes[0].set_ylabel('Amplitude')

axes[1].set_title('FFT of brain signal \nFrequency domain')
axes[1].set_xlabel('Frequency')
axes[1].set_ylabel('Amplitude')

plt.tight_layout()

确实如此。看到不同的峰值，正如预期的那样，由于 high_gamma 是我们信号中最强的，我们在 100Hz 处得到一个尖峰。我们还看到不同类型脑电波的不同峰值。但是，请注意，由于噪声，您的真实世界数据可能无法提供清晰的 FFT 频谱。您可以通过使这些数据更嘈杂来测试它。