我试着用Python语言对EEG信号进行FFT,然后根据带宽确定它是alpha信号还是beta信号。它看起来很好,但得到的曲线图与它们应该的完全不同,频率和震级值并不是我所期望的。如果有人愿意帮忙,代码如下:
from scipy.io import loadmat
import scipy
import numpy as np
from pylab import *
import matplotlib.pyplot as plt
eeg = loadmat("eeg_2013.mat");
eeg1=eeg['eeg1'][0]
eeg2=eeg['eeg2'][0]
fs = eeg['fs'][0][0]
fft1 = scipy.fft(eeg1)
f = np.linspace (fs,len(eeg1), len(eeg1), endpoint=False)
plt.figure(1)
plt.subplot(211)
plt.plot (f, abs (fft1))
plt.title ('Magnitude spectrum of the signal')
plt.xlabel ('Frequency (Hz)')
show()
plt.subplot(212)
fft2 = scipy.fft(eeg2)
f = np.linspace (fs,len(eeg2), len(eeg2), endpoint=False)
plt.plot (f, abs (fft2))
plt.title ('Magnitude spectrum of the signal')
plt.xlabel ('Frequency (Hz)')
show()
如果采样频率是},频率的后半部分对应于镜像频率范围
fs
,并且有N=len(eeg1)
个样本,那么fft
过程当然会返回一个N
值的数组。其中第一个N/2
对应于频率范围{-fs/2..0
。对于实际输入信号,镜像半部分只是正半部分的复共轭,因此在进一步的分析中可以忽略它(但在逆fft中则不能)。在所以本质上,你应该格式化
f=linspace(0,N-1,N)*fs/N
编辑:或者更简单,只需对初始代码进行最小的更改
f = np.linspace (0,fs,len(eeg1), endpoint=False)
因此
f
的范围从0
到fs
之前,忽略输出中fft结果的后半部分:plt.plot( f(0:N/2), abs( fft1(0:N/2) ) )
补充:您可以使用fftshift来交换两半,那么正确的频率范围是
f = np.linspace (-fs/2,fs/2,len(eeg1), endpoint=False)
为了得到fft频率的数组,您应该使用fftfreq;它提供了一个用作横坐标的频率数组:
抱歉,我无法在实际情况下测试此代码,因为您没有发布数据示例,但我希望这应该可以工作。在
关于如何确定带宽,据我所知,你想得到基频。有不同的方法,或多或少复杂,不管你的信号是否有噪音。。。在你的例子中,你只想知道基频f0是在8-13Hz(alpha)还是13-30Hz(beta)范围内;一个非常简单的方法是计算fft在8-13Hz范围内的最大值:
fft1[(f>8) & (f<13)].max()
,如果它超过,比如说1000,那就是α波,否则就是β波。如果你的信号不太相似,请张贴一些不同种类的样本和结果的例子,以便我们可以尝试更复杂的算法。在相关问题 更多 >
编程相关推荐