前言

前文近紅外數(shù)據(jù)的預(yù)處理和平均（上）提到fNIRS是一種評(píng)估氧和脫氧血紅蛋白濃度變化的方法，可與fMRI相媲美。fNIRS的不足是它的空間分辨率比fMRI差，但其優(yōu)點(diǎn)是有更高的時(shí)間分辨率，并允許測(cè)量無法通過fMRI掃描儀測(cè)試的人群的大腦活動(dòng)。使用多通道近紅外記錄可以讓研究人員考察所發(fā)現(xiàn)的與事件相關(guān)的大腦活動(dòng)是局部的(主要在一個(gè)通道中發(fā)現(xiàn))，還是多個(gè)通道中都能觀察到。雖然受限于測(cè)量的空間分辨率，但是還是可以確定效應(yīng)的位置。

通常，探測(cè)器可以檢測(cè)來自多個(gè)源光電二極管的光。NIRS系統(tǒng)的不同之處在于如何實(shí)現(xiàn)這一點(diǎn)：例如，有些系統(tǒng)對(duì)每個(gè)源光電二極管使用略有不同的源波長(zhǎng)，而另一些系統(tǒng)則改變每個(gè)源光電二極管的脈沖頻率。系統(tǒng)的靈活性也各不相同：有些只允許特定的探測(cè)器“觀察”特定的光源，另一些則允許更多不同的光源和光電極的組合。有些靈活的系統(tǒng)允許這樣的放置，例如源光電極位于x，y位置，而在x+1，y位置放置探測(cè)器以及在x+2，y位置放置另一個(gè)探測(cè)器。這樣的排列方式能夠使研究人員發(fā)現(xiàn)近端、遠(yuǎn)端效應(yīng)，即在鄰近的源端和探測(cè)器之間觀察到的效應(yīng)，以及在較深處的大腦活動(dòng)，即在距離較遠(yuǎn)的源端和探測(cè)器之間觀察到的效應(yīng)。將短通道分離與正?；蜉^長(zhǎng)的通道分離相結(jié)合(通道分離是指源和探測(cè)器之間的距離)是一種消除皮膚和運(yùn)動(dòng)等其他因素偽影的有效方法。在進(jìn)行這些更復(fù)雜的分析類型之前，首先需要熟悉分析的多個(gè)通道，這是本文的目標(biāo)。

多通道近紅外數(shù)據(jù)的預(yù)處理和平均

本文將處理由多個(gè)通道組成的功能性近紅外光譜(fNIRS)數(shù)據(jù)集。包括讀取原始數(shù)據(jù)，查看設(shè)置和數(shù)據(jù)，結(jié)合多通道設(shè)置的特定程序?qū)?shù)據(jù)進(jìn)行預(yù)處理，并探索使用不同的方法來可視化響應(yīng)時(shí)間和地形圖。注意，其他NIRS系統(tǒng)的文件格式與這里使用的不同，本文可能不會(huì)詳細(xì)的講解不同格式數(shù)據(jù)的轉(zhuǎn)換，因?yàn)檫@不是本教程的目標(biāo)，但是這種總體分析結(jié)構(gòu)同樣適用于其他近紅外系統(tǒng)。

本教程所使用的數(shù)據(jù)集

本文使用的數(shù)據(jù)可以從https://download.fieldtriptoolbox.org/tutorial/nirs_multichannel/上獲得。對(duì)于XML文件，需要右鍵單擊并使用另存為選項(xiàng)。

下載完成之后，你可以在文件夾中看到以下文件:

LR-01-2015-06-01-0002.oxy3LR-02-2015-06-08-0001.oxy3LR-03-2015-06-15-0001.oxy3LR-04-2015-06-17-0001.oxy3LR-05-2015-06-23-0001.oxy3nirs_48ch_layout.matoptodetemplates.xml
Oddball任務(wù)

參與者在主動(dòng)(目標(biāo)探測(cè))傾聽情境下完成事件相關(guān)的聽覺Oddball范式任務(wù)，同時(shí)記錄左、右顳葉皮層內(nèi)氧合(HbO2)和脫氧血紅蛋白(HbR)的變化。只要有聲音出現(xiàn)，就向數(shù)據(jù)采集系統(tǒng)發(fā)送一個(gè)TTL脈沖。這種脈沖能夠精確地對(duì)刺激進(jìn)行計(jì)時(shí)，并與fNIRS的反應(yīng)同步。刺激由標(biāo)準(zhǔn)的1000Hz的音調(diào)或偏差的1500Hz的音調(diào)組成。刺激間隔為150ms，偏差為連續(xù)一列刺激中的第3~12個(gè)刺激。

fNIRS測(cè)量fNIRS數(shù)據(jù)的記錄使用了Artinis公司的4個(gè)Oxymon系統(tǒng)，這些系統(tǒng)連接起來可以實(shí)現(xiàn)48通道的記錄。光電纖維的一半(n=16)放置在左顳葉皮層，另一半放置在右顳葉皮層。源極和探測(cè)器光電極的距離分別為3cm和1.5cm的深通道和淺通道。采樣率降至250Hz。觸發(fā)器事件記錄在ADC通道1(標(biāo)準(zhǔn))和通道2(偏差)中。

分析步驟

根據(jù)數(shù)據(jù)和實(shí)驗(yàn)設(shè)計(jì)的不同，分析可以有許多不同的方式和順序。本教程的步驟順序如下(圖2):1、讀取連續(xù)數(shù)據(jù)；2、優(yōu)先降采樣，以減少內(nèi)存需求；3、刪除不良通道；4、定義epochs；5、將光密度轉(zhuǎn)化為含氧血紅蛋白(oxyHb)和脫氧血紅蛋白(deoxyHb)濃度的變化；6、將功能響應(yīng)與系統(tǒng)反應(yīng)分開，可以使用以下方法其中之一，或者組合使用；

①濾波

②減去參考通道

③每個(gè)通道都有反相關(guān)的oxyHb/deoxyHb

7、將每種條件下的試次數(shù)據(jù)進(jìn)行平均；8、結(jié)果可視化。

讀取數(shù)據(jù)并降采樣

首先需要將數(shù)據(jù)讀入到MATLAB工作區(qū)中，然后執(zhí)行ft_preprocessing：

cfg = [];cfg.dataset = 'LR-01-2015-06-01-0002.oxy3';data_raw = ft_preprocessing(cfg);

你會(huì)在命令窗口中看到如下內(nèi)容：

processing channel { 'Rx1a-Tx1 [844nm]' ... 'ADC007' 'ADC008' }reading and preprocessingreading and preprocessing trial 1 from 1the call to "ft_preprocessing" took 10 seconds and required the additional allocation of an estimated 732 MB

結(jié)構(gòu)data_raw包含關(guān)于實(shí)驗(yàn)的所有數(shù)據(jù)和信息，都存儲(chǔ)在單獨(dú)的字段中。要從如上所示的數(shù)據(jù)文件中查看通道布局，可以使用：

cfg = [];cfg.opto = 'LR-01-2015-06-01-0002.oxy3';ft_layoutplot(cfg);

檢測(cè)觸發(fā)器

對(duì)于標(biāo)準(zhǔn)音調(diào)刺激的時(shí)間在通道ADC001中表示，對(duì)于偏差刺激，則在通道ADC002中表示。這些通道可以通過查看label字段找到:

find(strcmp(data_raw.label,'ADC001'))find(strcmp(data_raw.label,'ADC002'))

將ADC001和ADC002的數(shù)據(jù)繪制成下圖，顯示模擬電壓的TTL脈沖。接下來，當(dāng)對(duì)數(shù)據(jù)進(jìn)行分段時(shí)，可以使用一個(gè)自動(dòng)化程序來查找這些標(biāo)記的更改。

figure; hold on% plot the voltage of ADC001 and ADC002% increase the scale of ADC002 a little bit to make it more clear in the figureplot(data_raw.time{1}, data_raw.trial{1}(97,:)*1.0, 'b-')plot(data_raw.time{1}, data_raw.trial{1}(98,:)*1.1, 'r:')

FieldTrip自動(dòng)檢測(cè)ADC通道中的起始點(diǎn)，并將ADC通道中向上的一側(cè)表示為事件。

event = ft_read_event('LR-01-2015-06-01-0002.oxy3')

接下來，將使用這些事件來定義感興趣的部分，并從連續(xù)數(shù)據(jù)中剪切出標(biāo)準(zhǔn)試次和偏差試次。如前所述，當(dāng)前存儲(chǔ)的fNIRS數(shù)據(jù)為250Hz。你可以用以下代碼進(jìn)行查看：

data_raw.fsample

為了節(jié)省內(nèi)存并使后續(xù)的處理更快，可以使用ft_resampledata將數(shù)據(jù)降采樣至10Hz。

cfg = [];cfg.resamplefs = 10;data_down = ft_resampledata(cfg, data_raw);
在EEG或MEG中，每個(gè)試次的片段通常在一秒左右，試次不重疊。但在fNIRS中，需要非常長(zhǎng)的片段來進(jìn)行分析，因?yàn)镠RF比較緩慢。由于該實(shí)驗(yàn)中的聽覺刺激相互間的快速跟隨，導(dǎo)致實(shí)驗(yàn)的數(shù)據(jù)段重疊。這些重疊的段會(huì)影響內(nèi)存運(yùn)行效率，因此需要使用降采樣。另一種選擇是跳過對(duì)標(biāo)準(zhǔn)音調(diào)的反應(yīng)的處理，只處理偏差音調(diào)(因?yàn)槟壳胺治鲋粚?duì)偏差數(shù)據(jù)感興趣)。

接下來，繪制數(shù)據(jù)，讓我們來看看數(shù)據(jù)是什么樣子。在cfg.preproc中，可以為動(dòng)態(tài)預(yù)處理指定一些選項(xiàng)。這里，將對(duì)數(shù)據(jù)進(jìn)行去均值處理，即減去平均值。在cfg.preproc中指定的選項(xiàng)與ft_preprocessing中的選項(xiàng)基本相同，不同的是，在當(dāng)前指令ft_databrowser中，去均值化僅用于繪圖，數(shù)據(jù)本身保持不變。

cfg = [];cfg.preproc.demean = 'yes';cfg.viewmode = 'vertical';cfg.continuous = 'no';cfg.ylim = [ -0.003 ? 0.003 ];cfg.channel = 'Rx*'; % only show channels starting with Rxft_databrowser(cfg, data_down);

這...太噪了！但不要放棄，還有拯救的機(jī)會(huì)。在接下來的步驟中，可以去除大部分噪聲。由于我們對(duì)血流動(dòng)力學(xué)響應(yīng)中非常緩慢的變化也不感興趣，所以可以通過高通濾波安全地去除低頻信息。

cfg = [];cfg.hpfilter = 'yes';cfg.hpfreq = 0.01;data_flt = ft_preprocessing(cfg,data_down);

這一步驟消除了通道間血流動(dòng)力學(xué)響應(yīng)的一些可變性。接下來，繪制過濾后的數(shù)據(jù)，看看改善情況。

cfg = [];cfg.preproc.demean = 'yes';cfg.viewmode = 'vertical';cfg.continuous = 'no';cfg.ylim = [ -0.003 ? 0.003 ];cfg.channel = 'Rx*'; % only show channels starting with Rxft_databrowser(cfg, data_flt);

分段

對(duì)數(shù)據(jù)進(jìn)行分段，以獲得感興趣的時(shí)間段，然后再進(jìn)一步處理數(shù)據(jù)。在該實(shí)驗(yàn)中，興趣段是刺激開始前5s到刺激后20s之間的一段時(shí)間?？梢允褂胒t_redefinetrial函數(shù)刪除數(shù)據(jù)中的片段，和使用ft_definetrial來確定段。

event = ft_read_event('LR-01-2015-06-01-0002.oxy3');
adc001 = find(strcmp({event.type}, 'ADC001'));adc002 = find(strcmp({event.type}, 'ADC002'));
% get the sample number in the original data% note that we transpose them to get columnssmp001 = [event(adc001).sample]';smp002 = [event(adc002).sample]';
factor = data_raw.fsample / data_down.fsample
% get the sample number after downsamplingsmp001 = round((smp001-1)/factor + 1);smp002 = round((smp002-1)/factor + 1);
pre ? ?= ?round( 5*data_down.fsample);post = ?round(20*data_down.fsample);offset = -pre; % see ft_definetrial
trl001 = [smp001-pre smp001+post];trl002 = [smp002-pre smp002+post];
% add the offsettrl001(:,3) = offset;trl002(:,3) = offset;
trl001(:,4) = 1; % add a column with the condition numbertrl002(:,4) = 2; % add a column with the condition number
% concatenate the two conditions and sort themtrl = sortrows([trl001; trl002])
% remove trials that stretch beyond the end of the recordingsel = trl(:,2)<size(data_down.trial{1},2);trl = trl(sel,:);
cfg ? ? = [];cfg.trl = trl;data_epoch = ft_redefinetrial(cfg,data_down);

如果輸入data_epoch，應(yīng)該會(huì)在命令窗口中看到這個(gè):

data_epoch =
?struct?with?field ??????????hdr:?[1x1?struct]????????trial:?{1x597?cell}?????????time:?{1x597?cell}??????fsample:?10????????label:?{104x1?cell}?????????opto:?[1x1?struct]????trialinfo:?[597x1?double]???sampleinfo:?[597x2?double]??????????cfg:?[1x1?struct]

值得注意的是，trial和time字段現(xiàn)在都有1x597個(gè)單元陣列。這與呈現(xiàn)的597種刺激相對(duì)應(yīng)。在data_epoch.trialinfo中：存儲(chǔ)有關(guān)刺激類型的信息(事件1或事件2)。因此，可以找到哪個(gè)單元是第一個(gè)偏差刺激：

idx = find(data_epoch.trialinfo==2, 1, 'first')

運(yùn)行以上行代碼，將得到：

idx = ? ? 8
讓我們通過繪制平均光密度圖來看看第一個(gè)偏差附近的情況：

cfg = [];cfg.channel = 'Rx*';cfg.trials = 8;cfg.baseline = 'yes';ft_singleplotER(cfg, data_epoch)

刪除不良通道

首先移除與頭皮接觸不良、因此產(chǎn)生不良信號(hào)的通道。從光密度軌跡可以估計(jì)光電二極管和頭皮之間的耦合是否良好，因?yàn)閬碜悦總€(gè)光電二極管的兩個(gè)信號(hào)(對(duì)應(yīng)于兩個(gè)波長(zhǎng))應(yīng)該有一個(gè)正相關(guān)的心跳。如果相關(guān)性很小或?yàn)樨?fù)，那么將在進(jìn)一步處理中排除該光電極。這一步可以在ft_nirs_scalpcouplingindex中實(shí)現(xiàn)。

cfg = [];data_sci = ft_nirs_scalpcouplingindex(cfg, data_epoch);

可以看到在data_sci.label中刪除了一些通道，現(xiàn)在只有86個(gè)標(biāo)簽，而不是104個(gè)：

data_sci =
?struct?with?field ?????????hdr:?[1x1?struct]???????trial:?{1x597?cell}????????time:?{1x597?cell}?????fsample:?10???????label:?{86x1?cell}????????opto:?[1x1?struct]???trialinfo:?[597x1?double]??sampleinfo:?[597x2?double]?????????cfg:?[1x1?struct]
刪除偽影

這里已經(jīng)通過分段去除了主要的運(yùn)動(dòng)偽影，并刪除了耦合不良的光電極。因此，對(duì)于該數(shù)據(jù)集，可以忽略此步驟。

將光密度轉(zhuǎn)化為氧和脫氧血紅蛋白濃度的變化

通過ft_nirs_transform_ODs將光密度轉(zhuǎn)換為含氧和脫氧血紅蛋白的濃度。

cfg = [];cfg.target = {'O2Hb', 'HHb'};cfg.channel = 'nirs'; % e.g., one channel incl. wildcards, you can also use ?all? to select all NIRS channelsdata_conc = ft_nirs_transform_ODs(cfg, data_sci);

使用ft_singleplotER再次檢查數(shù)據(jù)。你可以在信號(hào)中看到一個(gè)明顯的心跳。

將功能響應(yīng)與系統(tǒng)反應(yīng)分開低通濾波

心跳并不是該數(shù)據(jù)集對(duì)應(yīng)的研究中感興趣的信號(hào)，所以將使用低通濾波對(duì)數(shù)據(jù)進(jìn)行過濾，使其低于心跳頻率(約1Hz)。

cfg = [];cfg.lpfilter = 'yes';cfg.lpfreq = 0.8;data_lpf = ft_preprocessing(cfg, data_conc);

這里的平均濃度變化是一個(gè)完美的血液動(dòng)力學(xué)響應(yīng)例子。信號(hào)在刺激開始時(shí)上升，在4s左右達(dá)到峰值，然后再次下降。

繪制結(jié)果

首先，運(yùn)行ft_timelockanalysis計(jì)算均值。ft_timelockanalysis默認(rèn)對(duì)所有試次進(jìn)行平均。若想對(duì)偏差試次和標(biāo)準(zhǔn)試次分別進(jìn)行平均，需要告訴代碼哪些試次屬于標(biāo)準(zhǔn)刺激哪些屬于偏差刺激。有關(guān)條件的信息存儲(chǔ)在trialinfo中，其中1表示標(biāo)準(zhǔn)，2表示偏差。

cfg = [];cfg.trials = find(data_lpf.trialinfo(:,1) == 1);timelockSTD = ft_timelockanalysis(cfg, data_lpf);

然后，使用ft_timelockbaseline進(jìn)行基線校正。刺激前5s將被用作基準(zhǔn)的時(shí)間窗。

cfg = [];cfg.baseline = [-5 0];timelockSTD = ft_timelockbaseline(cfg, timelockSTD);

對(duì)偏差刺激執(zhí)行同樣的步驟。

cfg = [];cfg.trials = find(data_lpf.trialinfo(:,1) == 2);timelockDEV = ft_timelockanalysis(cfg, data_lpf);cfg ? ? ? ? ? = [];cfg.baseline = [-5 0];timelockDEV = ft_timelockbaseline(cfg, timelockDEV);

該通道布局可以使用MATLAB函數(shù)load讀取nirs_48ch_layout.mat文件。該文件包含一個(gè)名為lay的結(jié)構(gòu)。

load('nirs_48ch_layout.mat')figure; ft_plot_layout(lay) % note that O2Hb and HHb channels fall on top of each other

有許多FieldTrip選項(xiàng)可用于結(jié)果可視化，例如ft_singleplotER(ER表示事件相關(guān))，該函數(shù)允許繪制單個(gè)通道，以及ft_multiplotER，該函數(shù)允許繪制多個(gè)通道。ft_multiplotER還可以在交互模式下使用，以選擇感興趣的數(shù)據(jù)片段(例如選擇特定的通道和特定的時(shí)間窗)。

但值得注意的是，要運(yùn)行ft_multiplotER，需要使用cfg.layout = lay將FieldTrip指向布局結(jié)構(gòu)。

cfg = [];cfg.showlabels = 'yes';cfg.layout = lay; ? ? ?% you could also specify the name of the mat filecfg.interactive = 'yes';cfg.linecolor = 'rb';cfg.colorgroups(contains(timelockDEV.label, 'O2Hb')) = 1; % these will be redcfg.colorgroups(contains(timelockDEV.label, 'HHb')) = 2; % these will be blueft_multiplotER(cfg, timelockDEV);

此外，還可以在某個(gè)時(shí)間點(diǎn)或時(shí)間窗上生成信號(hào)的空間地形圖。時(shí)間窗可以通過cfg.xlim = [5 7];命令進(jìn)行設(shè)置，響應(yīng)強(qiáng)度的范圍可以設(shè)置為-0.2到0.2，但這取決于你研究中的數(shù)據(jù)情況：例如，許多fNIRS研究人員使用block設(shè)計(jì)，根據(jù)block持續(xù)時(shí)間，響應(yīng)可能具有更大的振幅。

cfg = [];cfg.layout = lay; ? ? ?% you could also specify the name of the mat filecfg.marker = 'labels';cfg.xlim = [5 7];cfg.zlim = [-0.2 0.2];cfg.channel ?= '* [O2Hb]';figure; subplot(1,2,1);ft_topoplotER(cfg, timelockDEV);title('[O2Hb]');cfg.channel ?= '* [HHb]';subplot(1,2,2);ft_topoplotER(cfg, timelockDEV);title('[HHb]');

參考來源：Preprocessing and averaging of multi-channel NIRS data.https://www.fieldtriptoolbox.org/tutorial/nirs_multichannel/https://download.fieldtriptoolbox.org/tutorial/nirs_multichannel/

小伙伴們點(diǎn)星標(biāo)關(guān)注茗創(chuàng)科技，將第一時(shí)間收到精彩內(nèi)容推送哦～