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

Question

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

我的问题之后是我的previous question

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

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

delta: 0.5-4 hz theta: 4-8 hz alpha: 8-12 hz beta: 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中，高γ具有最高的振幅（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()

确实如此。看到不同的峰值，正如预期的那样，由于高伽马是我们信号中最强的，我们在100Hz处得到一个尖锐的峰值。我们也看到不同类型的脑电波有不同的峰值。但是，请注意，由于噪声的原因，您的真实世界数据可能不会给出像此数据那样清晰的FFT频谱。您可以通过使这些数据更嘈杂来测试这一点