devbox
木道寻08
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

EEG代码实践：数据集特征提取方法一览（以SEED为例）_seed数据集

作者：木道寻08 | 2024-07-27 15:14:37

赞

踩

seed数据集

2023/8/12 -  2023/8/15 脑机接口学习内容一览：

        本文主要针对上海交通大学的SJTU（SEED）数据集，在前文SEED数据集批量读取的基础上进行闭门造车的特征提取，旨在之后通过更加标准的流程实现分类，顺便对以前学到的东西做一下梳理。

        tips.文章中代码包含库如下：

  1. import mne
  2. import numpy as np
  3. # 引入前文中构建的文件读取函数
  4. from read_file import read_all_files

1. 频谱特征

        操作主要是将EEG信号转换为频谱，提取不同频带的能量，如"delta": [0.5, 4.5], "theta": [4.5, 8.5], "alpha": [8.5, 11.5], "sigma": [11.5, 15.5], "beta": [15.5, 30]等。

        提取的频谱特征矩阵shape如下：

read_all_files start...
共读取了1个文件
共有15段数据
read ended...
频谱特征矩阵的shape为：(15, 310)

        呈现出这个shape很好理解，这里我们选取了5个频段，每个raw数据有62个通道，每个通道的数据都能计算一次频段能量，因此得到数据5*62=310个，代码如下：

  1. def frequency_spectrum(raws):
  2.     # 提取EEG数据在五个频段的能量特征
  3.     # delta(0.5-4Hz) theta(4-8Hz) alpha(8-13Hz) beta(13-30Hz) gamma(30-100Hz)
  4.     # 特定频带
  5.     FREQ_BANDS = {"delta": [0.5, 4.5],
  6.                   "theta": [4.5, 8.5],
  7.                   "alpha": [8.5, 11.5],
  8.                   "sigma": [11.5, 15.5],
  9.                   "beta": [15.5, 30]}
  10.     # 特征矩阵
  11.     feature_matrix = []
  12.     # 遍历每个raw
  13.     for raw in raws:
  14.         # 生成频谱特征向量
  15.         feature_vector = []
  16.         # 遍历每个频段
  17.         for band in FREQ_BANDS:
  18.             # 提取每个频段的数据，不打印信息
  19.             raw_band = raw.copy().filter(l_freq=FREQ_BANDS[band][0], h_freq=FREQ_BANDS[band][1], verbose=False)
  20.             # 计算能量
  21.             power = np.sum(raw_band.get_data() ** 2, axis=1) / raw_band.n_times
  22.             # 添加到特征向量
  23.             feature_vector.extend(power)
  24.         # 添加到特征矩阵
  25.         feature_matrix.append(feature_vector)
  26.     # 返回特征矩阵
  27.     print("频谱特征矩阵的shape为：{}".format(np.array(feature_matrix).shape))
  28.     # print("频谱特征矩阵内容为：{}".format(np.array(feature_matrix)))
  29.     return np.array(feature_matrix)
  30.  
  31. raws, label = read_all_files("Preprocessed_EEG/", 1)
  32. frequency_spectrum(raws)

2. 功率谱密度

        我在实验中最常用的特征，计算每个频带的功率，反映在不同频率范围内的能量分布情况。其实这个特征和上面的频谱能量很类似，但是关于计算功率谱密度的函数，因为使用了加窗和fft方法，有几个新手经常疑惑的点需要注意，当时我查了大量资料才差不多搞明白。

        在计算功率谱的时候主要用到了mne中的一个函数compute_psd，这个函数的详细解释和加窗这块的知识在往期博文中可见，即mne库函数compute_psd()记录 。

        计算功率谱密度的结果如下：

read_all_files start...
共读取了1个文件
共有15段数据
read ended...
功率谱密度特征提取开始...
功率谱密度特征矩阵的shape为：(15, 310)
功率谱密度特征提取结束

        之前我计算功率谱的时候用的是比较短的epoch（不多于30s），而这次使用的数据段较长，不知道会有什么影响。两种情况的区别在于epochs是一个raw中的不同时间段，如果以时间为轴依次生成epoch，再对生成的epochs进行功率谱计算，能够得到一个以时间序列排序的功率谱数据，可以反映各个频段随时间的变化。但是以下代码仅对一整段raw计算：

  1. def power_spectrum(raws):
  2.     print("功率谱密度特征提取开始...")
  3.     # 计算每个频段的功率谱密度
  4.     # 特定频带
  5.     FREQ_BANDS = {"delta": [0.5, 4.5],
  6.                   "theta": [4.5, 8.5],
  7.                   "alpha": [8.5, 11.5],
  8.                   "sigma": [11.5, 15.5],
  9.                   "beta": [15.5, 30]}
  10.     # 特征矩阵
  11.     feature_matrix = []
  12.     # 遍历每个raw
  13.     for raw in raws:
  14.         # 生成频谱特征向量
  15.         feature_vector = []
  16.         # 遍历每个频段
  17.         for band in FREQ_BANDS:
  18.             # 计算功率谱密度
  19.             power = raw.compute_psd(picks='all', method='welch', fmin=FREQ_BANDS[band][0],
  20.                                            fmax=FREQ_BANDS[band][1], verbose=False)
  21.             # print(power.shape)
  22.             # 添加到特征向量，在第二个维度方向扩展
  23.             for i in range(power.shape[0]):
  24.                 feature_vector.extend(power[i])
  25.         # 添加到特征矩阵
  26.         # print(len(feature_vector))
  27.         feature_matrix.append(feature_vector)
  28.     # 将特征矩阵转换为numpy数组
  29.     feature_matrix = np.array(feature_matrix, dtype=object)
  30.     # 返回特征矩阵
  31.     print("功率谱密度特征矩阵的shape为：{}".format(feature_matrix.shape))
  32.     # print("功率谱密度特征矩阵内容为：{}".format(np.array(feature_matrix)))
  33.     print("功率谱密度特征提取结束")
  34.     return feature_matrix
  35.  
  36.  
  37. raws, label = read_all_files("Preprocessed_EEG/", 1)
  38. # frequency_spectrum(raws)
  39. power_spectrum(raws)

3. 连接性

        分析不同脑区域之间的相互关系，如相关性、相干性、互信息等。这部分的输出可以在图神经网络中应用。在mne的新版本中，与这块相关的函数已经另开了一个包，即mne_connectivity，需要另外install。这个包好像比较新，网上没有对应的内容，需要自己看文档。

        这里贴一下相干性和相关性的概念，其实我自己也不是很清楚这些：

        相干性(coherence )和相关性(correlation) 有什么区别和联系？

        下面代码使用了 mne_connectivity.spectral_connectivity_epochs 函数，按翻译来看计算的是频域和时频域的连接性。在设定只读取单个raw之后，使用该函数进行计算，得到了如下的输出：

共读取了1个文件
共有15段数据
read ended...
相干性特征提取开始...
event.shape:  (127, 3)
Using data from preloaded Raw for 127 events and 401 original time points ...
0 bad epochs dropped
epochs.shape:  (127, 62, 401)
con info:  <SpectralConnectivity | freq : [2.493766, 99.750623], , nave : 127, nodes, n_estimated : 62, 3844, ~5.8 MB>
相干性特征矩阵的shape为：(1, 3844, 196)
相干性特征提取结束

         以上可知，计算一段raw数据得到的相干性特征矩阵的维度是（3844，196），3844易知为62*62，即通道数的平方，而196则与设定计算的频域范围有关。

        更具体地说，频率分析在频率范围内进行，每个频率点会产生一组连接性值，而时间分析则是在一系列时间点上进行，因此这两个维度的组合构成了相干性特征矩阵的形状。

  • 频率维度：这是从频谱分析中获得的。假设你的频率范围是 [2.493766, 99.750623] Hz，你可能在这个频率范围内进行了若干个频率点的计算，每个频率点都计算出了一组相干性值。这些值被组织成了频率维度。

  • 时间维度：这是从你的数据时间序列中获得的。你可能在每个时间点上计算了一组相干性值，这些时间点的数量由你的原始数据的时间点数量决定。

        因此，196 表示了频率和时间点的组合，指的是在频率范围内的每个频率点上，你在一系列时间点上计算了相干性值，这些值被组织成了 196 维的时间维度。这个组合反映了你在频率和时间领域内进行的相干性分析，以便更好地理解信号的变化和交互。说实话不是很好理解，但是这个数据量较大，可能不太好处理，而且代码提取特征花费的时间也比较长久，如果能知道自己的任务需要用到哪个特定频域范围的话可以设置范围小一点的 fmin 和 fmax 参数来让计算速度快一点（提取一个文件中的15段raw数据计算特征就花了近一两个分钟左右）。

        代码如下：

  1. def spectral_connectivity(raws):
  2.     print("相干性特征提取开始...")
  3.     # 计算EEG数据的相干性特征
  4.     # 特征矩阵
  5.     feature_matrix = []
  6.     # 遍历每个raw
  7.     for raw in raws:
  8.         # 生成events, 2s为一个event，重叠0.2s
  9.         events = mne.make_fixed_length_events(raw, duration=2.0, overlap=0.2)
  10.         # print('event.shape: ', events.shape)
  11.         # 生成epochs
  12.         epochs = mne.Epochs(raw, events, tmin=0, tmax=2.0, baseline=None, verbose=False)
  13.         # print('epochs.shape: ', epochs.get_data().shape)
  14.         # 计算相干性矩阵
  15.         con = spectral_connectivity_epochs(epochs, method='coh', mode='multitaper', verbose=False)
  16.         # print('con info: ', con)
  17.         feature_matrix.append(con.get_data())
  18.     # 将特征矩阵转换为numpy数组
  19.     feature_matrix = np.array(feature_matrix, dtype=object)
  20.     # 返回特征矩阵
  21.     print("相干性特征矩阵的shape为：{}".format(feature_matrix.shape))
  22.     # print("相干性特征矩阵内容为：{}".format(np.array(feature_matrix)))
  23.     print("相干性特征提取结束")
  24.     return feature_matrix
  25.  
  26.  
  27. raws, label = read_all_files("Preprocessed_EEG/", 1)
  28. spectral_connectivity(raws)

 4. 时域特征

        这里提取的时域特征为均值、方差、标准差、最大最小值与差、斜率、偏度、峰度以及能量等，很简单。输出信息如下，其中15为raw个数，682 = 62（通道数） * 11（特征数）。

read_all_files start...
共读取了1个文件
共有15段数据
read ended...
时域特征提取开始...
时域特征矩阵的shape为：(15, 682)

        代码：

  1. def temporal_feature(raws):
  2.     # 提取EEG数据的时域特征
  3.     print("时域特征提取开始...")
  4.     # 特征矩阵
  5.     feature_matrix = []
  6.     # 遍历每个raw
  7.     for raw in raws:
  8.         # 生成时域特征向量
  9.         feature_vector = []
  10.         # 计算EEG数据的时域特征
  11.         # 1.均值
  12.         mean = raw.get_data().mean(axis=1)
  13.         # 2.方差
  14.         var = raw.get_data().var(axis=1)
  15.         # 3.标准差
  16.         std = raw.get_data().std(axis=1)
  17.         # 4.最大值
  18.         max = raw.get_data().max(axis=1)
  19.         # 5.最小值
  20.         min = raw.get_data().min(axis=1)
  21.         # 6.峰峰值
  22.         p_p = max - min
  23.         # 7.斜率
  24.         slope = np.diff(raw.get_data(), axis=1)
  25.         # 8.峭度
  26.         kurtosis = np.mean((raw.get_data() - mean[:, None]) ** 4, axis=1) / var ** 2 - 3
  27.         # 9.偏度
  28.         skewness = np.mean((raw.get_data() - mean[:, None]) ** 3, axis=1) / var ** (3 / 2)
  29.         # 10.能量
  30.         power = np.sum(raw.get_data() ** 2, axis=1) / raw.n_times
  31.         # 11.方根幅值
  32.         rms = np.sqrt(np.mean(raw.get_data() ** 2, axis=1))
  33.         # 12.脑电活动时域特征向量
  34.         feature_vector.extend(mean)
  35.         feature_vector.extend(var)
  36.         feature_vector.extend(std)
  37.         feature_vector.extend(max)
  38.         feature_vector.extend(min)
  39.         feature_vector.extend(p_p)
  40.         feature_vector.extend(slope)
  41.         feature_vector.extend(kurtosis)
  42.         feature_vector.extend(skewness)
  43.         feature_vector.extend(power)
  44.         feature_vector.extend(rms)
  45.         # 添加到特征矩阵
  46.         feature_matrix.append(feature_vector)
  47.     # 将特征矩阵转换为numpy数组
  48.     feature_matrix = np.array(feature_matrix, dtype=object)
  49.     # 返回特征矩阵
  50.     print("时域特征矩阵的shape为：{}".format(feature_matrix.shape))
  51.     # print("时域特征矩阵内容为：{}".format(np.array(feature_matrix)))
  52.     return feature_matrix
'
运行

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/木道寻08/article/detail/890656
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号