python读取BCI的gdf数据

        被老师坑了做这个项目,要用脑电分析,现在共享一下一些文档,有好心人告诉我一下怎么转码吗?

# coding=utf8
## 进行带通滤波,并进行小波变换,画出原始的小波图像
# 可以选择的是小波的数据的显示
# 还要注意小波图的大小
# CWT可以完成很多事情,但是效果不一定好
import mne
import matplotlib.pyplot as plt
import numpy as np
import pywt
import scipy.io
from scipy import signal
from sklearn import preprocessing
import os
import gc
import cv2
import glob
import scipy.io as scio

eventDescription = {'276': "eyesOpen", '277': "eyesClosed", '768': "startTrail", '769': "cueLeft",
                    '770': "cueRight", '781': "feedback", '783': "cueUnknown",
                    '1023': "rejected", '1077': 'horizonEyeMove', '1078': "verticalEyeMove",
                    '1079': "eyeRotation", '1081': "eyeBlinks", '32766': "startRun"}

ch_types = ['eeg', 'eeg', 'eeg', 'eog', 'eog', 'eog']
ch_names = ['EEG_Cz', 'EEG_C3', 'EEG_C4', 'EOG_ch01', 'EOG_ch02', 'EOG_ch03']
label_name = ["left", "right"]  # 标签
time_label = ["4_55", "55_7"]


def butter_bandpass_filtfilt(data, cutoff, f_low, f_high, order=8):  # 巴特沃得滤波
    # data 输入
    # cutoff  采集频率
    # f_low 下限频率
    # f_high 上限频率
    # order 阶数
    w_down = 2 * f_low / cutoff
    w_up = 2 * f_high / cutoff
    b, a = signal.butter(order, [w_down, w_up], 'bandpass', analog=False)
    output = signal.filtfilt(b, a, data, axis=0)
    return output


if __name__ == "__main__":
    wavlist = pywt.wavelist(kind='continuous')
    for j in range(1, 10):
        base_line = None

        # base_line = (None,None)
        #
        wavename = 'cmor1-1'
        sampling_rate = 250
        totalscal = 600
        fc = pywt.central_frequency(wavename)
        cparam = 2 * fc * totalscal
        scales = cparam / np.arange(totalscal, 1, -1)
        left_num = 0
        right_num = 0
        gc.enable()
        data_list = glob.glob(f"..\\source_data\\BCICIV_2b_gdf\\B0" + str(j) + "0*.gdf")  # 实验数据
        lable_path = "..\\source_data\\BCICIV_2b_gdf\\true_labels\\"
        save_path = "..\\PIC_2\\0" + str(j) + "\\CWT_PIC\\" + wavename + "_minmax_noeq_nobaseline_300x100" + "\\"
        for flie in data_list:
            print(flie)
            file_name = flie.split("\\")[-1].split(".")[0]
            file_type = file_name[-1]
            rawDataGDF = mne.io.read_raw_gdf(flie, preload=True, eog=['EOG:ch01', 'EOG:ch02', 'EOG:ch03'])
            # 创建数据的描述信息
            info = mne.create_info(ch_names=ch_names, sfreq=rawDataGDF.info['sfreq'], ch_types=ch_types)
            # 创建数据结构体
            data = np.squeeze(np.array(
                [rawDataGDF['EEG:Cz'][0], rawDataGDF['EEG:C3'][0], rawDataGDF['EEG:C4'][0], rawDataGDF['EOG:ch01'][0],
                 rawDataGDF['EOG:ch02'][0], rawDataGDF['EOG:ch03'][0]]))
            # 创建RawArray类型的数据
            rawData = mne.io.RawArray(data, info)
            # 获取事件
            event, _ = mne.events_from_annotations(rawDataGDF)
            event_id = {}
            for i in _:  # 整理event_id
                event_id[eventDescription[i]] = _[i]
            # 提取epoch
            epochs4_55 = mne.Epochs(rawData, event, event_id, tmin=4, tmax=5.5, baseline=base_line,
                                    event_repeated='merge')
            epochs55_7 = mne.Epochs(rawData, event, event_id, tmin=5.5, tmax=7, baseline=base_line,
                                    event_repeated='merge')
            ##得到epochs,可以得到事件的列表。
            for time in time_label:
                if time == "4_55":
                    epochs = epochs4_55
                if time == "55_7":
                    epochs = epochs55_7

                if file_type == "T":
                    ev_left = epochs['cueLeft']
                    ev_right = epochs['cueRight']
                    left_data = ev_right.get_data()
                    right_data = ev_left.get_data()
                if file_type == "E":
                    left_data = []
                    right_data = []
                    ev_Unknown = epochs['cueUnknown']
                    unknow_data = ev_Unknown.get_data()
                    labels = scio.loadmat(lable_path + file_name + ".mat")["classlabel"]
                    j = 0
                    for label in labels:
                        if label == 1:
                            left_data.append(unknow_data[j])
                        if label == 2:
                            right_data.append(unknow_data[j])
                        j = j + 1

                # 得到左手或者右手的数据。
                left_data = np.array(left_data)
                right_data = np.array(right_data)

                # 小波的系数
                for trip_label in label_name:
                    if trip_label == "left":
                        data = left_data
                        num = left_num
                    if trip_label == "right":
                        data = right_data
                        num = right_num
                    # 先对左手进行滤波加进行小波变换,然后把图存放在指定的文件夹中  \\left\\通道名称
                    [trip_num, channel_num, data_num] = data.shape
                    for i in range(trip_num):
                        num = num + 1
                        for j in range(channel_num):
                            print(trip_label, "prossing", num)
                            print("channle", ch_names[j])
                            tmp_data = data[i, j, :]  # 取出一个通道的数据
                            tmp_data_filter = butter_bandpass_filtfilt(tmp_data, 250, 7, 40)  # 进行带通滤波 7-30Hz 8阶滤波器
                            min_max_scaler = preprocessing.MinMaxScaler()  # 默认为范围0~1,拷贝操作
                            tmp_data_filter_min_max = np.squeeze(
                                min_max_scaler.fit_transform(tmp_data_filter.reshape(-1, 1)))  # 变成0-1之间进行归一化。

                            [cwtmatr, frequencies] = pywt.cwt(tmp_data_filter_min_max, scales, wavename,
                                                              1.0 / sampling_rate)  # 进行小波变换
                            # 参数需要选择,进行比较 小波出来是个复数,也需要继续进行分类
                            min_frequency = np.sum(frequencies >= 8)
                            max_frequency = frequencies.shape[0] - np.sum(frequencies <= 38)  # 要选频段,不能全部进去
                            # 计算需要的频带,我们需要的在[max_frequency:min_frequency]之间
                            abs_cwtmatr = abs(cwtmatr[max_frequency:min_frequency])  # 这里是求绝对值,如果小波基不一样就要改一下
                            abs_cwtmatr_min_max = min_max_scaler.fit_transform(abs_cwtmatr)  # 取最大最小值,让图像明显一点
                            pic_data = np.array(abs_cwtmatr_min_max * 255, dtype='uint8')  # 直接转换成图像
                            # pic_data_equalizeHist = cv2.equalizeHist(pic_data)  # 直方图均衡化,让图片更明显一点

                            path_saving = save_path + trip_label + "\\" + ch_names[j]  # 保存的目录
                            if not os.path.exists(path_saving):  # 创建目录
                                os.makedirs(path_saving)

                            pic_data = cv2.resize(pic_data, (300, 100), interpolation=cv2.INTER_AREA)

                            cv2.imwrite(path_saving + "\\" + time + "_" + str(num) + ".png", pic_data)  # 保存图

                    if trip_label == "left":
                        left_num = num
                    if trip_label == "right":
                        right_num = num

 

  • 0
    点赞
  • 12
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值