NIRS-KIT: a MATLAB toolbox for both resting-state and task fNIRS data analysis

Abstract. Significance: Functional near-infrared spectroscopy (fNIRS) has been widely used to probe human brain function during task state and resting state. However, the existing analysis toolboxes mainly focus on task activation analysis, few software packages can assist resting-state fNIRS studies. Aim: We aimed to provide a versatile and easy-to-use toolbox to perform analysis for both resting state and task fNIRS. Approach: We developed a MATLAB toolbox called NIRS-KIT that works for both resting-state analysis and task activation detection. Results: NIRS-KIT implements common and necessary processing steps for performing fNIRS data analysis, including data preparation, quality control, preprocessing, individual-level analysis, group-level statistics with several popular statistical models, and multiple comparison correction methods, and finally results visualization. For resting-state fNIRS analysis, functional connectivity analysis, graph theory-based network analysis, and amplitude of low-frequency fluctuations analysis are provided. Additionally, NIRS-KIT also supports activation analysis for task fNIRS. Conclusions: NIRS-KIT offers an open source tool for researchers to analyze resting-state and/or task fNIRS data in one suite. It contains several key features: (1) good compatibility, supporting multiple fNIRS recording systems, data formats of NIRS-SPM and Homer2, and the shared near-infrared spectroscopy format data format recommended by the fNIRS society; (2) flexibility, supporting customized preprocessing scripts; (3) ease-to-use, allowing processing fNIRS signals in batch manner with user-friendly graphical user interfaces; and (4) feature-packed data viewing and result visualization. We anticipate that this NIRS-KIT will facilitate the development of the fNIRS field.


Introduction
As a promising non-invasive neuroimaging technique, functional near-infrared spectroscopy (fNIRS) has seen increasing use in neuroscientific studies over the last decades. [1][2][3] Compared to other neuroimaging techniques, fNIRS has several specific advantages: acceptable spatial resolution; lower cost, silent, higher temporal resolution, lower sensitivity to head motion, and fewer physical constraints compared to functional magnetic resonance imaging (fMRI); portability/mobility and few restrictions on use location. These advantages make it an ideal alternative for detecting human brain activity in both healthy people (from newborns to the elderly) and patients with neurological or psychiatric conditions, in both laboratory and natural contexts. Validation studies using simultaneous fNIRS and fMRI scanning have shown fNIRS to be effective for both task-driven activation detection [4][5][6][7][8] and resting-state brain function research. 9,10 FNIRS has been successfully applied to various neuroscience domains, including areas such as motor-related studies with large bodily movements, social cognitive neuroscience, and neurodevelopment. [11][12][13][14][15] In addition to task-evoked activation studies, fNIRS has also seen increasing use in detecting spontaneous brain activity absent of external stimuli during resting state, [16][17][18][19][20] which is particularly useful for several special populations who cannot perform complicated tasks, such as infants and patients with certain deficits.
There are several fNIRS data analysis toolboxes currently available, such as NIRS-SPM, 21 Homer, 22 nirsLAB, 23 FieldTrip, 24 and POTATo. 25 (For more software, see the official website of the Society for fNIRS in Ref. 26.) These useful software packages mainly focus on task-driven cortical activation. However, few software packages can assist for resting-state fNIRS studies. FC-NIRS 27 was recently developed for resting-state fNIRS, supporting functional connectivity (FC) analysis and graph theory-based network analysis. However, it only supports few data sources, which limited its application for resting-state studies severely.
As a fundamental feature of the resting brain, low-frequency fluctuations are closely related to spontaneous neural activity and vary among different brain regions and individuals. Amplitude of low-frequency fluctuation (ALFF) 28 and fractional ALFF (fALFF) 29 quantify the amplitude of these low-frequency oscillations and reflect individual differences or dysfunction. [30][31][32][33] These indices have been supported by several fMRI analysis toolboxes, such as REST, 34 DPARSF, 35 and RESTplus, 36 but have not been supported by any existing fNIRS related toolboxes.
Researchers usually collect data from both task and resting-state modalities in their fNIRS studies to explore the same or related scientific questions. 19,37 For example, researchers have used significantly activated regions in task fNIRS analysis to define seed regions for restingstate FC analysis 38 or used resting-state FC strength to predict task activation. 39 In this situation, researchers must switch between at least two packages to perform task and resting-state fNIRS analyses. Additionally, different processing methods and parameters between packages make it more difficult to integrate, compare, and explain analysis results in a single study. Therefore, there is an urgent need to develop a comprehensive fNIRS toolbox to perform analysis for both modalities.
We developed and present here a MATLAB toolbox called NIRS-KIT that allows researchers to analyze resting-state and/or task fNIRS data in one suite. With a friendly graphical user interface (GUI), this open source package allows researchers to perform common and necessary fNIRS data analysis steps in an easy way. It has good compatibility and flexibility, supporting multiple fNIRS recording systems and customized preprocessing scripts. It can execute most processing steps in a batch processing model, which processes data from multiple participants simultaneously rather than one-by-one. NIRS-KIT can be freely downloaded from the NIH NeuroImaging Tools & Resources Collaboratory website in Ref. 40. In this tutorial paper, we first introduce the main functionality of NIRS-KIT and then provide step-by-step instructions for users. WA, USA) environment. NIRS-KIT has been successfully tested under a variety of operating systems with MATLAB installed, including Windows, Linux, and Mac OS.

General Overview of NIRS-KIT
NIRS-KIT has two main analysis modules: resting-state fNIRS module and task fNIRS module [Figs. 1(b) and 1(c)]. The fNIRS data analysis pipeline implemented in NIRS-KIT is illustrated in Fig. 2. This pipeline consists of common and necessary processing steps for fNIRS data analysis, including data preparation, quality control, preprocessing, individual-level analysis, grouplevel statistics, and results visualization.
For resting-state fNIRS individual-level analysis, FC analysis, ALFF and fALFF, and graph theory-based network analysis to investigate complex topological properties of brain networks (such as local or global efficiency) are supported. In individual-level analysis for task fNIRS, general linear model (GLM) is used to detect task activation.
The following sections will explain these functions and steps in detail, first for the restingstate module and then for specific parts of the task fNIRS module.

Data Preparation
A variety of data are obtained in an fNIRS experiment, and data first need to be imported into the analysis toolkit. The most important data, the raw fNIRS signal time series, are the relative concentration changes of oxyhemoglobin (HbO), deoxyhemoglobin (HbR), or total hemoglobin (HbT). In addition to signal data series, the fNIRS recording's spatial information, consisting mainly of probe and channel locations in two-dimensions (2D) or three-dimensions (3D), are also very important for some analyses, results visualization, data sharing, and publication. Thus we designed a NIRS-KIT data format (stored in MATLAB.mat files), which includes not only the time series signals, but also spatial information about probe set geometry and standard coordinates of source, detector, and channel positions in Montreal Neurological Institute (MNI) coordinates.

Preparation of temporal signals (hemoglobin concentrations)
FNIRS optical time course data are obtained from diverse commercial devices in different types or formats. Some devices only output the raw optical intensity data, whereas others support the calculation of the relative hemoglobin concentration change with in-house conversion functions. Different fNIRS recording systems use different output formats (such as .txt and .csv). These inconsistent data sources pose difficulty for following-up analysis, thus a unified format for hemoglobin concentration is need.
NIRS-KIT provides time course preparation functionality for diverse fNIRS data sources (Fig. 3). If the raw optical intensity data are obtained, such as from Hitachi ETG4000/7000 (.csv), or NIRX (.wl1 and .wl2), they can first be converted into optical density data and then converted to concentration changes of HbO and HbR via the modified Beer-Lambert law. 41 The converted hemoglobin concentration data will be saved in the NIRS-KIT data format. If the converted hemoglobin concentration change data are directly obtained, such as from Hitachi ETG4000/7000 (.csv) or Shimadzu LABNIRS (.txt), they will be reformed and saved in the NIRS-KIT supported format. If the recording system is not one of the above, NIRS-KIT also provides a manual input function to generate the required data format. Here users just need to reorganize their raw data (optical density data or hemoglobin concentration data) into a specific format (.csv) according to the format in the sample file included in the toolbox.
NIRS-KIT has a good compatibility with two widely used fNIRS data analysis packages. Both raw and preprocessed fNIRS time course data from NIRS-SPM (.mat) or HomER2 (.nirs) can be read and saved in the NIRS-KIT supported data format.
Recently, the Shared Near-Infrared Spectroscopy Format (SNIRF), a new fNIRS data format standard, has been proposed by the fNIRS community recently (see Ref. 26), and adopted by Homer3 and FieldTrip. NIRS-KIT also supports SNIRF as input.

Preparation of topospatial information (probe setup)
The probe setup, consisting of the configuration of the sources, detectors, and measurement channels, provides spatial information that needs to be included in the data file. The topographical geometry of the channels is very useful for checking time course signals for problems and result visualization, especially when adopting a complex or irregular probe design. NIRS-KIT provides several standard probe settings in the software package's sample folder (including standard 3 × 3, 3 × 5, 3 × 5 × 2, 3 × 11, and 4 × 4 probes). If an appropriate probe setup file does not already exist, users can use the topomaker module (see Fig. 3) to generate, in a simple and flexible way, a customized probe setup file with arbitrary arrangements of sources and detectors.

Preparation for spatial information in standard brain space
Three-dimensional information in standard brain space, consisting of the spatial locations of optodes and channels, is important for result interpretation, visualization in brain space, and publication so as to enable replication and meta-analysis. NIRS-KIT allows spatial registration of NIRS optodes or channels to standard MNI space with the NFRI toolbox developed by Singh et al. (2005). 42 Two text files for each participant are necessary. One contains fiducial (landmark) coordinates in real space (usually obtained using a 3D digitizer), in which at least four reference positions of the international 10-20 system 43 are given. The reference positions should be spatially separated, and the combination of five standard positions-Iz (inion), Nz (nasion), AL (left preauricular point), AR (right preauricular point), and Cz (central zero)-are commonly used with good effect. The other text file contains the optode and channel coordinates in real space. After loading the above files, a text file containing MNI coordinates of each channel and optode in standard brain space will be outputted. In addition, NIRS-KIT also reports the structural labels of the brain regions closest to the channel locations in three brain atlases (Brodmann, AAL, and LPBA40 atlases).

Data Previewing and Quality Control
After data preparation, the raw fNIRS data should be checked in the time-domain and frequency domain before formal data analysis to examine their characteristics and quality, which help inform choices regarding which participants' data to exclude due to noise and how parameters should be set in subsequent preprocessing. The Data Viewer ( Fig. 4) module of NIRS-KIT provides a variety of visualization functions for plotting individual raw and processed signals, in time or frequency domain.
The Data Viewer module offers abundant and versatile features which include: (1) selection of channels by typing or clicking on the channel in the probe setup; (2) display of both time series and periodogram for selected channels; (3) choosing specific or multiple hemoglobin signal types (HbO, HbR, or HbT); (4) display of signals from multiple channels simultaneously; (5) simultaneously display of raw signals and signals after processing by the preprocessing module (see the following section); and (6) various other features to assist in data examination, such as zoom, pan, data cursor, and change color scheme adjustment. Fig. 4 The Data Viewer module for resting-state fNIRS. A 7-min resting-state fNIRS dataset (n ¼ 9) using two standard 3 × 5 probes (total 44 channels) is shown as an example.

Preprocessing
The fNIRS data typically undergo a series of preprocessing steps prior to statistical analysis, to minimize the influence of noise and artifacts from exogenous and endogenous sources. The NIRS-KIT preprocessing module [ Fig. 5(a)] provides several common preprocessing functions for fNIRS signals, including time point trimming, detrending, motion correction, filtering, noise regression, resampling, and customized processing methods. Researchers can arbitrarily designate the parameters and the order of these steps.

Time point trimming
The first few time points of recordings sometimes need to be discarded due to participants' process of adaptation to the experimental environment. If this option is selected, NIRS-KIT will delete the specified initial and/or last time periods from the data of each participant.

Detrending
FNIRS signals may suffer from systematic increases or decreases over time due to long-term physiological shifts or instrumental instability. NIRS-KIT uses a polynomial regression model to estimate a linear or non-linear trend and subtracts it from the raw hemoglobin concentration signal [as illustrated in Fig. 5(b)]. The default order of the polynomial model is first order to remove linear detrends from the raw time course.

Motion correction
Two parameter-free motion correction methods are available in NIRS-KIT. One is correlation-based signal improvement (CBSI), 44 which mainly relies on the assumption that the neural components of HbO and HbR signals are strongly negatively correlated, whereas motion artifacts are positively correlated. The other new method is temporal derivative distribution repair (TDDR), 45

Filtering
The signal of interest usually occupies a frequency range different from interference, such as machinery and physiological noise (e.g., respiration, heart rate, Mayer wave, and very low-frequency oscillations). NIRS-KIT provides three filtering models (high-pass, low-pass, and bandpass filtering) to remove irrelevant low-frequency and/or high-frequency components [as illustrated in Fig. 5

Noise regression
Undesirable systemic artifacts (especially scalp blood flow) usually contaminate fNIRS functional signals, and how to minimize these sources of interference are still a challenge for fNIRS signal processing. A prominent approach is to use short-distance reference channels (that are sensitive only to signals in the superficial layers of the tissue outside the brain) to record superficial noise, [49][50][51][52][53] then use regression to remove that noise from neural recordings. When shortdistance channels are used, NIRS-KIT provides a noise regression functionality, in which the signals of short-distance channels can be used as the regressor and be removed. In some cases, however, short-distance channels are not available. In these cases, another promising signal separation method proposed by Yamada et al. (2012) 54 can be used to separate functional and systemic signals based on their hemodynamic differences. This method has been incorporated into the preprocessing module to remove the unwanted systemic noise (also see Sec. 4.3.7).

Resampling
For large scale, multi-site, or hyper-scanning studies, it is sometimes necessary to resample the fNIRS time series data recorded with different recording systems (with different sampling rates) into one sampling rate prior to subsequent processing. Also raw signals sometimes need to be downsampled for computational and memory efficiency. For these reasons, users can select the resample option, set the desired sampling rate, and NIRS-KIT will resample the time series data.

Customized processing methods
In fNIRS signal preprocessing, users might want to use their own processing methods which are adapted to specific needs. Also fNIRS signal processing methodology is continuously being improved by researchers around the world so that novel and modified processing methods are periodically published. If these methods are not included or have not been updated in NIRS-KIT, users can write customized processing methods as a MATLAB functions, and easily incorporate them into NIRS-KIT via the customized interface [see Fig. 5

Individual-Level Analysis
After preprocessing, individual-level (or first-level) analysis can be performed. Individual-level FC, ALFF, and fALFF can be calculated here for resting-state studies. In addition, NIRS-KIT also calculates FC matrices as input for subsequent graph theory-based network metrics analysis [ Fig. 6(a)].

Functional connectivity
Resting-state FC measures temporal correlations in spontaneous fluctuations among spatially distributed brain regions 55,56 and has been increasingly used for fNIRS. 17,57,58 NIRS-KIT supports fast ROI2ROI, ROI2Whole-brain, and whole brain (channel-wise) FC analyses.

Network metrics
Network metrics capture complex topological properties of brain networks, such as local and global efficiency 59,60 and modularity. 61 NIRS-KIT allows definition of the network of interest and then calculates the FC matrix of the desired network and saves the result in the file format of Gretna, 62 a widely used MATLAB-based toolbox for calculating graph-theoretic metrics. Prior to network topological analysis, a thresholding procedure is typically applied to exclude the confounding effects of spurious relationships in interregional connectivity matrices. Two commonly applied thresholding techniques are provided in Gretna: the absolute connectivity strength threshold and relative sparsity threshold (or proportional threshold). 9,62-64 For the former, users need to set a threshold value and the connections with weight greater than the given threshold are retained. For the latter, the threshold value is defined as the ratio of the number of actual edges divided by the maximum possible number of edges in a network. Note that the absolute threshold method may lead to different numbers of network edges across different datasets, which may confound between-group comparisons. 65

ALFF and fALFF
ALFF and fALFF, two common resting-state indices characterizing regional spontaneous brain activity, 28,29 were not available in any existing fNIRS analysis software. We implemented this functionality in NIRS-KIT [ Fig. 6(b)]. The time series signal can be approximated by a finite sum of weighted cosine and sine functions whose frequencies (f k ) depend on the sampling frequency (f s ) and the number of collected time points (N): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 1 6 2 where a 0 is the direct-current component [i.e., the mean of the time series xðtÞ]; a k , b k are the weights of the cosine and sine functions; f k ¼ kΔf ¼ kðf s ∕NÞ. The amplitude and phase of k'th frequency component can be expressed as The time series for each channel is first transformed to a frequency domain to generate the amplitude spectrum by Fourier transform (Fig. 7). Then the averaged amplitude across a specific low-frequency range (usually 0.01 to 0.08 Hz) is calculated as the ALFF of the channel [see Fig. 7 and Eq.
where N k is the number of frequency components within the specific frequency range. The fALFF is a ratio of the total amplitude within the specific low-frequency range (0.01 to 0.08 Hz) to the total fluctuation amplitude within the entire frequency range [see Fig. 7 and Eq. (3)]: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 3 1 5 fALFF ¼ X k∶f k ∈½0.01;0.08 For each individual, the standardized ALFF and fALFF (zALFF and zfALFF) of a channel (i) are generated by subtracting the global mean value across all channels and then dividing by the standard deviation across all channels [see Eqs. (4) and (5)] to improve the normality of data distribution across participants, 29,31 which is more compatible to normal-based group-level statistics: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 2 0 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 4 7

Group-Level Statistics
Group-level analyses consist of comparing differences or computing correlations across multiple participants and making population inferences. Several popular statistical parametric models are available in NIRS-KIT for group-level analysis [see Fig. 8 two-sample t-test, paired t-test, correlation analysis, ANOVA (independent or repeated measurement), and average of individual indices. Covariates other than condition (e.g., age, gender, or training time) can be added into these statistical models [ Fig. 8(b)].
To assess statistical significance in the context of multiple tests, the false discovery rate 66,67 and Bonferroni correction to control familywise error rate are supported [ Fig. 8(b)]. Multiple comparison correction can be set to be performed within certain channels of interest by applying a specified mask [ Fig. 8(b)].
Even though several statistical models are provided in NIRS-KIT, several external statistics packages are recommended for more complex statistical analysis, such as the Statistical Package for the Social Sciences (SPSS Inc., Chicago, IL, USA), Statistical Analysis System (SAS Institute Inc., Cary, NC, USA), and R. 68

Results Visualization
Finally, NIRS-KIT provides plentiful visualization functions in both 2D and 3D to visualize individual-level and group-level results. Because of different ways to organize results and different display requirements, the visualization module is introduced in terms of two approaches, channel-wise visualization, and connectivity matrix visualization.

Channel-wise visualization
In ALFF/fALFF and ROI2WholeBrain FC, each channel has a single scalar representing the group-level statistical value or individual-level index. Hence we use channel-wise visualization to display these results. For channel-wise 2D visualization, the probe setup (see Sec. 4.1) is required. A blank canvas with channel locations is generated according to the provided probe geometry. Then color values are mapped with interpolation [ Fig. 9(a)] or without [ Fig. 9(b)]. For ROI2Whole-brain FC, non-interpolated topographical visualization is recommended, and the toolbox will automatically identify the ROI channels and mark them with red circles [ Fig. 9(b)]. For channel-wise 3D visualization, the MNI coordinates of each channel are needed for 3D visualization in standard brain space. In NIRS-KIT, three non-interpolated and one interpolated channel-wise 3D visualizations are provided. The first non-interpolated 3D visualization is mapping the channel-wise statistical values on to the standard brain space using the nfri_m-ni_plot function in the NFRI toolbox. 42,69,70 The second is to project the analysis result of each measurement channel directly on to a standard surface template (e.g., ICBM152), set the display parameters (such as, the sphere size and colorbar limits), and then output the resulting figure. The last, more flexible approach, is to output a brain image in the widely used NIFTI format. Then users can load and plot this image with several other brain visualization toolboxes, such as MRIcroGL 71 and Surfice. 72 For interpolated channel-wise 3D visualization, the EasyTopo toolbox provided by Tian et al. (2013) 73 is used. Here the stereotaxic coordinates of the brain surface and fNIRS channels are converted into spherical coordinates, where 2D angular interpolation of the channel-wise data is implemented to obtain a topographic image of brain activation in the latitude-longitude space. Then the interpolated image is projected back onto the original brain surface. The first non-interpolated 3D visualization using the NFRI toolbox was illustrated in Fig. 9(c). The other three visualization types will be shown in Sec. 5.6 for displaying the task fNIRS analysis results.

Connectivity matrix visualization
The result of whole brain FC is saved as an n × n matrix (n is the number of channels), where each value in the matrix represents the FC between a paired channels or group-level statistic's value of each FC. For 2D visualization of whole brain FC matrix, users can choose to only  7  8  11  12  13  15  16  17  20  21  22  1  2  3  4  5  6  9  10  14  18  19  28  29  32  33  34  37  38  39  41  42  43  23  24  25  26  27  30  31  35  36 Fig. 9 Examples of visualization for resting-state analysis results: (a) group-level zfALFF map (interpolation mode); (b) ROI2Whole FC map (non-interpolation mode, and channel 12 is the seed channel); (c) non-interpolated channel-wise 3D visualization for group-level zfALFF result using the NFRI toolbox; (d) 2D whole brain FC matrix map; and (e) 3D visualization of whole brain connectivity matrix by BrainNet Viewer using the node and edge files generated by NIRS-KIT. Note: All above color bars represent group-level (n ¼ 9) one-sample t -test statistical values, and the color axis limits are set for illustration only. In A and B panels, only the left hemisphere probe (channels 1 to 22) is displayed to reduce clutter. display the significant edges by setting a threshold p-value, define the displayed subnetworks, and arbitrarily reorder the channels [ Fig. 9(d)]. For 3D visualization of whole brain FC matrix, an Excel (Microsoft Corp., Redmond, WA, USA) file containing MNI coordinates of each channel (see sample files in the package) is need to be loaded, then NIRS-KIT will output a node file and an edge file that can be loaded and visualized by the most existing brain connectome visualization toolkits [ Fig. 9(e)], such as BrainNet Viewer. 74 Several visualization options are available, such as setting display modes, color bar limits, and p-value thresholds.

Data Preparation
There is no difference between task fNIRS data preparation and resting-state fNIRS data preparation (see Sec. 4.1).

Data Previewing and Quality Control
A specific panel is provided in the Data Viewer module for task fNIRS (Fig. 10). This panel allows stacked display of time series and periodograms of reference signals from the experiment's protocol. The task-related information (stimuli onset and duration) can be defined by manual input or loading a corresponding task design information file. Users can select singleor multiple-experimental conditions/task types to be displayed.

Preprocessing
Task fNIRS data preprocessing is similar to that of resting-state fNIRS (see Sec. 4.3). For taskbased data preprocessing, the frequency components related to the task protocol should not be removed, regardless of which filter is used.

Individual-Level Analyses
After preprocessing for task fNIRS data, individual-level analysis to detect channel-wise taskevoked neural activation can be done via the widely used mass univariate statistical method based on GLMs. 75 The routine individual-level statistical analysis of task fNIRS data comprises of the following steps: (1) specification of a GLM, which considers the observed hemodynamic signal (dependent variable) as a linear combination of regressors of interest (task variables), nuisance covariates (such as the superficial noise measured by short-distance channels), and an error term. For GLM specification, the canonical hemodynamic response function implemented in SPM is used to construct the reference time series representation from task variables; (2) estimation of GLM parameters on a channel-by-channel basis, which finds the activation beta value (weight coefficient in the linear model) for each experimental condition; and (3) calculation of the condition-wise effects of interest using user-defined contrast vectors as the input for subsequent group-level inference.
NIRS-KIT supports all the above steps in an easy-to-use, batch processing manner in the Task Individual Analysis module (Fig. 11). There are two ways to input the experimental design information for model specification. If the onset time and duration of each event/block are the same across all participants, the task design information can be entered manually [ Fig. 11(a)]. If the design information varies across participants, the design matrix can be constructed by loading a prepared .mat file containing the required experimental design information (names of experimental conditions, onset times of each condition, and event or epoch durations) for each participant [ Fig. 11(b)]. Fig. 11 Individual-level analysis interfaces for task fNIRS. (a) The individual-level analysis interface for task fNIRS with manual input of task design (when design is the same for all participants). (b) The individual-level analysis interface for task fNIRS with different designs per participants; one column of randomly generated covariates is added into the GLM as an illustration.

Group-Level Statistics
There is no difference with resting-state fNIRS group-level statistical analysis (see Sec. 4.5).

Results Visualization
For task fNIRS activation results, each channel has a single-value representing the group-level statistic's value, thus channel-wise result visualization should be used for visualization (see Sec. 4.6). Here group-level statistical results from a finger tapping task fNIRS study (n ¼ 9) with two standard 3 × 5 probes are shown for illustration. One 2D [ Fig. 12

Discussion
We developed a comprehensive and versatile MATLAB toolbox called NIRS-KIT, suitable for both resting-state and task fNIRS data analysis. It covers common and necessary fNIRS data analysis steps, including data preparation, quality control, preprocessing, individual-level and group-level analyses, as well as results visualization, in a GUI-based, easy-to-use, batch-processing-supported software suite. We believe NIRS-KIT provides researchers an alternative toolbox for performing complex fNIRS data analysis and will facilitate the advancement of human brain function research.
NIRS-KIT is an integrated platform that supports analysis for both resting-state and task fNIRS data. In recent years, in addition to task-evoked activation studies, fNIRS has also been increasingly used to detect the spontaneous brain activity pattern in resting state without external stimuli. However, few packages support comprehensive analysis for resting-state fNIRS.  NIRS-KIT provides related researchers as powerful analysis to perform resting-state fNIRS data analysis. On the other hand, sometimes researchers will collect both task and resting-state fNIRS data in a single study. Prior to NIRS-KIT, the existing fNIRS analysis toolboxes handled one type of experiment or the other. Using NIRS-KIT, researchers can complete fNIRS data analysis of both modalities in one software suite. NIRS-KIT has good compatibility with multiple data sources and flexibility in usage. First, it supports a variety of the raw data obtained directly from several widely used fNIRS devices, and the shared data format (SNIRF) proposed by the fNIRS community. In addition, it also supports the raw or processed data output from two widely used analysis software (NIRS-SPM and Homer2) as the input. If the data format is not yet support, users can manually reformat the data to a relativity simple text format supported by NIRS-KIT in an easily way. Furthermore, this platform also allows add-ons to incorporate customized processing methods into the preprocessing procedure. Additionally, it supports batch data analysis, which decreases repetitive operations, time cost, and the probability of mistakes, especially when dealing with a large sample with complex processing protocol.
It has been noted that there are still no standard processing procedures in the field of fNIRS, and high heterogeneity of signal processing methods and parameter settings may lead to biased results and undermine the reproducibility of and comparability between studies. 76,77 In NIRS-KIT, several mainstream and validated preprocessing methods 53,78 are provided, which has the potential to reduce the impacts of heterogeneity to a certain extent. However, more efforts should be taken by researchers around the world to build standardized processing procedures, and processing methods and parameters should also be reported in more detail in publications.
Since NIRS-KIT development is at an early state, improvements are still required. For example, an alternative to FC, effective connectivity mainly depicts the causal influence that one neural system exerts over another. 79 Effective connectivity has been applied in several fNIRS studies. [80][81][82] In addition, although multiple popular group-level statistical models are provided, individual-level statistical tests are not supported in the current version of NIRS-KIT. We plan to incorporate these functionalities to future versions of NIRS-KIT. Furthermore, some of our research group's on-going work will be added to NIRS-KIT, e.g., blind source separation (BSS) for separating unrelated physiological noise (such as scalp blood flow) from neural related signals. 37,83,84 Six BSS algorithms (PCA, 85 SOBI, 86 JADE, 87 Fast-ICA, 88 Infomax-ICA, 89 and ERBM-ICA 90 ) are planned for inclusion.

Disclosures
The authors declare no conflicts of interest.