脑机接口基础学习18---求大佬拯救，如何生成模拟原始脑电数据

最新推荐文章于 2024-06-10 12:24:26 发布

多欢喜 

最新推荐文章于 2024-06-10 12:24:26 发布

阅读量509

点赞数 1

分类专栏：脑机接口基础文章标签：脑机接口

本文链接：https://blog.csdn.net/duohuanxi/article/details/109604889

版权

脑机接口基础专栏收录该内容

20 篇文章 50 订阅

订阅专栏

在实验中，有时需要原始脑电数据来进行模拟实验，但又限于实验条件的不足，需要构造原始的脑电数据
这篇文章通过多次重复所需的源激活来生成原始数据

#导入工具包

import numpy as np
import matplotlib.pyplot as plt

import mne
from mne import find_events,Epochs,compute_covariance,make_ad_hoc_cov
from mne.datasets import sample
from mne.simulation import (simulate_sparse_stc,simulate_raw,
                           add_noise,add_ecg,add_eog)

#导入数据
data_path=r'E:\脑机接口资料\MNE-sample-data\MEG\sample'
raw_fname=data_path+'/sample_audvis_raw.fif'
fwd_fname=data_path+'/sample_audvis-meg-eeg-oct-6-fwd.fif'

#加载真实数据作为模板,主要是为了使用电极位置等一系列相关信息
raw=mne.io.read_raw_fif(raw_fname)
#设置参考电极，默认为average，这里projection为True说明只是添加了参考，但没有进行运算
raw.set_eeg_reference(projection=True)

在这里插入图片描述

##生成偶极子时间序列

#设置偶极子的数量
n_dipoles=4
#每个epoch或event的时间窗口长度
epoch_duration=2.0
#谐波数
n=0
#随机状态（可复制），主要是为了使每次运行的结果相同，方便调试，
#如果不设置为固定状态，每次的结果会不一样
rng=np.random.RandomState(0)

#数据生成函数，在10Hz谐波下产生时间交错的正弦波，函数中的n为全局变量，
#这里有4个偶极子会调试4遍
#也就是生成10Hz，20Hz，30Hz，40Hz的正弦波
def data_fun(times):
    '''产生时间交错的正弦波，弦波为10Hz'''
    global n#全局变量
    n_samp=len(times)#时间序列长度
    window=np.zeros(n_samp)#时间窗口，主要为了后面生成时间交错的波形
    start,stop=[int(ii*float(n_samp)/(2*n_dipoles))
               for ii in (2*n,2*n+1)]#时间窗的起始和结束范围，谐波间隔出现
    window[start:stop]=1.#对想保留的范围赋值为1
    n+=1#全局变量加一
    data=25e-9*np.sin(2.*np.pi*10.*n*times)#正弦波函数为sin(2*pi*w*t)
    data*=window#保留时间窗内对应的数据，也可以理解为滤波
    return data#返回生成的波形数据
    

#获取时间序列，前面我们读取了sample_audivs_raw.fif作为模板，获取其采样频率，其对应的epoch_duration的（2秒钟）的时间序列如下
times=raw.times[:int(raw.info['sfreq']*epoch_duration)]
#获取位置正确信息，即信号源的位置信息
fwd=mne.read_forward_solution(fwd_fname)
#获取源位置
src=fwd['src']
#利用data_fun函数生成4个偶极子离散源的时间序列
stc=simulate_sparse_stc(src,n_dipoles=n_dipoles,times=times,
                       data_fun=data_fun,random_state=rng)



#Look at our source data
fig,ax=plt.subplots(1)
ax.plot(times,1e9*stc.data.T)
ax.set(ylabel='Amplitude (nAm)',xlabel='Time (sec)')
mne.viz.utils.plt_show()

在这里插入图片描述

#####得到模拟原始数据并绘制

'''
模拟原始数据

'''
# 模拟原始数据，[stc]*10指的是时间序列，上面设置的[stc]长度为2秒（epoch_duration）
# raw.info中的events有些是在两秒后的,这里相当于模拟了10个epoch长的的数据也就是20秒的数据
raw_sim = simulate_raw(raw.info, [stc] * 10, forward=fwd,verbose=True)
# 生成噪声的协方差
cov = make_ad_hoc_cov(raw_sim.info)
# 添加噪音
add_noise(raw_sim, cov, iir_filter=[0.2, -0.2, 0.04], random_state=rng)
# 添加心电（ECG）信号
add_ecg(raw_sim, random_state=rng)
# 添加眼电（EOG）信号
add_eog(raw_sim, random_state=rng)
# 绘制添加噪声后的数据
raw_sim.plot()
plt.show()

在这里插入图片描述
为什么我生成的图是这种呢，示例中的图应该是这样的

'''
绘制诱发数据

'''

# 提取事件
events = find_events(raw_sim)
# 根据事件生成epochs
epochs = Epochs(raw_sim, events, 1, tmin=-0.2, tmax=epoch_duration)
# 计算epoch的经验协方差
cov = compute_covariance(epochs, tmax=0., method='empirical', verbose='error')
# epoch求取均值得到evoke
evoked = epochs.average()
# 绘制evoke白化后的数据
evoked.plot_white(cov, time_unit='s')
plt.show()