脑地形图动态显示

• 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中的图片还有些是黑的

所以建议大家还是按照第一种方法来吧🤦‍,好奇怪，gif不能一直循环，那就看下图吧

Addition：这是将gif上传到网站上并可以获得一个Link的工具，但我一直上传失败GIPHY