用枕建立脑电波的时频表示

Question

我想要创建脑电波的时频表示(由MEG神经成像获得的神经活动数据)。

我的previous question跟进了我的问题

我需要的是脑电波不同频带的时频表示(振幅值或功率值)。

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

δ: 0.5-4 hzθ: 4-8 hzα: 8-12 hzβ: 12-30 hz low_gamma 30-80 hz high_gamma = 80 - 120 hz

我来自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

正如我们在对您的older question的注释中所讨论的，CWT的Y轴只是基于某个小波的输入波的缩放版本，其中X轴是原始波上的时间点。小波变换只是小波和原始波的卷积。正如Tim在评论中提到的，您可能正在寻找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是我们信号中最强的，我们在100 we时得到了一个尖峰。对于不同类型的脑电波，我们也看到了不同的峰值。然而，请注意，您的真实世界的数据可能不会给出一个清晰的FFT频谱，因为这是因为噪音。您可以通过增加这些数据的噪音来测试这些数据。