python np fft_用python对航弈生物BCIduino脑电

用python对航弈生物BCIduino放大器脑电/肌电数据进行滤波及频谱计算（一）

这篇稿子介绍如何用pythonmne对fif格式保存的脑电数据进行读取和简单的滤波，并用numpy对上述数据进行FFT分析，绘制相应的图像。本稿不涉及BCIduino数据如何保存，关于数据如何保存会在之后的文章给出。本稿先给出一个已经保存好的.fif格式的BCIduino肌电数据，数据保存了大约40s的肌电信号，信号内容是正常人模拟帕金森患者震颤时的肌电信号，采样频率250Hz，一共采集了5次，每次4s，对应fif文件中的标签11.下文会给出代码并进行解释。

1. 读取数据及利用mne简单滤波——

关键代码：

importmne

raw = mne.io.read_raw_fif('2019_12_26_09_44_22-raw.fif')

前提是安装了mne库，安装方法是：在dos窗口里运行pip installmne，如下：

raw这个变量就是获取到的原始数据，如果需要了解raw的信息，可以多运行一行

print(raw)，可以看到dos窗口会打印出下列信息：

可以看出这段数据一共45.6s，一共9个通道（其中第0个是event通道，第1-8个是原始数据通道），一共11401个采样点，采样频率我们已经知道是250Hz，可以算出来确实是45.6s的数据。

把原始数据存到变量里之后，我们可以找一下events（即标签），运行：

events= mne.find_events(raw,output='step',consecutive=True, shortest_event=1, uint_cast=True )

此时dos窗口输出：

我们关心的数据是标签11对应的数据，其他标签对应的数据对我们意义不大。因此要找出相应的数据段。

epochs = mne.Epochs(raw,events,event_id = 11,tmin = 0,tmax = 4)

上述代码语句将raw（原始数据）中event_id=11的数据提取出（一共提取了4s，从标签11开始出现开始截取数据，一共截取4s，这4s的数据是有意义的）。

如果将epoch打印出，即print（epochs），可以看到一共有五段数据，因为我们一共采集了5次数据，每次4s，如下：

然后将epoch数据转化为evoked数据（方便处理）：

evoked = epochs.average()

相当于将epochs的5段数据进行了平均。

然后对evoked数据进行带通滤波：

evoked.filter(l_freq=4,h_freq=8.0)

并对evoked数据绘图：

evoked.plot(titles='EMG',window_title='EMG(Hangyi Biotech)')

效果如下图：

2. 对evoked数据进行FFT计算

该部分借鉴https://blog.csdn.net/qq_39516859/article/details/79794549博客，以下内容带有详细注释，有问题可以加入文末微信群或QQ群讨论交流。

x = evoked._data[0]

sr_fft = {250:128,500:256}

fft_size = sr_fft[sampling_rate] #FFT处理的取样长度

t = np.arange(0, 1.0, 1.0/sampling_rate)#np.arange(起点，终点，间隔)产生1s长的取样时间

xs = x[:fft_size]# 从波形数据中取样fft_size个点进行运算

xf = np.fft.rfft(xs)/fft_size# 利用np.fft.rfft()进行FFT计算，rfft()是为了更方便对实数信号进行变换，由公式可知/fft_size为了正确显示波形能量

# rfft函数的返回值是N/2+1个复数，分别表示从0(Hz)到sampling_rate/2(Hz)的分。

#于是可以通过下面的np.linspace计算出返回值中每个下标对应的真正的频率：

freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)

# np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)

#在指定的间隔内返回均匀间隔的数字

xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))

#最后我们计算每个频率分量的幅值，并通过 20*np.log10()将其转换为以db单位的值。为了防止0幅值的成分造成log10无法计算，我们调用np.clip对xf的幅值进行上下限处理

arr_