Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
tutorial:nirs_multichannel [2018/07/11 12:56]
robert
tutorial:nirs_multichannel [2018/09/19 10:16] (current)
robert [Separate functional from systemic responses]
Line 1: Line 1:
 {{tag>​tutorial nirs preprocessing nirs-multichannel}} {{tag>​tutorial nirs preprocessing nirs-multichannel}}
 +
  
 ====== Preprocessing and averaging of multi-channel NIRS data ====== ​ ====== Preprocessing and averaging of multi-channel NIRS data ====== ​
 +
 +<note important>​
 +This tutorial is still under development.
 +</​note>​
  
 =====Introduction===== =====Introduction=====
-At the end of this tutorial, you will be able to process a functional near-infrared spectroscopy (fNIRS) data set consisting of multiple channels ​and you will be able to visualize ​the temporal and spatial aspects of the response.+In this tutorial, you will process a functional near-infrared spectroscopy (fNIRS) data set consisting of multiple channels. We will read in the raw data, have a look at the setup and the data, preprocess the data incorporating specific procedures for multichannel setups, and explore different methods of visualizing ​the temporal and spatial aspects of the response.
  
-We suggest to first read the [[http://​www.fieldtriptoolbox.org/​getting_started/​artinis|Getting started with Artinis]] page to get the details on the recording system. Furthermore,​ we suggest you follow the [[/​tutorial/​nirs_singlechannel|NIRS single channel]] tutorial to learn how to pre-process a simpler fNIRS data set with only a a single channel. ​+We suggest to first read the [[/​getting_started/​artinis|Getting started with Artinis]] page to get the details on the recording system. Furthermore,​ we suggest you follow the [[/​tutorial/​nirs_singlechannel|NIRS single channel]] tutorial to learn how to pre-process a simpler fNIRS data set with only a a single channel. ​
  
 Note that other NIRS systems will have a different file format than the ones used here and that this tutorial might not directly translate. The global structure however applies to other NIRS systems as well. Note that other NIRS systems will have a different file format than the ones used here and that this tutorial might not directly translate. The global structure however applies to other NIRS systems as well.
  
 ====Background==== ​ ====Background==== ​
-In this tutorial we will read in the unprocessed raw datahave a look at the setup and the datapreprocess ​the data incorporating specific procedures for multichannel setups, and explore different methods of visualizing ​and analyzing the data+ 
 +fNIRS is a method to assess changes in oxygenated and deoxygenated hemoglobin concentrations and as such is comparable to fMRI. Disadvantage of fNIRS over fMRI is that it has a (much!) poorer spatial resolution than fMRI, but, the advantage is that it has a better temporal resolution and allows the brain measurements of populations that typically cannot easily be scanned using fMRI. More about the method can be read on the page of the [[tutorial:​nirs_singlechannel|single channel tutorial]]. 
 + 
 +Using multiple channels allows the researcher to investigate whether the found event-related brain activity is local, found in mainly one channel, or more widespread and can be observed in multiple channels. Moreover, the location of the effect can then be determined, though of course one is limited by the spatial resolution of the measurement.  
 + 
 +Typically, a detector optode can detect light that originates from multiple source optodes. NIRS systems differ in how this is implemented:​ for instance, some use slightly different source wavelengths per source optode, others vary the pulse frequency per source optode. Systems also vary in their flexibility:​ some allow only specific detectors to "look" ​at specific sources, others allow many more different combinations of sources ​and optodes. A flexible system can allow you to put for instance a source optode on location x,y, a detector optode on location x+1,y and another detector optode on location x+2,y. Such an arrangement allows ​the researcher to find proximalshallow, effects, namely effects observed between ​the source and detector pair that are positioned next to each otheras well as brain activity somewhat deeper in the brain, namely effects observed between the source ​and detector pair that are further apart. Combining short channel separations with normal or longer channel separations (channel separation being the distance between source ​and detector) is a helpful way to cancel out artifacts from among others the skin and motion. Before you can turn to these more sophisticated types of analyses, you need to familiarize yourself with analyzing ​multiple channels, which is the aim of the current tutorial.
  
 ====Dataset information==== ====Dataset information====
Line 34: Line 44:
 ===Oddball task=== ===Oddball task===
  
-The participant was engaged in a basic event-related auditory oddball paradigm in an active (target detection) listening situation, while recording changes in oxygenated (HbO2) and deoxygenated hemoglobin (HbR) within left and right temporal cortices. Whenever a tone was presented, a transistor–transistor logic (TTL) pulse was sent to the data acquisition system (the ADC channels of the Oxymon system, see below). This pulse enables us to exactly time and synchronize the stimuli with the fNIRS responses. The stimulus consisted of either a standard 1000 Hz tone or a deviant 1500 Hz tone. Stimuli were presented at intervals of 150 ms, and the deviant was presented as the 3rd to 12th stimuluss ​in a row.+The participant was engaged in a basic event-related auditory oddball paradigm in an active (target detection) listening situation, while recording changes in oxygenated (HbO2) and deoxygenated hemoglobin (HbR) within left and right temporal cortices. Whenever a tone was presented, a transistor–transistor logic (TTL) pulse was sent to the data acquisition system (the ADC channels of the Oxymon system, see below). This pulse enables us to exactly time and synchronize the stimuli with the fNIRS responses. The stimulus consisted of either a standard 1000 Hz tone or a deviant 1500 Hz tone. Stimuli were presented at intervals of 150 ms, and the deviant was presented as the 3rd to 12th stimulus ​in a row.
  
 ===fNIRS Measurement=== ===fNIRS Measurement===
 The fNIRS data was recorded using 4 connected Oxymon systems from Artinis to enable a 48 channel recording. ​ The fNIRS data was recorded using 4 connected Oxymon systems from Artinis to enable a 48 channel recording. ​
  
-Figure 1. Oxymons and cap+FIXME add a photo of the cap as Figure 1
  
 Half of the optode fibers (n = 16) were placed over left temporal cortex, the other half over the right temporal cortex (Fig. 1 grey text). Of the optodes, half were detectors (or receivers, Rx), and the other half were sources (or transmitters,​ Tx). The source and detector optodes were positioned such that there were deep and shallow channels of 3 cm and 1.5 cm, respectively (black text indicating both receivers and transmitters).  ​ Half of the optode fibers (n = 16) were placed over left temporal cortex, the other half over the right temporal cortex (Fig. 1 grey text). Of the optodes, half were detectors (or receivers, Rx), and the other half were sources (or transmitters,​ Tx). The source and detector optodes were positioned such that there were deep and shallow channels of 3 cm and 1.5 cm, respectively (black text indicating both receivers and transmitters).  ​
Line 62: Line 72:
 {{tutorial:​nirs_tut2_multichannel_analysis_steps.png?&​400}} {{tutorial:​nirs_tut2_multichannel_analysis_steps.png?&​400}}
  
-**//​Overview of the fNIRS analysis procedure for the data set of this tutorial.//​**+**//Figure2; ​Overview of the fNIRS analysis procedure for the data set of this tutorial.//​**
  
 ====Read data and downsample==== ====Read data and downsample====
-We first need to read in the data into the MATLAB workspace, by executing ft_preprocessing:​+We first need to read in the data into the MATLAB workspace, by executing ​**[[/​reference/​ft_preprocessing]]**:
  
 <​code>​ <​code>​
-fname                = '​LR-01-2015-06-01-0002.oxy3';​ 
 cfg          = []; cfg          = [];
-cfg.dataset ​        ​ = ​fname;+cfg.dataset ​        ​ = ​'​LR-01-2015-06-01-0002.oxy3'​;
 data_raw ​           = ft_preprocessing(cfg);​ data_raw ​           = ft_preprocessing(cfg);​
 </​code>​ </​code>​
Line 87: Line 96:
 The structure **data_raw** contains all data and information about the experiment, all stored in separate fields. The structure **data_raw** contains all data and information about the experiment, all stored in separate fields.
  
-Note: for information about FieldTrip data structures and their fields, see this [[faq/​how_are_the_various_data_structures_defined|frequently asked question]].+<note note> 
 +For information about FieldTrip data structures and their fields, see this [[faq/​how_are_the_various_data_structures_defined|frequently asked question]]. 
 +</​note>​
  
-The layout ​representation ​as shown above can be obtainedby executing+ 
-  cfg.optofile ​           fname;+To retrieve the layout ​from the data file as shown above, ​you can use 
 + 
 +  cfg           = []; 
 +  cfg.optofile ​ '​LR-01-2015-06-01-0002.oxy3'​;
   ft_layoutplot(cfg);​   ft_layoutplot(cfg);​
  
 +{{tutorial:​nirs_tut2_optodepositions.png?&​400}}
 +**//Figure 3; Layout of the optode positions//​**
 +  ​
 +
 +===Trigger channels===
 Important information about the timing of the stimuli, is stored in channel ADC001 for the standard tones and in channel ADC002 for the deviant tones. The labels of the channels can be found in the field label: Important information about the timing of the stimuli, is stored in channel ADC001 for the standard tones and in channel ADC002 for the deviant tones. The labels of the channels can be found in the field label:
  
Line 100: Line 119:
 which should you tell you that the data in channel ADC001 is the 97th row in the data matrix, which is stored in data_raw.trial{1}. which should you tell you that the data in channel ADC001 is the 97th row in the data matrix, which is stored in data_raw.trial{1}.
  
-Plotting the data from ADC001 and ADC002 will yield the figure below, showing the TTL pulses as analog voltages. ​At this stage, we do not have exact timings (in seconds or in samples) ​of the stimulus onsetsWe will obtain those later in the tutorial ​when epoching the data. +Plotting the data from ADC001 and ADC002 will yield the figure below, showing the TTL pulses as analog voltages. ​Stimulus onset is marked by an abrupt increase ​in one of the analog channelsLater on, when epoching the data, we will use an automatic routine to find these marked changes.
 <​code>​ <​code>​
 figure; hold on figure; hold on
 % plot the voltage of ADC001 and ADC002 % plot the voltage of ADC001 and ADC002
 % ADC002 is scaled up a little bit to make it more clear % ADC002 is scaled up a little bit to make it more clear
-plot(data_raw.time{1},​ data_raw.trial{1}(97,:​)*1.0, ​b-+plot(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:)+plot(data_raw.time{1},​ data_raw.trial{1}(98,:​)*1.1, ​'r:')
 </​code>​ </​code>​
 +
  
 {{tutorial:​nirs_tut2_datatrigger.png?&​400}} {{tutorial:​nirs_tut2_datatrigger.png?&​400}}
  
-**//Oddball paradigm. All stimuli onsets are indicated by the blue lines. Red dotted lines indicate onsets of the deviants. You can see that there are four blocks of events.//**+**//Figure 4; Oddball paradigm ​trigger. All stimuli onsets are indicated by the blue lines. Red dotted lines indicate onsets of the deviants. You can see that there are four blocks of events.//**
  
-We can look into a smaller time window (in the code snippet below, we look at 355 and 365 seconds in the experiment) to better see what is going on: 
  
-<code+<note exercise
-idx1 = find(data_raw.time{1} > 355,1); +** Exercise ​** 
-idx2 = find(data_raw.time{1} > 365,1); +Zoom in on 355 to 365 seconds to better see what is going on.  All stimuli onsets are indicated by the blue lines Red dotted lines indicate onsets of the deviants ​(the oddballs). Can you now better spot the oddball? 
-figure; hold on +</note>
-plot(data_raw.time{1}(idx1:​idx2),​ data_raw.trial{1}(98,​idx1:​idx2)*1.1, '​r:'​) +
-plot(data_raw.time{1}(idx1:​idx2),​ data_raw.trial{1}(97,​idx1:​idx2)*1.0,​ '​b-'​) +
-</code>+
  
-{{tutorial:​nirs_tut2_oddballparadigm.png?&​400}} 
- 
-**//Oddball paradigm (zoomed). All stimuli onsets are indicated by the blue lines. ​ Red dotted lines indicate onsets of the deviants. 
-//** 
  
 FieldTrip can detect the onset in the ADC channels automatically and represent the upward going flank in the ADC channels as “event”. FieldTrip can detect the onset in the ADC channels automatically and represent the upward going flank in the ADC channels as “event”.
  
-  event = ft_read_event(fname) +  event = ft_read_event('​LR-01-2015-06-01-0002.oxy3')
-  adc001 = find(strcmp({event.type}, ‘ADC001’)) +
-  adc002 = find(strcmp({event.type},​ ‘ADC002’))+
  
-  figure; hold on +<note exercise>​ 
-  ​plot([event(adc001).sample][event(adc001).value],​ '​bx'​) +** Exercise 2 ** 
-  plot([event(adc002).sample], [event(adc002).value], 'ro')+Explore the information in the event structureHow many stimuli were played and how many oddballs? As not all events are stimuli onsetsit might help to find the oddballs by running adc002 = find(strcmp({event.type}, 'ADC002')); 
 +</​note>​
  
 We will use these events later to define segments of interest and to cut the standard and deviant trials out of the continuous data.  We will use these events later to define segments of interest and to cut the standard and deviant trials out of the continuous data. 
Line 151: Line 161:
   data_down ​           = ft_resampledata(cfg,​ data_raw);   data_down ​           = ft_resampledata(cfg,​ data_raw);
    
-Note: it is better to resample multiple times if the resampling factor is larger than 10, see [[https://​allsignalprocessing.com/​very-low-frequency-filtering/​|here]]+<note note> 
 +It is better to resample multiple times if the resampling factor is larger than 10, see [[https://​allsignalprocessing.com/​very-low-frequency-filtering/​|here]] 
 +</​note>​
  
-Note: the resampling also includes low-pass filtering of the data. As the new sampling rate is 10 Hz, we will lose data with frequencies larger than 5 Hz. This means we will lose a lot of information from the standards in our experiment, as they are presented near 6.7 Hz, but we keep the deviant information,​ which is presented near 0.6 Hz. For the current analysis, we are only interested in the deviant data. Just remember: be wary of filtering!+<note note> 
 +The resampling also includes low-pass filtering of the data. As the new sampling rate is 10 Hz, we will lose data with frequencies larger than 5 Hz. This means we will lose a lot of information from the standards in our experiment, as they are presented near 6.7 Hz, but we keep the deviant information,​ which is presented near 0.6 Hz. For the current analysis, we are only interested in the deviant data. Just remember: be wary of filtering! 
 +</​note>​
  
-We can now plot the data and see what it looks like. In cfg.preproc we can specify some options for on-the-fly preprocessing. Here, we will demean the data, i.e. subtract the mean value. The options you can specify in cfg.preproc are largely the same as the options for ft_preprocessing. ​  +We can now plot the data and see what it looks like. In cfg.preproc we can specify some options for on-the-fly preprocessing. Here, we will demean the data, i.e. subtract the mean value. The options you can specify in cfg.preproc are largely the same as the options for **[[reference:​ft_preprocessing|ft_preprocessing]]** with as a difference that in our current command, namely ft_databrowser,​ the demeaning is only applied for plotting, the data itself remains the same.   
  
 <​code>​ <​code>​
-cfg = [];+cfg                = [];
 cfg.preproc.demean = '​yes';​ cfg.preproc.demean = '​yes';​
-cfg.viewmode = '​vertical';​ +cfg.viewmode ​      ​= '​vertical';​ 
-cfg.continuous = '​no';​ +cfg.continuous ​    ​= '​no';​ 
-cfg.ylim = [ -0.003 ​  0.003 ]; +cfg.ylim ​          ​= [ -0.003 ​  0.003 ]; 
-cfg.channel = Rx*; % only show channels starting with Rx+cfg.channel ​       'Rx*'; % only show channels starting with Rx
 ft_databrowser(cfg,​ data_down); ft_databrowser(cfg,​ data_down);
 +</​code>​
 +
 +{{tutorial:​nirs_tut2_fig5_databrowser.png?​400}}
 +
 +**//Figure 5; Optical density traces for down-sampled data before high-pass filtering.//​**
 +
 +This is very noisy! Do not give up hope. In the next steps, you will remove most of the noise.
 +
 +As we are also not interested in very slow changes (and/or a constant offset/ DC) in the hemodynamic response, we can ‘safely’ throw away low-frequency information by high-pass filtering.
 +
 +<​code>​
 +cfg                 = [];
 +cfg.hpfilter ​       = '​yes';​
 +cfg.hpfreq ​         = 0.01; 
 +data_flt ​           = ft_preprocessing(cfg,​data_down);​
 +</​code>​
 + 
 +This step has removed some of the variability in the hemodynamic response between channels. Let’s plot the filtered data to see how things have improved.
 +
 +<​code>​
 +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 Rx
 +ft_databrowser(cfg,​ data_flt);
 </​code>​ </​code>​
  
 {{tutorial:​nirs_tut2_opticaldensitytracesafterhighpass.png?&​400}} {{tutorial:​nirs_tut2_opticaldensitytracesafterhighpass.png?&​400}}
  
-**//Optical density traces for down-sampled data after high-pass filtering. Note that the DC (offset) has been largely removed by this step.//**+**//Figure 6; Optical density traces for down-sampled data after high-pass filtering. Note that the DC (offset) has been largely removed by this step (cf. Fig. 5).//**
  
 ====Epoch==== ====Epoch====
-In the single channel tutorial, after initial preprocessing we continued with removing ‘bad’ data. In this tutorial, we will first segment the data to get the time segments of interest. The motivation here to first segment and then detect artifacts is that the largest artifacts in the data are due to motion artifacts that occur between the experimental blocks. By segmenting the data in trials, these non-relevant sections in the data are ignored and we obtain a cleaner data set already+In the single channel tutorial, after initial preprocessing we continued with removing ‘bad’ data as there were no pieces of the data that were both irrelevant (say, during a break) and very noisy. In this tutorial, we will first segment the data to get the time segments of interest ​before we move on to cleaning the data further. The motivation here to first segment and then detect artifacts is that the largest artifacts in the data are due to motion artifacts that occur between the experimental blocks. By segmenting the data in trials, these non-relevant sections in the data are ignored and we obtain a cleaner data set already.
- +
-In this experiment, the segment of interest is a period of 5 s before and 15 s after stimulus onset. The stimuli consist of pure tones of 1000 or 1500 Hz, of which the standard 1000 Hz tones are presented most often and are referred to as event 1. The deviant 1500 Hz stimuli are the most important stimuli, and are referred to as event 2.+
  
-We will cut out the segments in the data using the function **[[/​reference/​ft_redefinetrial]]**. Normally we would use **[[/​reference/​ft_definetrial]]** to determine the segments, but due to the resampling the sample indices have changed and hence we will do it by hand. +In this experiment, the segment of interest is a period of 5 s before and 20s after stimulus onset. ​We will cut out the segments in the data using the function **[[/​reference/​ft_redefinetrial]]**. Normally we would use **[[/​reference/​ft_definetrial]]** to determine the segments, but due to the resampling the sample indices have changed and hence we will do it by hand. 
  
 <​code>​ <​code>​
-event = ft_read_event(fname);+event = ft_read_event('​LR-01-2015-06-01-0002.oxy3'​);
  
 adc001 = find(strcmp({event.type},​ ‘ADC001’));​ adc001 = find(strcmp({event.type},​ ‘ADC001’));​
Line 189: Line 228:
 smp002 = [event(adc002).sample]’;​ smp002 = [event(adc002).sample]’;​
  
-factor = data_raw.fsample / data_flt.fsample+factor = data_raw.fsample / data_down.fsample
  
 % get the sample number after downsampling % get the sample number after downsampling
Line 195: Line 234:
 smp002 = round((smp002-1)/​factor + 1); smp002 = round((smp002-1)/​factor + 1);
  
-pre    =  round( 5*data_flt.fsample);​ +pre    =  round( 5*data_down.fsample);​ 
-post   ​= ​ round(20*data_flt.fsample);+post   ​= ​ round(20*data_down.fsample);
 offset = -pre; % see ft_definetrial offset = -pre; % see ft_definetrial
  
Line 213: Line 252:
  
 % remove trials that stretch beyond the end of the recording % remove trials that stretch beyond the end of the recording
-sel = trl(:,​2)<​size(data_flt.trial{1},​2);​+sel = trl(:,​2)<​size(data_down.trial{1},​2);​
 trl = trl(sel,:); trl = trl(sel,:);
  
-cfg = []+cfg     ​= []
 cfg.trl = trl cfg.trl = trl
-data_epoch = ft_redefinetrial(cfg,​data_flt);+data_epoch = ft_redefinetrial(cfg,​data_down);
 </​code>​ </​code>​
  
Line 239: Line 278:
 </​code>​ </​code>​
  
-Notably, both trial and time fields will now have 1x597 cell array (compare this to data_flt). This corresponds to the 597 stimuli that were presented. In data_epoch.trialinfo the information about the type of stimulus is stored (event 1 or event 2). Thus, we can find which of those cells belongs to the first deviant:+Notably, both trial and time fields will now have 1x597 cell array (compare this to data_down). This corresponds to the 597 stimuli that were presented. In data_epoch.trialinfo the information about the type of stimulus is stored (event 1 or event 2). Thus, we can find which of those cells belongs to the first deviant:
  
   idx = find(data_epoch.trialinfo==2,​1,'​first'​)   idx = find(data_epoch.trialinfo==2,​1,'​first'​)
Line 253: Line 292:
  
 <​code>​ <​code>​
-cfg = []; +cfg          = []; 
-cfg.channel = '​Rx*';​ +cfg.channel ​ = '​Rx*';​ 
-cfg.trials = 8;+cfg.trials ​  ​= 8;
 cfg.baseline = '​yes';​ cfg.baseline = '​yes';​
 ft_singleplotER(cfg,​ data_epoch) ft_singleplotER(cfg,​ data_epoch)
Line 262: Line 301:
 {{tutorial:​nirs_tut2_epocheddata.png?&​400}} {{tutorial:​nirs_tut2_epocheddata.png?&​400}}
  
-**//Epoched optical density data around the first deviant stimulus.//​**+**//Figure 7; Epoched optical density data around the first deviant stimulus.//​**
  
 The most obvious thing you should see, is the heartbeat. This is great! It means that your subject is alive and has some blood flowing through his/her brain (or skin). Importantly,​ this is an indicative sign of a good measurement. If you would not see this, you could throw this data in the bin (see next paragraph). The most obvious thing you should see, is the heartbeat. This is great! It means that your subject is alive and has some blood flowing through his/her brain (or skin). Importantly,​ this is an indicative sign of a good measurement. If you would not see this, you could throw this data in the bin (see next paragraph).
  
-If you look carefully, you will also see that the signal ​actually increases after sound onset at 0 sreaches a peak around 4 s, and then decays again. ​Could this be a functional response? We have to do a few more additional analysis steps, before we know for sure.+<note exercise>​ 
 +** Exercise 3 ** 
 +Inspect ​the signal ​carefully! When does it increase/​decreasewhen does it peakCould this be a functional response? We have to do a few more additional analysis steps, before we know for sure. 
 +</​note>​
  
 ====Remove bad channels==== ====Remove bad channels====
-First, we will remove the optode channels that make poor contact with the skin of the scalp yielding bad signal because of that. From the optical density traces we can estimate whether there is a good coupling between optode and scalp, because the two signals from each optode (corresponding to the two wavelengths) should have a heartbeat that is positively correlated. If the correlation is small or negative, we exclude that optode from further processing. This is implemented in ft_nirs_scalpcouplingindex. ​ [For more details see Polloniniet al. (2014), Auditory cortex activation to natural speech and simulated cochlear implant speech measured with functional near-infrared spectroscopy. Hearing Research, 309, 84-93. doi:​10.1016/​j.heares.2013.11.007]+First, we will remove the optode channels that make poor contact with the skin of the scalp yielding bad signal because of that. From the optical density traces we can estimate whether there is a good coupling between optode and scalp, because the two signals from each optode (corresponding to the two wavelengths) should have a heartbeat that is positively correlated. If the correlation is small or negative, we exclude that optode from further processing. This is implemented in **[[/​reference/​ft_nirs_scalpcouplingindex]]**.  [For more details see Polloniniet al. (2014), Auditory cortex activation to natural speech and simulated cochlear implant speech measured with functional near-infrared spectroscopy. Hearing Research, 309, 84-93. doi:​10.1016/​j.heares.2013.11.007]
  
   cfg                 = [];   cfg                 = [];
Line 292: Line 334:
  
 ====Remove artefacts==== ====Remove artefacts====
-We already removed major motion artefacts by epoching, thus removing the periods in between blocks, and by throwing away poorly coupled optodes. Therefore, this step can be ignored for this dataset ​(check ​this yourself!)+We already removed major motion artefacts by epoching, thus removing the periods in between blocks, and by throwing away poorly coupled optodes. Therefore, this step can be ignored for this dataset
 + 
 +<note exercise>​ 
 +** Exercise 4 ** 
 +We just wrote "​Therefore,​ this step can be ignored."​ Check this yourself, are there indeed no artifacts? Hint: you can use cfg.artfctdef.zvalue.interactive = '​yes';​ and [cfg, artifact] = ft_artifact_zvalue(cfg,​ data_epoch); like in Exercise 2 of the single channel tutorial. 
 +</​note>​
  
 ====Transform optical densities to oxy- and deoxyhemoglobin concentration changes==== ====Transform optical densities to oxy- and deoxyhemoglobin concentration changes====
-Like in the [[http://​www.fieldtriptoolbox.org/​tutorial/​nirs_singlechannel|single channel tutorial]], we will now convert the optical densities into oxygenated and deoxygenated hemoglobin concentrations by using **[[/​reference/​ft_transform_ODs]]**.+Like in the [[/​tutorial/​nirs_singlechannel|single channel tutorial]], we will now convert the optical densities into oxygenated and deoxygenated hemoglobin concentrations by using **[[/​reference/​ft_nirs_transform_ODs]]**.
  
 <​code>​ <​code>​
 cfg                 = []; cfg                 = [];
 cfg.target ​         = {'​O2Hb',​ '​HHb'​};​ cfg.target ​         = {'​O2Hb',​ '​HHb'​};​
-cfg.channel ​        = '​nirs';​ % e.g. one channel incl. wildcards, you can also use ?all?all NIRS channels ​are selected then+cfg.channel ​        = '​nirs';​ % e.g. one channel incl. wildcards, you can also use ?all? to select ​all nirs channels
 data_trans ​         = ft_nirs_transform_ODs(cfg,​ data_sci); data_trans ​         = ft_nirs_transform_ODs(cfg,​ data_sci);
 </​code>​ </​code>​
Line 307: Line 354:
 {{tutorial:​nirs_tut2_hemoglobinovertime.png?&​400}} {{tutorial:​nirs_tut2_hemoglobinovertime.png?&​400}}
  
-**//​Hemoglobin concentration as a function of time, averaged over all channels for the epoch around the first deviant ​(cf. previous fig.)+**//Figure 8; Hemoglobin concentration as a function of time, averaged over all channels for the epoch around the first deviant
 //** //**
  
Line 317: Line 364:
 cfg.lpfilter ​       = '​yes';​ cfg.lpfilter ​       = '​yes';​
 cfg.lpfreq ​         = 0.8; cfg.lpfreq ​         = 0.8;
-cfg.bpfilttype ​ = '​fir';​+cfg.bpfilttype ​         = '​fir';​
 data_lpf ​           = ft_preprocessing(cfg,​data_trans);​ data_lpf ​           = ft_preprocessing(cfg,​data_trans);​
 </​code>​ </​code>​
Line 324: Line 371:
 {{tutorial:​nirs_tut2_hemoglobinovertimeafterlowpass.png?&​400}} {{tutorial:​nirs_tut2_hemoglobinovertimeafterlowpass.png?&​400}}
  
-**//​Low-pass filtered hemoglobin concentrations (cf. two previous figs.).//**+**//Figure 9; Low-pass filtered hemoglobin concentrations (cf. two previous figs.).//**
  
-===Layout=== 
-To retrieve the layout from the data file, you can use:  
- 
-  cfg                  = []; 
-  cfg.optofile ​   = '​LR-01-2015-06-01-0002.oxy3';​ 
-  ft_layoutplot(cfg);​ 
- 
-{{tutorial:​nirs_tut2_optodepositions.png?&​400}} 
- 
-**//Layout of the optode positions//​** 
  
 ====Plot results==== ====Plot results====
 Now that we obtained the functional responses, the next step is to average over trials and to visualize the results. Now that we obtained the functional responses, the next step is to average over trials and to visualize the results.
-First, we will run ft_timelockanalysis to compute the average. The default behavior of the **[[/​reference/​ft_timelockanalysis]]** is to average across all trials. We want to make a separate average for the deviant and for the standard trials, hence we need inform the code which trials belong to the standard stimuli and which belong to the deviants. Information about the conditions is stored in the tiralinfo, where 1 represents the standards, and 2 represents the deviants.+First, we will run **[[/​reference/​ft_timelockanalysis]]** to compute the average. The default behavior of the **[[/​reference/​ft_timelockanalysis]]** is to average across all trials. We want to make a separate average for the deviant and for the standard trials, hence we need inform the code which trials belong to the standard stimuli and which belong to the deviants. Information about the conditions is stored in the trialinfo, where 1 represents the standards, and 2 represents the deviants.
  
   cfg               = [];   cfg               = [];
-  cfg.trials ​       = find(data_func.trialinfo(:,​1) == 1); +  cfg.trials ​       = find(data_lpf.trialinfo(:,​1) == 1); 
-  ​data_sta ​         ​= ft_timelockanalysis(cfg, ​data_func);+  ​timelockSTD ​      = ft_timelockanalysis(cfg, ​data_lpf);
  
-Then, we will apply a baseline correction using ft_timelockbaseline. The five seconds preceding the stimulus will be used as time window for the baseline. ​In+Then, we will apply a baseline correction using **[[/​reference/​ft_timelockbaseline]]**. The five seconds preceding the stimulus will be used as time window for the baseline.
  
   cfg                 = [];   cfg                 = [];
   cfg.baseline ​       = [-5 0];   cfg.baseline ​       = [-5 0];
-  ​data_sta ​           ​= ft_timelockbaseline(cfg, ​data_sta);+  ​timelockSTD ​        = ft_timelockbaseline(cfg, ​timelockSTD);
  
 In the previous steps, you averaged over all standard trials and baseline corrected the result. The same can be done for the deviants. In the previous steps, you averaged over all standard trials and baseline corrected the result. The same can be done for the deviants.
-  + 
-  cfg         = []; +  cfg           ​= []; 
-  cfg.trials = find(data_func.trialinfo(:,​1) == 2); +  cfg.trials ​   = find(data_lpf.trialinfo(:,​1) == 2); 
-  ​data_dev ​     ​= ft_timelockanalysis(cfg, ​data_func);+  ​timelockDEV ​  = ft_timelockanalysis(cfg, ​data_lpf);
   ​   ​
   cfg           = [];   cfg           = [];
-  cfg.baseline = [-5 0]; +  cfg.baseline ​ = [-5 0]; 
-  ​data_dev = ft_timelockbaseline(cfg, ​data_dev);+  ​timelockDEV ​  = ft_timelockbaseline(cfg, ​timelockDEV);
    
-To visualize the data in spatial terms (‘where on the head do we find functional brain activity in response to my different conditions?​’),​ FieldTrip requires information about the spatial layout about the location of the channel on the head. For this tutorial a layout file is provided, which is called ‘nirs_48ch_layout.mat’. In case you would like to get an idea of how to create your own layout file, the following page might be informative:​ [[http://​www.fieldtriptoolbox.org/​tutorial/​layout]]. ​+To visualize the data in spatial terms (‘where on the head do we find functional brain activity in response to my different conditions?​’),​ FieldTrip requires information about the spatial layout about the location of the channel on the head. For this tutorial a layout file is provided, which is called ‘nirs_48ch_layout.mat’. In case you would like to get an idea of how to create your own layout file, the following page might be informative:​ [[/​tutorial/​layout]]. ​
  
 The layout can be loaded using the standard MATLAB function ‘load’. The file ‘nirs_48ch_layout.mat’ contains a structure called ‘lay’. Just for clarity, we will rename the O2Hb channels representing changes in oxygenation concentration ‘functional’. We do this as follows: The layout can be loaded using the standard MATLAB function ‘load’. The file ‘nirs_48ch_layout.mat’ contains a structure called ‘lay’. Just for clarity, we will rename the O2Hb channels representing changes in oxygenation concentration ‘functional’. We do this as follows:
  
-  load(nirs_48ch_layout.mat'​)+  load('nirs_48ch_layout.mat'​)
   label               = lay.label;   label               = lay.label;
   label               = strrep(label,​ '​O2Hb',​ '​functional'​);​   label               = strrep(label,​ '​O2Hb',​ '​functional'​);​
   lay.label ​          = label;   lay.label ​          = label;
    
-There are a number of FieldTrip options available for visualizing the results, such as ft_singleplotER (ER stands for Event Related), which allows you to plot a single channel, and ft_multiplotER,​ which allows you to plot multiple channels on a schematic representation of the head. The ft_multiplotER can also be used in interactive mode to select pieces of the data of interest (for instance specific channels and a specific time window). For now, we are interested to first see whether we find the typical hemodynamic response, hence the changes in oxygenated concentration. Therefore, we select all channels which contain ‘functional’ in their label. This is done by cfg.channel = ‘* [functional]’;​ If you want to see what other options the plotting functions has, you can look at the documentation:​+There are a number of FieldTrip options available for visualizing the results, such as **[[/​reference/​ft_singleplotER]]** (ER stands for Event Related), which allows you to plot a single channel, and **[[/​reference/​ft_multiplotER]]**, which allows you to plot multiple channels on a schematic representation of the head. The **[[/​reference/​ft_multiplotER]]** can also be used in interactive mode to select pieces of the data of interest (for instance specific channels and a specific time window). For now, we are interested to first see whether we find the typical hemodynamic response, hence the changes in oxygenated concentration. Therefore, we select all channels which contain ‘functional’ in their label. This is done by cfg.channel = ‘* [functional]’;​ If you want to see what other options the plotting functions has, you can look at the documentation:​
  
   doc ft_multiplotER   doc ft_multiplotER
  
-Important to remember is that for ft_multiplot ​to run, you need to point FieldTrip to the layout structure using cfg.layout = lay;+Important to remember is that for **[[/​reference/​ft_multiplotER]]** ​to run, you need to point FieldTrip to the layout structure using cfg.layout = lay; 
 + 
 +  cfg                   = []; 
 +  cfg.showlabels ​       = '​yes';​ 
 +  cfg.layout ​           = lay; 
 +  cfg.interactive ​      = '​yes';​ 
 +  cfg.channel ​          = '* [functional]';​ 
 +  cfg.graphcolor ​       = '​r';​ 
 +  ft_multiplotER(cfg,​ timelockDEV);​ 
 + 
 +{{tutorial:​nirs_tut2_multiploter.png?&​800}}
  
-  cfg              = []; +**//Figure 10A so-called multiplot of the data: the average time course displayed per channel.//**
-  cfg.showlabels ​ = '​yes';​ +
-  cfg.layout ​     = lay; +
-  cfg.interactive = '​yes';​ +
-  cfg.channel ​     = '[functional]';​ +
-  cfg.graphcolor ​ = '​r';​ +
-  ft_multiplotER(cfg,​ data_dev);+
  
-{{tutorial:nirs_tut2_multiploter.png?&​400}}+You can also generate a spatial representation of the signal at a certain time point, or averaged over a time window. To plot the response that was found during a specific time window, you will need to specify this by setting limitations to the time dimension. In this case, time is the first dimension, and therefore, the time window can be set by using '​cfg.xlim = [5 7];'. The third dimension here is the strength of the response. We set the scale here from -0.2 to 0.2, but this depends on your datamany fNIRS researchers use block designs, and depending on the block duration, the response may gain a larger amplitude. In the current data, the scale can be derived from the previous figure, which was generated without setting the response scale (zlim in this case).
  
-You can also generate a spatial representation of the signal at a certain time point, or averaged over a time window. To plot the response that was found during a specific time window, you will need to specify this by setting limitations to the time dimension. In this case, time is the first dimension, and therefore, the time window can be set by using '​cfg.xlim = [5 7];'. The third dimension here is the strength of the response. We set the scale here from -0.2 to 0.2, but this depends on your data: many fNIRS researchers use block designs, and depending on the block duration, the response may gain a larger amplitude. In the current data, the scale can be derived from the previous figure. Without setting the response scale (zlim in this case), ​FieldTrip uses the minimum and the maximum in the selected part of the data. Setting the scale manually has the advantage that you can set zero as the middle point in the scale, which can be helpful for the interpretation of the color-coded graph.+<​note>​ Per default ​FieldTrip uses the minimum and the maximum in the selected part of the data for the zlim parameter. Setting the scale manually has the advantage that you can set zero as the middle point in the scale, which can be helpful for the interpretation of the color-coded graph.</​note>​
  
 <​code>​ <​code>​
Line 392: Line 433:
 cfg.layout ​  = lay; cfg.layout ​  = lay;
 cfg.channel ​ = '* [functional]';​ cfg.channel ​ = '* [functional]';​
 +cfg.marker ​  = '​labels';​
 cfg.xlim ​    = [5 7]; cfg.xlim ​    = [5 7];
 cfg.zlim ​    = [-0.2 0.2]; cfg.zlim ​    = [-0.2 0.2];
-ft_topoplotER(cfg, ​data_dev);+ft_topoplotER(cfg, ​timelockDEV);
 title('​[functional]'​);​ title('​[functional]'​);​
 </​code>​ </​code>​
  
-{{tutorial:nirs_tut2_colorcodedtopoplot.png?&​400}}+{{tutorial:nirs_tut2_fig11_topoplot.png?&​400}}
  
 +**//Figure 11; Topographical representation of the measured signal.//**
 =====Summary and conclusion===== =====Summary and conclusion=====
 FIXME  FIXME