Hibiscus的脑机接口学习周报（2022/12/26~2022/1/1）

2022/12/26（星期一）

在学习使用python mne库读取.set类型数据时，遇到了类似如下报错：buffer is too small

（来源于Python读取.set格式的脑电数据出现buffer is too small的问题）

博客内提供可行方法为：使用eeglab将set格式的数据转化为edf或者mat类型数据再进行读取。然而再次遇到问题：在下载plugin插件biosig之后，set文件类型似乎无法正确转化为edf类。

from mne.io import concatenate_raws, read_raw_edf
import matplotlib.pyplot as plt
import mneraw=read_raw_edf("third.edf", preload=False)events_from_annot, event_dict = mne.events_from_annotations(raw)
print(event_dict)
print(events_from_annot)

可以将数据读入，但输出似乎存在异常，如图，输出为空字段。

C:\Users\86136\anaconda3\envs\bci\python.exe C:/Users/86136/Desktop/innovation/python_mne/eeglab-read-set/import-edf.py
Extracting EDF parameters from C:\Users\86136\Desktop\innovation\python_mne\eeglab-read-set\third.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
{}
[]

在尝试过后得到了解决原始问题的方法（为何不能读入set类型数据）

似乎在被eeglab处理过后的数据就不能读入mne中，在使用eeglab自带数据集eeglab-data.set时程序运行正常。eeglab将set转化为edf数据时的异常留待解决，可能由插件版本问题造成。

202212/27（星期二）

进行之前一些操作的整理

主要着重于raw、epoch以及evoked三种数据结构以及之间的转化关系：

python_mne数据结构

对信号空间投影相关概念以及注释做了学习。

信号空间投影SSP

一旦知道了噪声向量，就可以创建一个与其正交的超平面，并构造一个投影矩阵，将实验记录投影到该超平面上。这样，测量中与环境噪声相关的部分就可以被移除。同样，应该清楚的是，投影降低了数据的维数-你仍然会有相同数量的传感器信号，但它们不会都是线性独立的-但通常有数十或数百个传感器，而你要消除的噪声子空间只有3-5维，因此自由度的损失通常是没有问题的。

注释连续数据(Annotating continuous data)

MNE-Python中的注释是一种存储关于原始对象的时间跨度的短字符串信息的方法。
注释是类似列表的对象，其中每个元素包含三部分信息:起始时间(以秒为单位)、持续时间(也以秒为单位)和描述(文本字符串)。
此外，annotation对象本身也跟踪orig_time，它是一个POSIX时间戳，表示相对于注释开始的真实时间。

2022/12/28（星期三）

对大部分样例进行实操。

对空间信号投影（SSP）进行了实操

在进行注释类型数据操作的时候出现了报错

执行以下代码：

# 构建事件戳
print(raw.info['meas_date'])
meas_date = raw.info['meas_date'][0] + raw.info['meas_date'][1] / 1e6
orig_time = raw.annotations.orig_time
print(meas_date == orig_time)

报错显示meas_date不能使用下标的方式访问？

Traceback (most recent call last):
File "C:\Users\86136\Desktop\innovation\python_mne\Annotation.py", line 21, in <module>
meas_date = raw.info['meas_date'][0] + raw.info['meas_date'][1] / 1e6
TypeError: 'datetime.datetime' object is not subscriptable

查询文档可知meas_date属于日期时间对象或元组类型，上述采用元组类型的提取方法（似乎）失效，故查询日期时间对象的提取方法。

如果需要小批量对数据进行注释，可视化标注似乎是一个不错的方法呢。

2022/12/29（星期四）

对机器学习算法随机森林判断睡眠类型进行进一步学习，主要掌握如何使用机器学习算法对eeg信号进行处理。

在学习的过程中发现自己对python的部分基础语法和库的用法极其不熟悉，在mne与plt结合实现可视化的过程中遇到了较大的阻碍，于是查阅了较多的资料。

plt: subplot()、subplots()详解及返回对象figure、axes的理解

matplotlib 是从MATLAB启发而创建的，这俩的命令几乎一模一样，所以概念基本完全相同，既然python教程看的那么麻烦，何不看看MATLAB官方文档，给你解释的明明白白的还附带一堆例子。

在绘图代码中这一行最令人费解，在查询过后没有找到让我满意的答案。

plt科研绘图 plt科研绘图，也许这一篇博客能有一些参考价值。
# 保留颜色代码以便进一步绘制
stage_colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
在进一步查找过程中发现这可能是返回了一个颜色的序列，是一个带有默认颜色参数的列表（？），在睡眠类型特征提取的绘图部分中，stage_colors这样被调用:
for ax, title, epochs in zip([ax1, ax2],['Alice', 'Bob'],[epochs_train, epochs_test]):for stage, color in zip(stages, stage_colors):epochs[stage].plot_psd(area_mode=None, color=color, ax=ax,fmin=0.1, fmax=20., show=False,average=True, spatial_colors=False)
这段有点绕的循环画出了各个event在各自epoch中的波形，在可视化后表现非常直观。

在频率谱psd分析方面，参考学习了mne中psd分析这一篇博客，里面的链接还包括了功率谱相关的详解，对0基础的人很友好。

PSD的意义就在于将不同频率分辨率下的数据归一化，排除了分辨率的影响，得到的PSD曲线趋势是一致的

在实现以下代码的时候发现弹出警告：

[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.8s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.8s finished
<decorator-gen-201>:4: FutureWarning: NOTE: psd_welch() is a deprecated function. Function ``psd_welch()`` is deprecated; for Raw/Epochs/Evoked instances use ``spectrum = instance.compute_psd(method="welch")`` instead, followed by ``spectrum.get_data(return_freqs=True)``.

这段报错指向这一段代码：

psds, freqs = psd_welch(epochs, picks='eeg', fmin=0.5, fmax=30.)

意思是mne中的函数psd_welch已经被弃用了，所以我们需要使用别的方法进行替代。

在mne文档的更新部分查找，我们能看到这样的提示：

The PSD functions that operate on Raw/Epochs/Evoked instances (mne.time_frequency.psd_welch and mne.time_frequency.psd_multitaper) are deprecated; for equivalent functionality create Spectrum or EpochsSpectrum objects instead and then run spectrum.get_data(return_freqs=True) (#10184 by Daniel McCloy)

psd_welch和用来计算能量和的psa_multitaper函数都已经被淘汰了，需要创造Spectrum或者EpochsSpectrum对象，并且使用spectrum.get_data(return_freqs=True)函数作为替代。

在csdn上好像找不到相关的操作。好吧，又要重新看mne文档了。

两个可选对象中，前者Spectrum与evoked有关，但是我已经懒得再转化成evoked了，操作也并不是很熟练，所以直接选用了后者EpochsSpectrum。

在它的method那一栏看到了熟悉的welch，说明离目标更近了一步。但是上面出现了一行提醒：

The preferred means of creating Spectrum objects from Epochs is via the instance method mne.Epochs.compute_psd(). Direct class instantiation is not supported.

好像直接的类实例化是不支持的，对于epoch数据结构的首选方式是通过compute_psd()函数来获取光谱对象。好吧，那去查找一下这个函数。查询得compute_psd()返回的是每一个epoch的光谱表示，而且刚好是EpochsSpectrum的频谱实例，那正好。

经过很多繁琐的查找操作后把上面的代码改写为以下几行，成功实现了同样的效果，并且程序没有报任何错误：

    # psds, freqs = psd_welch(epochs, picks='eeg', fmin=0.5, fmax=30.)spectrum = epochs.compute_psd(method='welch', picks='eeg', fmin=0.5, fmax=30.)psds = spectrum.get_data()freqs = spectrum.freqs

自己查询文档写代码的成就感还是挺大的，今天算是没有白学。

晚上继续看了接下来的一部分内容。

今天的学习有很多收获，同时也暴露了很多的不足之处。

2022/12/30（星期五）

中午学习evoked可视化部分。

（下午一点开始选课，选完已经快两点了，感觉今天没啥学的状态了，学一点是一点吧）

evoked诱发电位在我的印象里是由特定事和通道的epochs平均化得到的事件刺激特征。“原来在这种刺激下会出现这种反映”，就是这样的感觉。evoked相当于是对一个事件刺激的反应总结。

在print之后可以看见数据集提供了很完整的基于四个event的evoked数据:

Left Auditory、Right Auditory、Left visual、Right visual。

[<Evoked | 'Left Auditory' (average, N=55), -0.1998 – 0.49949 sec, baseline -0.199795 – 0 sec, 376 ch, ~4.5 MB>, <Evoked | 'Right Auditory' (average, N=61), -0.1998 – 0.49949 sec, baseline -0.199795 – 0 sec, 376 ch, ~4.5 MB>, <Evoked | 'Left visual' (average, N=67), -0.1998 – 0.49949 sec, baseline -0.199795 – 0 sec, 376 ch, ~4.5 MB>, <Evoked | 'Right visual' (average, N=58), -0.1998 – 0.49949 sec, baseline -0.199795 – 0 sec, 376 ch, ~4.5 MB>]

接下来是evoked.plot()环节，注意到与教程不一样的是似乎不用自己另外设置参数，视窗中显示的图片就已经是彩色的了。

感觉如下的subplots绘图方法还是有必要熟练掌握的：

fig, ax = plt.subplots(1, 5, figsize=(8, 2))
kwargs = dict(times=0.1, show=False, vmin=-300, vmax=300, time_unit='s')
evoked_l_aud.plot_topomap(axes=ax[0], colorbar=True, **kwargs)
evoked_r_aud.plot_topomap(axes=ax[1], colorbar=False, **kwargs)
evoked_l_vis.plot_topomap(axes=ax[2], colorbar=False, **kwargs)
evoked_r_vis.plot_topomap(axes=ax[3], colorbar=False, **kwargs)
for ax, title in zip(ax[:4], ['Aud/L', 'Aud/R', 'Vis/L', 'Vis/R']):ax.set_title(title)
plt.show()

evoked.plot_joint()函数可以结合头皮图和频谱图，看起来也更加直观。

对于移动高级绘图使用mne.viz.plot_compare_evokeds(),这个方法的使用在我看来还是有些难以理解，尤其是这里的pick：

pick = evoked_dict["Left/Auditory"].ch_names.index('MEG 1811')

教程中的pick是从字典中选取的通道，指定在左听觉部分选取通道MEG1811，但显示的图片却包括了所有的事件，一开始不免令人疑惑。后来查询文档然后实践过后可得以下代码是等价的：

# pick = evoked_dict["Left/Auditory"].ch_names.index('MEG 1811')
pick = ['MEG 1811']

重点在于通道的名称，而不是通道是在哪里选取的。

按照这个进度，这个星期应该就可以初步学完基础的部分

以下是几个未解决的问题，留给未来的自己吧

1.evoked专属的数据集（_ave.fif）应该怎么创建

2.evoked.plot()之后显示的图该如何进行具体分析，我们能从图中得到什么有用的信息？

3.磁力计和磁强在脑机接口中的意义

2022/12/31（星期六）

今天学习脑电图处理和ERP（事件相关电位）有关内容，其实感觉就是前面部分的总结，主要重在复习方面。

EEG（Electroencephalography，脑电图）是指脑神经细胞无时无刻不在进行自发性、节律性、综合性的电活动，将这种电活动的电位作为纵轴，时间特征作为横轴，记录下来的电位与时间相互关系的平面图即为脑电图。

ERP（Event Related Potential，事件相关电位）是指外加一种特定的刺激作用于感觉系统或脑的某一部位，给予刺激或撤销刺激时，在脑区所引起的电位变化。

那么EEG与ERP之间的关系如何？从概念中可知，EEG属于大脑的自发性电活动，主要表现为在广泛频谱上占主导地位且具有某些特征的波形。

ERP是指当人接收到与特定感觉、认知或运动事件相关的刺激时，自发性EEG活动会受到干扰。这种由事件诱发的神经响应会淹没在自发性EEG活动中，但可以借由简单的平均叠加技术或更复杂的单试次分析和时频分析等技术，将这些响应从自发性EEG活动中提取出来。这些经过平均叠加所获得的脑电响应被称为事件相关电位，表示它们是与特定事件相关的电位。

ERP与EEG思维导图

来自这位博主的思维导图看起来很直观

这里简要记录一下教程中对脑电图进行操作的处理步骤：

1.加载数据文件，以raw格式记录。

2.使用pick_types()筛选数据

3.使用mne.io.Raw.plot_sensors()绘制通道位置，数据没有位置，则可以使用MNE随附的Montages来设置

4.从原始对象中删除参考

5.定义Epochs，并计算左听觉状态的ERP

6.设置平均参考电极/自定义参考

7.创建一个包含4个条件的Epochs对象

8.创建左刺激和右刺激试验的平均值（evoked）

9.构建和绘制不同的ERP

tips：在字典中或者list中储存不同事件的evoked对象有利于检索

今天接下来的时间用于学习协方差矩阵相关方面。

首先提出一个问题：

协方差矩阵是什么？有什么意义？

协方差矩阵相关概念、性质、应用意义及矩阵特征向量的用处

在大致浏览完教程之后，感觉还是有许多困惑的地方没有解决，因为教程基本上通过文档机翻润色而来，很多地方还是感觉颇为生硬。

所以先把这个放一下，暂且回去复习以下机器学习判断睡眠类型的内容。

在回顾的时候发现自己对频谱分析这块的概念还不是很了解，写代码也就是照葫芦画瓢罢了，缺少对现象本质的理解。

在查找资料过程中发现了这个挺好的教程，把时频域分析讲的很清楚：脑电信号时频域分析

现在想来睡眠分类这块的函数mne.time_frequency.Spectrum()不就是时频域谱的意思吗...

感觉自己还是个半吊子

这个教程讲的很详细，看完感觉自己的收获很大，很多之前自己一知半解的问题都得到了解答，比如为什么要将波形分解为αβ等几个频段进行分析，傅里叶变换的意义等。之前导师好像都在上课的时候讲过类似的问题，但是没怎么记住。想来光是听课也学不到什么有用的东西，只有在课后实践主动想要了解的知识才能被自己记住。

现在重新来看这段代码：

def eeg_power_band(epochs):"""脑电相对功率带特征提取该函数接受一个""mne.Epochs"对象，并基于与scikit-learn兼容的特定频带中的相对功率创建EEG特征。Parameters----------epochs : EpochsThe data.Returns-------X : numpy array of shape [n_samples, 5]Transformed data."""# 特定频带FREQ_BANDS = {"delta": [0.5, 4.5],"theta": [4.5, 8.5],"alpha": [8.5, 11.5],"sigma": [11.5, 15.5],"beta": [15.5, 30]}# psds, freqs = psd_welch(epochs, picks='eeg', fmin=0.5, fmax=30.)spectrum = epochs.compute_psd(method='welch', picks='eeg', fmin=0.5, fmax=30.)psds = spectrum.get_data()freqs = spectrum.freqs# 归一化 PSDspsds /= np.sum(psds, axis=-1, keepdims=True)X = []for fmin, fmax in FREQ_BANDS.values():psds_band = psds[:, :, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)X.append(psds_band.reshape(len(psds), -1))return np.concatenate(X, axis=1)

在这个函数中：

spectrum = epochs.compute_psd(method='welch', picks='eeg', fmin=0.5, fmax=30.)
psds = spectrum.get_data()
freqs = spectrum.freqs

这一段代码的意义在之前便已经探讨过，是用welch法对epochs数据进行频谱分析，并且选取eeg通道的数据，截取波段0.5s到2s的长度。

对于频谱分析采用的方法异同，这篇文章中有较为详细的解释：eeg频谱分析方法，当然前一个讲时频域分析的博客讲的更好一些，值得多看几遍,其中提到welch方法的主要优点是能得到较为平滑的频谱曲线。

而compute_psd这个函数所返回的是spectrum数据结构，即“数据的频谱表示”（不知道为什么老是被网页翻译成光谱），也就是把时域表示的信号转化为频谱格式，具体如图：

在get_data()中psds得到了numpy数组格式的频谱数据，freqs得到的是一个频率数组。

最后重点在这儿一段：

for fmin, fmax in FREQ_BANDS.values():psds_band = psds[:, :, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)X.append(psds_band.reshape(len(psds), -1))

psds_band读取的到底是什么？为什么要这样提取psds中的数据？

这取决于psds的数组结构到底是什么样的，频谱数据是以什么形态存在里面。文档里似乎也没有给出说明，就是这里看的我迷迷糊糊的，然后我print了psds，似乎看不出什么端倪。

（2022/1/1 ps:之后查阅文档发现在后面的例子操作里面给出了说明(2802, 2, 75)，2802指epoch数量，2为通道数量（eog、misc），75为频率点数量）

this is psds
(2802, 2, 75)
[[[2.69509693e-10 1.36308256e-10 1.72018742e-10 ... 1.13069360e-131.20029249e-13 1.22991691e-13][1.34306665e-11 6.11940888e-12 3.24737467e-12 ... 3.55286266e-131.49511584e-13 2.10815700e-13]][[6.72512425e-10 6.53372564e-10 2.62354442e-10 ... 2.38603381e-132.34371591e-13 2.67583329e-13][2.27644810e-11 1.87908513e-11 5.12818638e-12 ... 3.98176515e-134.98784339e-13 6.07921963e-13]][[3.67234583e-10 2.76748578e-10 1.39229012e-10 ... 1.20067444e-131.28366221e-13 1.54891622e-13][1.76619334e-11 1.70534925e-11 2.88135066e-12 ... 1.29066707e-132.63699801e-13 1.95785116e-13]]...[[7.41855993e-10 2.10180223e-10 7.93061474e-11 ... 1.18051585e-131.89930521e-13 1.73559386e-13][1.60005625e-11 1.28544385e-11 3.40883735e-12 ... 2.34571839e-125.55404011e-12 4.28265328e-12]][[2.79403818e-10 1.22617958e-10 4.33599555e-11 ... 1.89819586e-132.13947811e-13 1.57405119e-13][1.61491040e-11 4.62766959e-12 3.50591136e-12 ... 3.45079821e-125.79890837e-12 4.63585865e-12]][[3.78707842e-10 2.23923702e-10 1.64358990e-10 ... 1.92693116e-131.51643092e-13 6.72470066e-14][2.31340631e-11 1.26137243e-11 7.73311628e-12 ... 5.38547188e-124.45856052e-12 2.57746185e-12]]]

但是直接从代码出发看这个结构，那psds中存储的就是所有的频率。

yysy找不到头绪还挺绝望的，不过特征提取部分的教程给了我答案：既然函数的目的是提取数据特征，那只要知道所谓的特征是什么，是如何表示的，就能够知道上面的计算原理！

于是我们得到了这样的答案，for循环是通过遍历来平均化区间内频率点，用来计算几个频率范围内最重要的频率点，即带限功率。每一轮计算过后，都将相应的结果存储在X中。

数组布尔索引高级

psds_band = psds[:, :, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)

至于这一句代码的细节，则在于此：

psds是经过归一化的频谱数据，即功率，而freqs可以说是未经处理的频率数据

在这一句代码中，每一个功率点所在的频率由freqs来表示，如果该功率点的所在频率满足当前的fmin与fmax条件，则将该功率点的功率加入功率和中，最后取平均。

至psds的shape.......我不是很懂为什么是长这样的，所以以上的猜测除了理论之外都可能是错误的。。。。。。

今天就到此为止了，明天再看吧。

明天就是元旦了，现在是晚上21:51分，想来也没有朋友会私发给我新年快乐，距离那个年头已经很久啦。那时候大家什么也不在乎，只需要朝着一个目标努力就好，全世界的人好像都在为自己加油。

提前祝自己新年快乐。

2022/1/1（星期天）

元旦快乐

昨天太累了就没有怎么看，今天起来print了freqs的shape

this is freqs shape
(75,)

是75，那就意味着昨天的猜想没错，0.5hz到30hz被分为了75个频率点，使psds的功率与freqs的频率一一对应，接下来进一步验证这75个点是怎么分出来的

this is freqs
[ 0.78125   1.171875  1.5625    1.953125  2.34375   2.734375  3.1253.515625  3.90625   4.296875  4.6875    5.078125  5.46875   5.8593756.25      6.640625  7.03125   7.421875  7.8125    8.203125  8.593758.984375  9.375     9.765625 10.15625  10.546875 10.9375   11.32812511.71875  12.109375 12.5      12.890625 13.28125  13.671875 14.062514.453125 14.84375  15.234375 15.625    16.015625 16.40625  16.79687517.1875   17.578125 17.96875  18.359375 18.75     19.140625 19.5312519.921875 20.3125   20.703125 21.09375  21.484375 21.875    22.26562522.65625  23.046875 23.4375   23.828125 24.21875  24.609375 25.25.390625 25.78125  26.171875 26.5625   26.953125 27.34375  27.73437528.125    28.515625 28.90625  29.296875 29.6875  ]

虽然不是很清楚这样分布的原理，但是基本每频率点前后差了大约0.4hz的距离。

根据这一段对频谱分析的描述可以看出频率点为什么会这样分布！

频率点的间隔为Fs/M

Fs显然是指0.5hz到30hz的间隔，而数据段的长度M是指什么？

在之前的代码中有一段其实让我很困惑，而我们刚学习的知识能为我们很好地解答这个问题：

tmax = 30. - 1. / raw_train.info['sfreq']  # tmax in included
"""
所创建的是时间从tmin=0开始，到tmax为止的epochs
"""

为什么之前取的tmax需要再减去 1/sfreqs ，为什么是这个数字？

根据频率点间隔=Fs/M 这个公式，我们能带入已知的频率点间隔0.390625hz以及Fs算出M等于75，即当前频率长度，则数据段长度M为75，即welch方法的windows大小为75个采样点，然而原数据的采样频率为100hz，epoch中的时间为0~29.99s

this is times here
[0.000000e+00 1.000000e-02 2.000000e-02 ... 7.949997e+04 7.949998e+047.949999e+04]
this is sfreq here
100.0

那么可以得知epoch中存在3000个时间采样点。

3000/75=400

即数据点在0~29.99s时，刚好可以被M被分为400个数据段，每段数据段存在75个数据点。

在compute_psd()使用method welch时，默认时间窗口重叠部分为0。

然而这个75其实是我们从数据倒推出来的，实际上为什么取的是75呢？

或者说为什么epoch要被分为400段呢,这个数字是怎么计算出来的？

查了好久资料还是一无所获，晚上再说吧......

晚上对这个的研究终于有了突破，另写了一篇稍微总结了一下compute_psd()的计算思路，今天算是这一周的圆满结束了:

脑电分析mne库函数compute_psd()记录