numpy微实践 - 生物电信号的叠加平均

海洋饼干叔叔

发布于：2024年3月23日

微实践

人体是一个复杂的带电体，大脑与人体其他部分之间依赖于生物电信号协同工作，肌肉有肌电信号、心脏有心电信号、大脑有脑电波。记录并提取这些与人体运作有关的电信号可以帮助医生诊断疾病。

版权声明

本文可以在互联网上自由转载，但必须：注明出处(作者：海洋饼干叔叔)并包含指向本页面的链接。

本文不可以以纸质出版为目的进行改编、摘抄。

本文节选自《Python编程基础及应用》，高等教育出版社，作者：陈波，刘慧君

图1 视觉诱发电位记录的工作原理示意

图1展示了一种名为“视觉诱发电位（Visual Evoked Potential，简称VEP）”的人体生物电信号的记录过程。位于人眼前的图形刺激器产生棋盘格图像对人眼进行“刺激”，通过角膜、晶状体、房水等构成的屈光系统，棋盘格在视网膜上成像。视网膜上的感光细胞进行光电反应，产生出电信号。然后，这些源自视网膜感光细胞的电信号通过视神经，最终到达大脑的视觉皮层，引起视觉皮层的兴奋，经大脑的运算推理，人才得以“看”到了棋盘格影像。

为了记录这种因为视觉刺激而诱导出来的视觉皮层电信号，如图1所示，我们将生物放大器的正电极安装在视觉皮层所在位置的颅骨外侧来采集这些电信号，然后通过模数转换，传入计算机，就会得到随时间t变化的电压值曲线。

遗憾的是，生物电放大器采集得到的信号并非完全来自视觉皮层，可能还包括来自于心脏、胃、肌肉、听觉系统、大脑中枢的各种电信号噪声。为了提高信噪比，需要使用叠加平均的技术，即以约500ms为周期，使用棋盘格图像多次重复刺激人眼，在记录到随时间t变化的原始完整曲线后，再将各次刺激相关的信号段相加，再平均，从而将与刺激有因果和时间关系的“视觉诱发电位”信号保留下来，而把那些与刺激无因果和时间关系的噪声去除。

在随书代码的CH10目录，vep_online.dat存储了完整的原始信号，下述程序用于从文件中读取原始信号序列至数组：

#addaverage.py
import json,numpy as np,matplotlib.pyplot as plt

samples = np.fromfile("vep_online.dat",dtype=np.float32)
print("Shape of samples:",samples.shape)

上述程序片段第5行对应的输出结果为：

1	Shape of samples: (65572,)

np.fromfile()函数用于从文件直接数据至numpy数组。由于关键字参数dtype指定元素类型为np.float32（32位浮点数，占4个字节），函数将会对vep_online.dat的文件内容进行拆分，每4个字节转换成1个浮点数。程序第5行的输出结果证实，结果数组包含65572个样本（sample）浮点数，每个样本表示一个时间轴上的电压值，对应4个字节，如果去操作系统中查看vep_online.dat文件的尺寸，一定是65572 × 4 = 262288字节。

文件vep_stimuli.json以json格式存储了一个整数列表，包含了每次棋盘格刺激的发生位置，对应samples样本序列中的下标。下述程序片段从文件中写出该列表，并切片打印了前6个值。json文件读写的具体细节，请回顾本书8.4节，此处不再赘述。

1
2
3

with open("vep_stimuli.json","rt") as f:
    stiPoses = json.load(f)
print("stimuli positions:",stiPoses[:6])

上述程序片段print()对应的输出结果为：

1	stimuli positions: [1545, 2589, 3591, 4594, 5637, 6641]

上述输出表明，在样本序列的下标1545，2589等处，对人眼施加了图形刺激。在接下来的下述程序片段中，程序遍历各次刺激的位置下标，叠加平均出结果信号。

signal = np.zeros(600)         #创建0值数组用于保存结果信号
plt.figure(figsize=(10,5))     #创建10英寸宽，5英寸高的图
ax = plt.subplot(121)          #在图中的左侧创建子图
ax.set_title("Raw VEP Signal") #左子图标题为"原始VEP信号"
for p in stiPoses:             #遍历全部刺激位置下标
    t = samples[p:p+600]       #从刺激位置开始切片600个样本
    signal += t                #累加样本片段至结果数组
    ax.plot(t)                 #在左子图绘出样本片段
signal /= len(stiPoses)        #对累加的结果数组进行平均

🚩第1行：创建包括600个0值的数组signal，用于存储结果信号。原始信号序列对应的采样频率为2000Hz，即每秒钟放大器以相等时间间隔测量正负极电压差2000次，两个相邻样本的时间间隔为0.5ms，600个样本的对应时间跨度为300ms。

🚩第5行：遍历stiPoses列表，其中包含各次刺激在原始信号序列中的对应下标p。

🚩第6行：以下标p从原始信号序列samples中切片，得与刺激相关的信号段序列t，其长度也是600个样本。

🚩第7行：将与“本次刺激”相关的信号段序列t累加至结果信号数组signal。加法操作所涉及的两个数组均包含600个浮点数，形状相同，加法操作事实上通过ufunc函数完成，其实际动作是将两个数组的对应下标元素相加，得到另一个形状相同的新数组。

🚩第8行：将与“本次刺激”相关的信号段序列t绘制在左子图中。ax.plot(t)仅为plot()函数提供了一个实参t，函数在绘图过程中会将t中元素的下标视为x坐标，而将t中元素的值视为y坐标。

🚩第9行：在累加循环结束后，将结果数组signal中的全部元素值除以刺激次数（len(stiPoses)），得到单次刺激对应的结果信号的均值。

ax = plt.subplot(122)          #在图中创建右子图
ax.plot(signal)                #在右子图中绘制结果VEP信号
ax.set_title("Add Adveraged VEP Signal") #右子图标题为"叠加平均后的VEP信号"
plt.subplots_adjust(0.05,0.05,0.95,0.95,hspace=0.1)
plt.show()

🚩第2行：将叠加平均后的结果信号signal绘制在右子图中。

🚩第4行：调整图的边距以及子图间距。关于subplots_adjust()函数中相关参数含义，请回顾本书8.3节。

图2 视觉诱发电位记录的叠加平均

图2展示了前述程序的绘图结果。左子图中那些包含噪声信号的杂乱无章的与刺激相关的信号段，在经过叠加平均后，露出了纯净的视觉诱发电位(VEP)信号的真容。从右子图可见，VEP信号在样本序号200处呈现一个明显的峰，在样本序号150和300处呈现两个明显的谷，识别和测量这些信号成份的时间和振幅可以帮助医生诊断视觉系统疾病。一个典型的应用场景与法医鉴定有关，在一个故意假装眼盲的“病人”身上如果可以诱发出可信的VEP信号，即可证明他事实上能看见。

示例完整程序：

#addaverage.py
import json,numpy as np,matplotlib.pyplot as plt

samples = np.fromfile("vep_online.dat",dtype=np.float32)
print("Shape of samples:",samples.shape)

with open("vep_stimuli.json","rt") as f:
    stiPoses = json.load(f)
print("stimuli positions:",stiPoses[:6])

signal = np.zeros(600)         #创建0值数组用于保存结果信号
plt.figure(figsize=(10,5))     #创建10英寸宽，5英寸高的图
ax = plt.subplot(121)          #在图中的左侧创建子图
ax.set_title("Raw VEP Signal") #左子图标题为"原始VEP信号"
for p in stiPoses:             #遍历全部刺激位置下标
    t = samples[p:p+600]       #从刺激位置开始切片600个样本
    signal += t                #累加样本片段至结果数组
    ax.plot(t)                 #在左子图绘出样本片段
signal /= len(stiPoses)        #对累加的结果数组进行平均

ax = plt.subplot(122)          #在图中创建右子图
ax.plot(signal)                #在右子图中绘制结果VEP信号
ax.set_title("Add Adveraged VEP Signal") #右子图标题为"叠加平均后的VEP信号"
plt.subplots_adjust(0.05,0.05,0.95,0.95,hspace=0.1)
plt.show()