标签：动态显示 maxlevel 地形图 iter path import eeg data
- OpenBCI将脑地形图以动态的gif方式显示出来，首先运行topographic.py将png保存下来，然后运行gif_generate.py将png做成gif，这是第一种方法，刚刚思考了一下，应该还有第二种方法，直接在gif_generate.py中直接调用RealtimeExtract.create_epoch()然后生成png，一生成png就把它保存下来，这样就不会占文件空间了。
- 这里顺带提一句matplotlib.use('agg')如果没有这句代码，plt.savefig()是无法将脑地形图保存下来了，真是惨痛的经历！
Method One
gif_generate.py
from PIL import Image
import os
path='image/'
im=Image.open(path+"0.png")
images=[]
for root,dir,files in os.walk(path):
for file in files:
images.append(Image.open(os.path.join(root,file)))
im.save('gif.gif', save_all=True, append_images=images,loop=1,duration=1)
topograpich map.py
import time
import numpy as np
import matplotlib
import mne
import brainflow
from brainflow.data_filter import DataFilter, FilterTypes, AggOperations, WindowFunctions
from brainflow.board_shim import BoardShim, BrainFlowInputParams, BoardIds
import mne
matplotlib.use('agg')
import matplotlib.pyplot as plt
import pywt
from RealtimeExtract import create_epoch
def plot_sth():
for i in range(10):
epochs, eeg_data = create_epoch(Time=2, Duration=1., Overlap=0.0)
epochs.set_montage('standard_1020')
epochs.plot_psd_topomap()
plt.savefig('image/'+str(i))
plot_sth()
RealtimeExtract.py
import time
import numpy as np
import matplotlib
import mne
import brainflow
from brainflow.data_filter import DataFilter, FilterTypes, AggOperations, WindowFunctions
from brainflow.board_shim import BoardShim, BrainFlowInputParams, BoardIds
import mne
import matplotlib.pyplot as plt
import pywt
import pandas as pd
from scipy.fftpack import fft
'''
Update:2020/12/02
Extract the alpha delta and so on of each second
'''
def create_epoch(Time,Duration,Overlap):
BoardShim.enable_dev_board_logger()
# use synthetic board for demo
params = BrainFlowInputParams()
board = BoardShim(BoardIds.SYNTHETIC_BOARD.value, params)
board.prepare_session()
board.start_stream()
time.sleep(Time)
data = board.get_board_data()
board.stop_stream()
board.release_session()
#不知为啥这里是16通道，ValueError: len(data) (16) does not match len(info["ch_names"]) (8)
eeg_channels = [1, 2, 3, 4, 5, 6, 7, 8]
eeg_data = data[eeg_channels, :]
eeg_data = eeg_data / 1000000 # BrainFlow returns uV, convert to V for MNE
# Creating MNE objects from brainflow data arrays
# 记得对应的通道 TP9，TP10，AF7，AF8
print(BoardShim.get_eeg_names(BoardIds.SYNTHETIC_BOARD.value))
ch_types = ['eeg', 'eeg', 'eeg', 'eeg', 'eeg', 'eeg', 'eeg', 'eeg']
ch_names = ['T7', 'CP5', 'FC5', 'C3', 'C4', 'FC6', 'CP6', 'T8']
sfreq = BoardShim.get_sampling_rate(BoardIds.SYNTHETIC_BOARD.value)
print(sfreq)
info = mne.create_info(ch_names=ch_names, sfreq=sfreq, ch_types=ch_types)
# 创建raw对象
#print('raw info:',eeg_data)
raw = mne.io.RawArray(eeg_data, info)
#raw.plot_psd(average=True)
plt.show()
'''
2021.1.30
保存为csv
'''
# raw = mne.io.read_raw_fif('raw.fif')
# 创建event对象,重叠时间overlap:2s
# 创建等距事件，默认id为1
events = mne.make_fixed_length_events(raw, duration=Duration, overlap=Overlap)
# 将事件和原始数据一起绘制
epochs = mne.Epochs(raw, events)
return epochs,eeg_data
# 小波包计算各个频段的能量分布
def WPEnergy(data, fs, wavelet, maxlevel):
# 小波包分解
wp = pywt.WaveletPacket(data=data, wavelet=wavelet, mode='symmetric', maxlevel=maxlevel)
db1 = pywt.Wavelet('db4')
print(len(pywt.wavedec(data, db1)))
# 频谱由低到高的对应关系，这里需要注意小波变换的频带排列默认并不是顺序排列，所以这里需要使用’freq‘排序。
freqTree = [node.path for node in wp.get_level(maxlevel, 'freq')]
# 计算maxlevel最小频段的带宽
freqBand = fs / (2 ** maxlevel)
# 定义能量数组
energy = []
# 循环遍历计算四个频段对应的能量
for iter in range(len(iter_freqs)):
iterEnergy = 0.0
for i in range(len(freqTree)):
# 第i个频段的最小频率
bandMin = i * freqBand
# 第i个频段的最大频率
bandMax = bandMin + freqBand
# 判断第i个频段是否在要分析的范围内
if (iter_freqs[iter]['fmin'] <= bandMin and iter_freqs[iter]['fmax'] >= bandMax):
# 计算对应频段的累加和
iterEnergy += pow(np.linalg.norm(wp[freqTree[i]].data, ord=None), 2)
# 保存四个频段对应的能量和
energy.append(iterEnergy)
# 绘制能量分布图
plt.plot([xLabel['name'] for xLabel in iter_freqs], energy, lw=0, marker='o',color='black')
plt.title('能量分布')
plt.show()
# 小波包变换-重构造分析不同频段的特征
def TimeFrequencyWP(data, fs, wavelet, maxlevel):
# 小波包变换这里的采样频率为250，如果maxlevel太小部分波段分析不到
wp = pywt.WaveletPacket(data=data, wavelet=wavelet, mode='symmetric', maxlevel=maxlevel)
# 频谱由低到高的对应关系，这里需要注意小波变换的频带排列默认并不是顺序排列，所以这里需要使用’freq‘排序。
freqTree = [node.path for node in wp.get_level(maxlevel, 'freq')]
# 计算maxlevel最小频段的带宽
freqBand = fs/(2**maxlevel)
#######################根据实际情况计算频谱对应关系，这里要注意系数的顺序
# 绘图显示
fig, axes = plt.subplots(len(iter_freqs)+1, 1, figsize=(10, 7), sharex=True, sharey=False)
# 绘制原始数据
axes[0].plot(data,color='black')
axes[0].set_title('原始数据')
color=['#33ffff','red','blue','green','black']
for iter in range(len(iter_freqs)):
# 构造空的小波包
new_wp = pywt.WaveletPacket(data=None, wavelet=wavelet, mode='symmetric', maxlevel=maxlevel)
for i in range(len(freqTree)):
print(len(freqTree))
# 第i个频段的最小频率
bandMin = i * freqBand
# 第i个频段的最大频率
bandMax = bandMin + freqBand
# 判断第i个频段是否在要分析的范围内
if (iter_freqs[iter]['fmin']<=bandMin and iter_freqs[iter]['fmax']>= bandMax):
# 给新构造的小波包参数赋值
new_wp[freqTree[i]] = wp[freqTree[i]].data
# 绘制对应频率的数据
axes[iter+1].plot(new_wp.reconstruct(update=True)[:2560],color=color[iter])
# 设置图名
axes[iter+1].set_title(iter_freqs[iter]['name'])
plt.show()
if __name__ == "__main__":
# plot_sth()
# create_raw_and_events()
print("-------------")
'''
CWT to find delta theta alpha and beta
'''
# 获取四个通道的数据并用小波变换的方式找
iter_freqs = [
{'name': 'Delta', 'fmin': 0, 'fmax': 4},
{'name': 'Theta', 'fmin': 4, 'fmax': 8},
{'name': 'Alpha', 'fmin': 8, 'fmax': 13},
{'name': 'Beta', 'fmin': 13, 'fmax': 35},
{'name': 'Gamma', 'fmin': 35, 'fmax': 100},
]
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
mne.set_log_level(False)
#epochsCom = mne.read_epochs(r'raw-epo.fif')
#while(True):
epochsCom, eeg_data = create_epoch(Time=10, Duration=1., Overlap=0.0)
# 三层list get_data()[0][i]其中i是通道口
print(epochsCom)
# print(epochsCom[0].get_data()[0])
# dataCom = epochsCom[0].get_data()[0][0]
# 这里获取的是第一个通道
dataCom = eeg_data[0]
#print(dataCom)
#WPEnergy(dataCom, fs=250, wavelet='db4', maxlevel=4)
TimeFrequencyWP(dataCom, fs=256, wavelet='db4', maxlevel=8)
'''
先面临以下问题：
① 256Hz获取用户1s的脑电信号并进行小波变化后，里面有256个样本，那如何得到1s中每个波段的值呢，每个样本取平均吗？
② 针对多通道的脑机接口，1s内的各个波值又该如何确定呢？
'''
Method Two
from PIL import Image,ImageFile
import os
from RealtimeExtract import create_epoch
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
global im
path='image/'
images=[]
'''
0305 have ready uploaded to cnblog
Method One
for root,dir,files in os.walk(path):
for file in files:
images.append(Image.open(os.path.join(root,file)))
'''
'''
Method Two
'''
new_path=path+'0.png'
for i in range(4):
epochs, eeg_data = create_epoch(Time=2, Duration=1., Overlap=0.0)
epochs.set_montage('standard_1020')
epochs.plot_psd_topomap()
plt.savefig(new_path)
if i==0:
global im
im = Image.open(new_path)
else:
images.append(Image.open(new_path))
im.save('1.gif', save_all=True, append_images=images,loop=1,duration=1)
但这个方法行不通，因为图片损坏了会报错OSError: unrecognized data stream contents when reading image file
加上以下两句代码
from PIL import Image,ImageFile
ImageFile.LOAD_TRUNCATED_IMAGES = True
加了代码之后虽然可以运行了，但还是有点小问题，不知道是不是图片产生和被替代的时间太快了，保存下来的gif中的图片还有些是黑的
所以建议大家还是按照第一种方法来吧
标签：动态显示,maxlevel,地形图,iter,path,import,eeg,data 来源： https://www.cnblogs.com/vincyicy/p/Vincy_King.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。