功能性近红外光谱技术在脑科学领域被广泛应用，市面上也已经有了许多基于MATLAB的优秀工具包及相关教程，如：homer、nirs_spm等。而本次教程将基于Python的MNE库对fNIRS数据进行处理。
本次教程基于：https://mne.tools/stable/auto_tutorials/preprocessing/70_fnirs_processing.html
教程所用代码、数据可通过链接下载，也可在茗创科技公众号后台回复关键词：MNE-FNIRS获取本次教程所用的工具、代码、数据。
课前准备
Anaconda下载安装（windows，64位系统）
主链接：https://repo.anaconda.com/archive/Anaconda3-2021.11-Windows-x86_64.exe
备用链接（该链接下载限速1000kb/s，仅在主链接无法访问时使用）：
http://1.14.108.54:8888/#s/7t2FaqPA
请按照默认选项安装，需要注意：
一定要勾选以下选项：
（如果忘记勾选，卸载软件后重新安装即可）
安装完毕后在开始菜单（左下角的windows图表
）找到
和
，两个都可以顺利打开证明安装成功。
打开Spyder
代码讲解实操
对于有Python基础的同学此步并不复杂，但此次课程面向零基础小白，将会逐步运行代码及讲解代码含义。
（1）载入软件库及fnirs数据
运行此次脚本前应加载好所需模块。
import numpy as np如运行以上语句，左下角提示无某模块时（类似下图），
则通过anaconda powershell prompt载入相应模块。
可通过conda或pip或wheel语句载入模块。
成功后再运行代码无报错提示。
如果你先前未下载示例数据，可通过代码下载，运行代码：
fnirs_data_folder = mne.datasets.fnirs_motor.data_path()
右下角Python控制台将出现一个进度条下载数据，如下图为下载完成：
如果没有成功下载，出现提示某某路径下不存在mne_data文件夹，则在相应路径（一般是在C盘相应用户文件夹下）创建mne_data文件夹即可。此时打开：
Python提示路径下mne_data文件夹，将会看到本次脚本所用示例数据。
运行此段落剩下三句，将会指定数据所在路径、读取数据、载入数据。
fnirs_cw_amplitude_dir = fnirs_data_folder / 'Participant-1'
右上角变量区将会出现相应变量，如下图所示。
如果已经下载好数据，想通过路径读取，或者想读取自己的fnirs数据，可参考链接：https://mne.tools/stable/auto_tutorials/io/30_reading_fnirs_data.html
如果读取NIRX数据，则使用函数：
mne.io.read_raw_nirx(fname,saturated='annotate', preload=False, verbose=None)如果读取Hitachi日立设备数据，则使用函数：
mne.io.read_raw_hitachi(fname,preload=False, verbose=None)如果读取SNIRF (.snirf)数据（注：.nirs格式数据可通过homer3转换成.snirf格式数据），则使用函数：
mne.io.read_raw_snirf(fname,optode_frame='unknown', preload=False, verbose=None)fname为数据所在路径，因此上述读取数据的代码可改为（数据所在路径应根据自身情况修改）：
raw_intensity=mne.io.read_raw_nirx('C:/Users/Administrator.DESKTOP-J86OEI2/mne_data/MNE-fNIRS-motor-data/Participant-1',saturated = 'annotate' , preload = False , verbose = None)运行后右上角同样出现变量
（2）指定duration，及修改mark
首先解释一下本次示例数据的实验mark，数据中共有四个mark，其中15.0为实验开始/结束mark，1.0为控制条件mark，2.0为Tapping/左侧条件mark，3.0为Tapping/右侧条件mark。每个刺激mark持续时间为5s。
代码中使用如下语句进行指定duration、修改mark、删除无用mark操作。
raw_intensity.annotations.set_durations(5)
运行raw_intensity.annotations.set_durations(5)语句后，数据的所有刺激持续时间都变成5秒。我们可以在变量区打开raw_intensity-->annotations-->set_durations，查看该语句的使用方法。
当set_durations后括号内只有一个数字时，所有刺激的持续时间都将改为该数字。
如果需要将不同刺激修改成不同持续时间，则参考示例语句：{'ShortStimulus' : 3, 'LongStimulus' : 12} 。该语句意为将Mark名为'ShortStimulus'的刺激的持续时间改为3，Mark名为'LongStimulus'的刺激的持续时间改为12。
raw_intensity.annotations.rename语句作用为修改Mark名字，示例语句为将Mark名为1的Mark改名为control
unwanted = np.nonzero(raw_intensity.annotations.description == '15.0')和raw_intensity.annotations.delete(unwanted)两句为删掉Mark名为‘15’的无用Mark。在本例数据中15代表实验开始和结束，与分析无关，故舍弃。
（3）删除短通道
研究者认为短通道（光电二极管之间的距离小于1厘米）无法有效检测神经反应，因此保留其他有效通道（非短通道）。
picks = mne.pick_types(raw_intensity.info, meg=False, fnirs=True)
运行完raw_intensity.pick(picks[dists > 0.01])后点开raw_intensity的ch_name变量，可见通道数变为40，删去了距离小于1cm的16个通道。
运行以下语句，可看到被保留的通道及数据整体情况。
raw_intensity.plot(
随后将原始光强度转换为光密度。
raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
变量区增加raw_od为光密度信息。
（此图为光密度信息）
（4）SCI法评估数据质量
研究者使用使用头皮耦合指数（scalp coupling index）来量化头皮与光电器件之间的耦合质量。
（此方法参考文献为：Luca Pollonini, Cristen Olds, Homer Abaya, Heather Bortfeld, MichaelS Beauchamp, and John S Oghalai. Auditory cortex activation to natural speech and simulated cochlear implant speech measured with functional near-infrared spectroscopy. Hearing Research, 309:84–93, 2014. doi:10.1016/j.heares.2013.11.007.）
sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)
（此图为SCI可视化）
其中sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)为计算各通道的头皮耦合指数，可以点开SCI变量查看各通道的头皮耦合指数。
随后使用以下语句标记SCI指数小于0.5的通道为坏通道。
raw_od.info["bads"] = list(compress(raw_od.ch_names, sci < 0.5))由于本次示例数据质量较高，没有通道被标记。（可自行尝试其他值）
也有文献依据SCI以一定比例去除坏通道，SCI阈值可根据需求或文献建议调整。
mne.preprocessing.nirs.scalp_coupling_index函数还可有其他参数设置，具体见网址：https://mne.tools/dev/generated/mne.preprocessing.nirs.scalp_coupling_index.html#examples-using-mne-preprocessing-nirs-scalp-coupling-index
sci = mne.preprocessing.nirs.scalp_coupling_index(raw)完整语句（默认参数设置）为：
mne.preprocessing.nirs.scalp_coupling_index(raw, l_freq=0.7, h_freq=1.5, l_trans_bandwidth=0.3, h_trans_bandwidth=0.3, verbose=False)可根据分析需求改动l_freq、h_freq等参数计算SCI指数。
（5）数据转换并移除心率噪声（梅耶尔波）
使用修正后的比尔-朗伯定律将光密度数据转换为血红蛋白浓度。
raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od, ppf=0.1)
（此图为血红蛋白浓度）
注：如前面通过SCI指数标记坏通道，在此步坏通道会显示灰色，如下所示：
研究者认为近红外研究关注的血流动力学响应的频率低于0.5Hz，而1Hz左右的梅耶尔波噪声可通过低通滤波器删除，同时通过高通滤波器移除缓慢的信号漂移。
raw_haemo_unfiltered = raw_haemo.copy()raw_haemo.filter后的滤波参数可根据分析需求改动。
出图可看到滤波效果：
可看到噪声频段显著降低。
同时右下角控制台输出FIR滤波器参数等。
（6）提取分段
提取事件相关分段：
在进行后续操作前，先对各Mark进行可视化，可通过后两句代码检查Mark数量、时间点是否有误。
检查无误后定义分段时长、拒绝标准、基线校正并提取分段。
reject_criteria = dict(hbo=80e-6)
右下角控制台将输出被拒绝分段情况。
使用以下语句，可视化被拒绝的分段。
epochs.plot_drop_log()
结合mne.Epochs具体用法可知reject_criteria = dict(hbo=80e-6)指定了当某分段内HBO信号值最大差值超过80e-6则剔除该分段。
tmin, tmax = -5, 15指定分段范围（此句为Mark前5s到Mark后15s）。
其中mne.Epochs具体用法可见链接：
https://mne.tools/stable/generated/mne.Epochs.html#mne.Epochs
可改动内部参数尝试不同标准对数据的影响。
（7）可视化分段后血氧浓度变化
检查跨分段中血氧响应的一致性。
使用以下语句，可视化实验条件（两个tapping Mark）的血氧浓度变化信息。
epochs["Tapping"].plot_image(
可通过上图观察血氧变换的峰值等。
上彩图图横坐标为时间，纵坐标为分段，可视化了血氧浓度值。
下波形图横坐标为时间，纵坐标为血氧单位，可视化了各分段的血氧浓度变化（中间黑线为均值）。
以下语句用于可视化控制（control Mark）条件数据。
epochs["Control"].plot_image(
通过以下语句，检查跨通道中血氧响应的一致性。
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(15, 6))可视化控制条件及实验条件下血氧响应变化，检查跨通道一致性。
（8）绘制HRF图、地形图等
使用以下语句绘制各条件下的HRF均值图，可以直观看到血氧浓度随时间的变换。
evoked_dict = {
使用以下语句绘制各通道HBO波形图及各时间节点地形图。
times = np.arange(-3.5, 13.2, 3.0)
其中times = np.arange(-3.5, 13.2, 3.0)指定画地形图从-3.5时间节点开始每隔3秒，到13.2为止，可自行改变参数绘制不同时间点地形图。
（9）对比左右手trapping条件
通过可视化左右手trapping条件地形图，对比不同时间点差异。使用以下语句分别可视化通道HBO和HRO均值地形图，出图顺序与代码语句顺序一致。
times = np.arange(4.0, 11.0, 1.0)
绘制差值图。
其中ts = 9.0语句指定绘制时间点。
（10）绘制波形图
可绘制所有通道HRF图。
点击波形可查看大图。
本文末尾特别感谢茗创科技机器学习授课李老师、核磁部分小齐老师、助教洋芋老师、编辑RR老师的帮助~助力本教程的诞生！
小伙伴们关注茗创科技，将第一时间收到精彩内容推送哦～
