使用 scipy

编程问答 2022-09-07

问题描述

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

我的问题由我的 previous question

跟进

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

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

增量：0.5-4赫兹 θ：4-8 赫兹 阿尔法：8-12赫兹 测试版：12-30 赫兹 低伽玛 30-80 赫兹 高伽玛 = 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,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 轴上的值是否正确显示了频率信息？如果没有，我该如何解决？

更新： 上面代码的结果：

enter image description here

提前致谢

解决方法

就像我们在对您的 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()

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

numpypythonscipyscipysignal-processingwavelet-transform