Reduction of global interference in functional multidistance near-infrared spectroscopy using empirical mode decomposition and recursive least squares : a Monte Carlo study

Y. Zhang zyhit@hit.edu.cn Dept. Automatic Meas. & Cont., Harbin Institute of Technology, No.92, West Da-zhi Street, Nangang District, Harbin 150001, China J. Sun Dept. Automatic Meas. & Cont., Harbin Institute of Technology, No.92, West Da-zhi Street, Nangang District, Harbin 150001, China P. Rolfe peterrolfe@aol.com Dept. Automatic Meas. & Cont., Harbin Institute of Technology, No.92, West Da-zhi Street, Nangang District, Harbin 150001, China


INTRODUCTION
Near-infrared spectroscopy (NIRS) has been developed to allow noninvasive measurement of optically absorbing and scattering molecules in medicine and biology [1].Although NIRS may be useful in a variety of tissues and organs there has been a major interest in its use for studies of the brain and it has been applied in adults [2], the fetus [3] and the newborn [4].Of particular importance has been the ability of the method to monitor changes in cerebral oxyhaemoglobin and deoxyhaemoglobin (HbO 2 , HHb), total haemoglobin (Hb tot ) and the redox state of the respiratory enzyme cytochrome oxidase (cyt aa3 ) [5,6].The quantitative measurement of the target chromophores in terms of absolute concentrations has always been a challenge and this has directed NIRS research, over several decades, particularly with respect to instrumentation and in the underlying theoretical considerations.This has led to the use of continuous wave, time-resolved, intensity modulated, and spatially-resolved instruments [7] based on theories of photon propagation through tissue described by the Radiative Transport Equation, the diffusion approximation and diffuse reflectance [8].
The use of stimuli to evoke physiological responses has been one particular area of relevance to NIRS, leading to what has become the well-recognised method of functional nearinfrared spectroscopy (fNIRS).This has been demonstrated to be able to determine cerebral concentration changes of HbO 2 and HHb during functional activation of the cerebral cortex [9].fNIRS may be compared with other techniques, such as functional magnetic resonance imaging (fMRI), magnetoencephalography (MEG), positron emission tomography (PET), and electroencephalography (EEG) [10,11].It does appear to have several advantages over these other methods, such as portability, fewer physical restrictions and greater practicality, good temporal resolution, safety, and inexpensive instrumentation.In fNIRS studies, measurement of the concentration changes of HbO 2 and HHb, the two primary absorbing chromophores in the brain tissue that vary dynamically with a functional task, is achieved.This then provides a useful description of the cerebral haemodynamics.However, there are problems in using fNIRS due to the presence of physiological interference, which is mainly from perturbations caused by cardiac events, such as changes in blood flow, blood pressure and blood volume following cardiac contractions, as well as from breathing [5].All of these forms of interference are located both in the vasculature of the superficial layer of the brain and inside the brain, and are generally termed "global interference" or "systemic physiological interference".Severe contamination of global interference remains a serious problem for fNIRS interpretation and analysis.This has meant that without appropriate global interference reduction the full potential of fNIRS has not yet been realised.
Previous attempts have been made to suppress global interference and thereby improve the brain activity measurements.The use of low pass filtering techniques to remove the interference caused by cardiac oscillations has been reported [12].However, low pass filtering is not effective at removing the specific physiological noise signals such as respiratory and blood pressure variation since these fluctuations are difficult to be distinguished from the haemodynamic response to brain activity by frequency characteristics alone.Saager et al. [13] utilized a dual-detector system and least squares method to remove top-layer-only fluctuations and validated the effect of the methodology by performing Monte Carlo simulations based on a two-layer turbid media model.Zhang et al. [11] proposed a spatial eigenfiltering algorithm to separate an activity-evoked response from systemic physiological interference in diffuse optical imaging data.They assumed that the first several spatial eigenvectors of the baseline concentration correlation matrices were dominated by interference patterns and hold most of the total energy.Morren et al. [14] and Zhang et al. [15] adopted adaptive filtering to remove global interference .The difference between the two approaches is that the former used signals from a pulse oximeter as a reference and the latter used signals from short interoptode distance as a reference.A Kalman filtering model has also been used to analyze interference components [16,17].Abdelnour et al. proposed an adaptive general linear model based on a Kalman filtering algorithm for real-time assessment of brain function [18].Despite the fact that the general linear model in their papers exhibited the potential for removal of physiological signals from cardiac, respiratory, and Mayer wave fluctuations, the authors expressly indicated that prior distributions must be assumed for both the process and observed noise in the Kalman filtering model.
Empirical mode decomposition (EMD) and Hilbert spectral analysis (HSA) have been shown to be able to separate, identify and remove interference arising from some cardiac events and breathing [19].This technique has the advantage of simplicity in instrument design because it only needs one source and one detector in the NIRS probe.However, the EMD-HSA technique may not be appropriate for removing other physiological fluctuations such as blood pressure and other lower frequency variations, although these fluctuations in superficial layers are highly correlated with those in deep layers.
In our study presented here we adopt the spatially-resolved or multidistance measurement method and a theoretical analysis of global interference reduction based on EMD and the least squares criterion.The short-distance fNIRS measurement is treated as the reference channel comprising of superficial haemodynamic changes induced by physiological fluctuations and the long-distance fNIRS measurement is treated as the signal channel containing both the functional haemodynamic response and global interference.We aim to remove global interference that is correlated with superficial haemodynamic fluctuations, evoked by cardiac contractions, breathing, blood pressure, etc.By decomposing the short-distance measurement with the EMD algorithm, we separated the in-terference into different intrinsic mode functions (IMFs) possessing distinct frequency characteristics.The least squares criterion was then used to adjust the corresponding weighting coefficients to estimate global interference with the obtained IMFs.To accelerate the computation we adopted the recursive least squares (RLS) to decrease the computation complexity.Monte Carlo simulations of a five-layer human model were used to investigate the performance of the EMD-RLS for removing global interference in brain activity measurement.

Empirical Mode Decomposition
EMD is an adaptive signal processing technique, which aims at extracting all the fundamental modes embedded in a signal without any requirement of stationarity or linearity of the data [20].EMD is based on an iterative algorithm in which the highest frequency oscillation is removed from the data with each repetition.Following each iteration a residue remains containing lower frequency information.The process is repeated until only a trend remains, leading to a local, adaptive decomposition, in which intrinsic oscillations within the data may be identified.
EMD is derived from the simple assumption that any signal consists of different coexisting modes of oscillation with a separated time-scale.Each of these oscillatory modes is represented by an intrinsic mode function (IMF) and the IMF must satisfy two criteria.Firstly, the number of extrema and the number of zero-crossings must either equal or differ at most by one and, secondly, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima must be zero.Thus, for a given signal x(k) the appropriate EMD algorithm can be achieved according to these criteria.The upper envelope is created by cubic spline interpolation of the local maxima and the lower envelope is created by interpolation of the local minima.The first component can be obtained by subtracting the mean m 1 (k) calculated with the upper and lower envelopes from the original series x(k).
The resulting component c 1 (k) is the first IMF if the mode fulfils the definition of an IMF.If c 1 (k) is not an IMF, the sifting process has to be repeated as many times as is required to reduce the extracted signal to an IMF.Once this condition is reached, the first component c 1 (k) is separated from the original series x(k) to leave the first residue r 1 (k): The residue r 1 (k) now contains information about the components for a longer period, and therefore it is treated as the new data and is resifted to obtain additional components.The procedure will be continued until finally it meets a stopping criterion.By this procedure, lower order IMFs capture fast oscillation modes while higher order IMFs typically disclose slow oscillation modes.The original signal x(k) can thus be expressed as follows: where N is the number of IMFs, and IMFs c j (k) (j = 1, . . ., N) are nearly orthogonal to each other and all have zero means, r N (k) is the final residue.The result of the EMD produces N fundamental components and a residue signal.The residue itself is a very low frequency component or a trend, and can also be regarded as the last IMF.

Multidistance Measurement and EMD-RLS Algorithm to Remove Global Interference
A five-layer model of the adult human head, consisting of scalp, skull, cerebrospinal fluid (CSF) , gray matter, and white matter was used here (see Figure 1).The use of particular source-detector separation allows us to distinguish between processes occurring at different tissue depths, as indicated by the "banana-shaped" regions that encompass statistically determined photon paths [21].The true haemodynamic changes related to the specific stimuli present in the gray matter layer.Thus we use short interoptode distance with S-D1 to probe the superficial tissue layer, and long interoptode distance with S-D2 to probe the deep tissue layer.In this approach, the near detector is positioned at 5 mm to 15 mm from the source and the far detector is positioned at more than 30 mm from the source [22].In our study here, the short interoptode distance is chosen at 10 mm and the far interoptode distance is chosen at 40 mm.To remove global interference in fNIRS measurement, short-distance measurements were therefore used as the reference channel to estimate global interference present with the long-distance measurement according to the following analysis.∆A(λ 1 ) where ε HbO 2 (λ) and ε HHb (λ) are the molar extinction coefficients of HbO 2 and HHb, respectively.Hereinto, L is denoted as: where DPF(λ) is the differential path length factor at a wavelength λ, and r is the linear distance from the source to the detector.The concentration changes of oxyhaemoglobin, ∆[HbO 2 ], and deoxyhaemoglobin, ∆[HHb], can be derived from Equation (4): where L −1 denotes the inverse matrix of L. When we have measured the change of optical density, the time series of ∆[HbO 2 ] and ∆[HHb] obtained with S-D1 and S-D2 can be calculated using Equation (6).The objective is to subtract the estimated global interference from the long-distance measurement.EMD was used here to decompose short-distance measurement into several IMFs according to the frequency characteristic.
where x(k) is the ∆[HbO 2 ] or ∆[HHb] acquired from the shortdistance measurement at instant k, N is the number of IMF.Note that for convenience the residue term is absorbed in the summation as the last term.Then we consider the matter of mapping the IMFs of x(k) in a fashion that approximates global interference i(k), that is related to x(k).The IMFs can be combined utilizing various possible methodologies under various objective functions designed to match i(k).Mathematically, the linear fit is adopted here, as given by Equation ( 8): where î(k) is the estimation of i(k) and the coefficient w i (k) is the weighting assigned to the ith IMF at instant k.
The haemodynamic changes acquired with S-D2 is modelled as a sum of the true brain activity response y(k) and global interference i(k).The true brain activity response can then be estimated as e(k): where s(k) is the ∆[HbO 2 ] or ∆[HHb] acquired from the longdistance measurement at instant k.To produce an estimation of global interference in the brain activity analysis, the following cost function of weighted least squares is chosen for the optimization criterion: The parameter 0 < χ ≤ 1 is the forgetting factor that controls the memory span of the algorithm.Using Equation (8) and Equation ( 9), the following equation can be obtained from Equation (10): Here, the weight w i (k) was used to adjust the proportion of each specific IMF in global interference.The optimal coefficient w i (k) for minimizing J k can be obtained by differentiating Equation (11) with respect to w i (k) and setting the derivatives to zero.This yields the following equation: by defining Then Equation ( 13) can be rewritten as Considering all the values of the index i and j, we can obtain the following equation with the matrix form.
This can be briefly denoted as The optimal coefficients can be acquired by the following equation.
Since the straightforward computation of the inverse of R −1 (k) results in an algorithm with computational complexity of O[N 3 ], we implement the computation of the inverse matrix by means of the Matrix Inversion Lemma [23]. where We define the a prior error as By expressing s(k) as a function of the a prior error [23] and the weight vector w(k) can be expressed as From the above analysis, e(k) is the estimated value of the true brain activity response y(k).The critical aspect of using EMD-RLS here is to optimize the estimation of global interference utilizing linear mapping between the decomposed IMFs and the long-distance measurement.

Method
In the present study we use one light source (S) and two light detectors (D1 and D2) shown in Figure 1.One light source with two wavelengths, 750 and 830 nm, and two detectors were positioned on the surface of the medium to collect diffuse reflectance data.The light source and the detectors are all positioned perpendicular to the surface of the medium.
A five-layer slab model consisting of scalp (sc), skull (sk), cerebrospinal fluid (CSF), gray matter (gm), and white matter (wm) was used to represent a human adult head.

Monte Carlo Simulations
To determine the effectiveness of our algorithm, Monte Carlo simulations of the five-layered human head model were performed.The simulations were carried out using the typical fNIRS probe arrangement illustrated in Figure 1.
The Monte Carlo code used here is an extension of the general multi-layer, three dimensional, weighted photon Monte Carlo codes developed by Wang et al [24].Scattering anisotropies were assumed to be 0.9 and Fresnel reflection at the tissueair boundary was also considered.We assumed the same refractive index of n = 1.4 for all layers [25].The standard parameters used in the simulation are given in Table 1.Thickness, transport scattering coefficient µ s , HbO 2 and HHb baseline concentrations were taken from published data [26][27][28].
The baseline concentration of HbO 2 and HHb assumed that the oxygen saturation in the head was 70%.Absorption coefficients were calculated with the HbO 2 and HHb baseline concentrations and the molar extinction coefficients.The molar extinction coefficients at 750 and 830 nm were obtained from the literature [29]. Baseline

Generation of Time Series of Haemodynamic Changes and Physiological Interference
The haemodynamic changes were simulated as a combination of the functional haemodynamic responses and the physiological interference.The functional haemodynamic responses in gray matter were defined as the convolution of the stimulation s(t), [s(t) = 0 for rest period, and 1 for stimulation] and a prototypical haemodynamic impulse h(t) [15,26,30].
The parameters b and d in the gamma variate function were set at 8.6 and 0.56 respectively, which corresponds to recent findings [26,30].
The evoked haemodynamic response u(t) was the convolution of s(t) and h(t): where K is a simple scaling factor in our model.In order to make the simulation as realistic as possible we generated physiological interference by a combination of cardiac fluctuation c(t),respiratory fluctuations r(t), low frequency oscillation m(t), very low frequency oscillation v(t),and independent fluctuation τ(t) induced by the temperature changes and the sweat on the skin.The haemodynamic changes in each layer are denoted as follows: C 2,3,4,5 HbO 2 (t) = HbO 2,3,4,5 2base + α C 2,3,4,5 HHb (t) = HHb 2,3,4,5 base + α 2,3,4,5 HHb c(t) + β 2,3,4,5 HHb r(t) where C 1 HbO 2 (t), C 1 HHb (t), C 2,3,4,5 HbO 2 (t), and C 2,3,4,5 HHb (t) represent the concentration of HbO 2 and HHb in each layer as a function of time, with the superscripts 1 to 5 indicating the layer index for scalp, skull, CSF, gray and white matters, respectively.HbO 2base and HHb base represent average or baseline concentrations.The cardiac fluctuation c(t), respiratory fluctuation r(t), low frequency oscillation m(t), and the very low frequency oscillation v(t) were all simulated as the summation of sinusoidal wave and white Gaussian noise.The form of c(t), r(t), m(t), and v(t) is defined as c(t) = sin(2π The frequency values of each sine wave were defined according to the mean frequencies of real signals ( f c = 1 Hz, f r = 0.25 Hz, f m = 0.1 Hz, f v = 0.04 Hz) [5,[31][32][33].The coefficients α, β, γ, ζ and ι with layer index as superscript and HHb or HbO 2 as subscript are the haemodynamic variation amplitude control parameters, which can be found in Table 2.The independent interference τ(t) is generated by biased and low pass filtered Gaussian white noise.The independent interference τ(t) in the scalp layer and the different weight on cardiac, respiratory, and low-frequency and very low frequency fluctuations in each layer were to simulate a certain amount of uncorrelated changes in the superficial layers compared with deep layers.The parameters used are based on Scholkmann et al. [33] and Zhang et al [15].
The simulated haemodynamic changes were used to calculate the optical measurement by Monte Carlo method.The sampling rate was set to 10 Hz and the whole time series for the changes of optical density were acquired under the assumption that the scattering properties of the head do not vary with time.The experiment is designed as a 5-epoch block and each individual epoch consisted of a series of 400 points, 200 points of rest and 200 points of stimulation.

Intrinsic Mode Functions of Haemodynamic Changes with Short Interoptode Distance
Partial pathlength (PPL) is the mean path length of the photons travelling in a specified region, which could reflect the effective transmitting path for the photons in each layer.According to the MLBL, the partial pathlength is the key factor that determines the light attenuation under the assumption of constant absorption and scattering coefficients.tions for the short interoptode distance, the partial pathlength in the gray matter (PPL gm ) is very small, the value of PPL gm being about 10 −2.6 × that in the superficial layer.Thus, the optical measurement made at a short interoptode distance could be very appropriate to be used to estimate global interference, where the superficial variation or non-activation signals dominate the measurement.
The simulated optical measurements at 750 and 830 nm were acquired with Monte Carlo simulations.The changes of optical density for the 10 mm source-detector distance are shown in Figure 3.The shaded regions in the figure indicate the periods of evoked stimulation.The changes of optical density should be zero under ideal conditions if no interference is present.As shown in Figure 3, the fluctuations are evident in the changes of optical density for both 750 and 830 nm.
The concentration changes of HbO 2 and HHb (∆[HbO 2 ] and ∆[HHb]) were calculated with the MLBL, which has been described in Equation ( 6).The calculated ∆[HbO 2 ] and ∆[HHb] contain several components of interference, including those induced by events related to the cardiac cycle and to breathing, as well as by spontaneous physiological low-frequency oscillations, and very low frequency oscillations.The EMD method was utilized to produce the IMF components for ∆[HbO 2 ] and ∆[HHb].As shown in Figure 4 and Figure 5, five IMF components and a residue were generated with EMD for ∆[HbO 2 ], and five IMF components and a residue were generated with EMD for ∆ [HHb].
It is clearly demonstrated in Figure 4 and Figure 5 that EMD can successfully decompose the signals of interest into components whose frequency relationship is fully consistent with the known sources of interference.The cardiac pulse is clearly evident in the first IMF, breathing is found in the second IMF, the low frequency oscillation is found in the third IMF, and the very low frequency oscillation is found in the fourth IMF.The additional IMF components presented in the decomposition process are mainly due to the introduction of a slowly varying random signal as the independent interference, which was simulated as the uncorrected superficial response compared with the intracerebral response.

Application of EMD-RLS for Global Interference Reduction
After the IMFs are acquired, the data analysis described in Section 2.2 is applied.From the results as shown in Figure 4 FIG. 5 The IMF components of concentration changes of HHb acquired with shortdistance optode.The units of the y-axis for IMFs are micromole (¯M) and Figure 5, the number of IMFs N = 6.The forgetting factor χ in our study was set to 0.99999, which is relatively stable in our tests.The total signal length is 3000 samples and the first 1000 samples are used to train the signal and obtain the weighting coefficients assigned to the specified component (see Table 3).After the optimal coefficients are determined, the IMFs are applied to approximate global interference.The remaining 2000 samples are then applied for interference cancellation by subtracting the estimated global interference.As seen in Table 3, the coefficients for each component are different.It can be seen that the magnitudes of the first four coefficients are larger than those of the last two coefficients.This can be interpreted that the first four IMFs are the interference components corresponding to heartbeat, breathing, low frequency oscillation and very low frequency oscillation, which are correlated with global interference.
Concentration changes of HbO 2 during the simulated measurement are shown in Figure 6.Thereinto, Figure 6a shows the ∆[HbO 2 ] time series results from a 10 mm source-detector distance that were calculated with the MLBL and Figure 6b shows the block average result.Similarly, Figure 6c and 6d show the ∆[HbO 2 ] time series and block average results from a 40 mm source-detector distance, again calculated with the MLBL.However, neither the results from Figure 6c nor 6d can clearly reveal the true activation signal, since it is obvious that the raw time series has been contaminated by global interference.After signal processing with the EMD-RLS, as seen in Figure 6e and 6f, the evoked haemodynamic changes are clear  and most of global interference has been removed.However, the magnitude of ∆[HbO 2 ] is underestimated when calculating with the MLBL because the partial path length in the gray matter is always smaller than the differential path length.This is generally referred to as the partial volume effect (PVE) [34].By means of Monte Carlo simulations, the ratio of the differential path length to the partial path length in the gray matter can be achieved to compensate for the PVE effect.After we have compensated for the PVE effect in the recovered results we can compare the recovered results with the real evoked haemodynamic changes in the gray matter.This comparison can be seen in Figure 6g and 6h.In these two figures, the solid line denotes the recovered results with PVE compensation and the dashed line denotes the true evoked haemodynamic changes used in the simulation.Although some fluctuations still remain, the recovered ∆[HbO 2 ] processed with the EMD-RLS algorithm provide an obvious evoked response.By calculating ensemble average for the whole time series, the results demonstrate that the proposed methodology could remove approximately 90% of global interference.IMF order i 1 2 3 4 5 6 w i for HbO 2 0.5437 0.8034 0.8130 0.8384 0.0541 0.0685 w i for HHb 0.1756 0.7280 0.8723 1.1101 0.0394 0.0511 TABLE 3 Optimal weighting coefficients assigned to the interference components.was recovered.We have also run the simulations for 5 mm and 45 mm short and long spacing respectively and the results are just as clear.

Discussion
Global interference can greatly degrade the performance of fNIRS measurement of evoked brain activity response.An approach to help identify and separate the interference components in fNIRS is to include auxiliary physiological measurements and to analyze the signal with an adaptive filtering algorithm [14,17].Many instruments such as the pulse oximeter, electrocardiogram (ECG), chest band respirometer, spirometer and capnograph can be used to achieve auxiliary physiological measurements.This method is effective in reducing global interference, but the indispensability of additional equipments is the limitation of its application.Prince et al. [16] and Adbelnour et al. [18] utilized sine/cosine terms and Kalman filtering to model the specific physiological noise.However, this method needs to assume a prior frequencies (e.g. the cardiac, ∼1 Hz; respiratory, ∼0.25 Hz; and Mayer wave frequency, ∼0.1 Hz) for the sine/cosine terms.In our previous work, the methodology based on EMD-HSA has also been presented good performance in removing cardiac and respiratory interference [19].However, this methodology is not appropriate for the interference induced by low fre-quency oscillations, very low frequency oscillations, etc.Thus, a novel method based on EMD-RLS algorithm and multidistance probe configuration was proposed in this paper.
In this paper, EMD algorithm adaptively decomposes the superficial haemodynamic changes into several interference components and RLS algorithm dynamically adjusts the weighting of the corresponding interference to estimate global interference.Although EMD has been widely used to decompose complex signals and the RLS method has been broadly adopted for denoising in many domains, the approach of combination of EMD and RLS has not been applied for denoising in fNIRS brain activity measurement.
We An important matter that should be considered here is the traditional EMD algorithm adopts the cubic spline interpolation as an effective tool processing non-stationary signal and several issues may occur during the decomposition process.
One of the issues is the presence of overshoot and undershoot found during the decomposition of the second IMF in ∆ [HHb].In order to solve the problem, many researchers attempted to overcome the overshoot [35], for example Qin et al. [36] used a Segment Slide Theory and Fan et al. [37] presented the piecewise linear fractal interpolation as the spline interpolating.Another issue is the modal distortion and mixture [38,39] that often influence the performance of the EMD.Although this issue is not observed in our simulation study, the resolution of this issue may be useful in our future study for the EMD-RLS algorithm.The effective EMD methods based on independent component analysis [38] and the masking signal [39] were proposed to overcome this problem.Thus, some deeper investigations of EMD should be able to improve the performance of the EMD-RLS further.
One of the additional issues in fNIRS is that there can be a further source of interference, which is observed synchronously with a heart rate increase during the specific task (e.g.finger opposition task) [40,41].This kind of interference is not independent of the functional response itself and thus using short-distance measurement as the reference channel may over-correct the target measurements.For these cases, methods considering the spatial locality difference between the brain functional response and global interference have been proposed based on principle component analysis [11,41] or independent component analysis [10].The EMD-RLS algorithm could still be applicable for these cases if we put the target channel on the activated area and the reference channel on the unactivated area.
A truly rigorous evaluation of this method in vivo requires an uncontaminated evoked brain activity response signal, which, unfortunately, is unavailable.Meanwhile, the PVE effect in vivo can not be exactly compensated and the quantitative comparison of the recovery response and the true response of brain activity is difficult.Thus, Monte Carlo simulation is implemented for quantitative analysis here as the preliminary study.As for the test of this method in vivo with practical NIRS measurement, other signal processing methods such as power spectrum density will be utilized as auxiliary evaluation tool.This is the subject of our on-going research and will be published in due course.

Conclusions
The empirical mode decomposition is a promising tool for analyzing fNIRS signals through the frequency characteristic of specified IMFs.However, the optimality for the reconstruction of the given signal for IMFs is not considered.Thus, EMD-RLS algorithm based on multidistance probe configuration is used to study global interference reduction in fNIRS brain activity measurement.The evaluation of this method is carried out through simulating the evoked haemodynamic response and physiological interference.The method is effective to remove global interference induced not only by heartbeat and respiration but also by low frequency oscillation and very low frequency oscillation, and other correlated interference between superficial and deep layer.The advantage of EMD-RLS in multidistance fNIRS measurement compared with other possible methods also arises from its convenient implementation; it neither requires an auxiliary measurement instrument nor the dependence on a prior knowledge of the global interference frequency.Thus, the methodology has clear potential for use in fNIRS brain activity measurement.

FIG. 1
FIG. 1 Schematic illustration of five-layered slab human head model and multidistance optode configuration.S is the photon source, and D1 and D2 are detectors.

Figure 2
Figure2shows the block diagram of the reduction of global interference for fNIRS brain activity measurement in our study.The objective is to subtract the estimated global interference from the long-distance measurement.EMD was used here to decompose short-distance measurement into several IMFs according to the frequency characteristic.

Fig. 2 FIG. 3
Fig.2The simulated measurements of optical density for th

Fig. 3 FIG. 4
Fig.3 IMF components of concentration change of oxyhaemoglobin acquired with short-distance optode.The units of the y-axis for IMFs are micromole ( M).

Fig. 4
Fig.4 IMF components of concentration change of deoxyhaemoglobin acquired with short-distance optode.The units of the y-axis for IMFs are micromole ( M).

Fig. 5
Fig.5 RLS adaptive filtering to remove physiological interference.(a) The time series of HbO2 calculated from S FIG. 6 EMD-RLS to remove global interference.(a) The time series of concentration changes of HbO 2 calculated from S-D1 with 10 mm source-detector separation, and (b) its block averaged result.(c) and (d) The time series of concentration changes of HbO 2 calculated from S-D2 with 40 mm source-detector separation.(e) and (f) EMD-RLS result for the target measurement.(g) and (h) EMD-RLS result with PVE compensation (solid line), together with the true evoked brain activity (dashed line) used for quantitative comparison and evaluation of the performance of EMD-RLS algorithm.

Figure 7
Figure 7 shows the equivalent results of Figure 6 but for ∆[HHb].Comparing Figure 6a and 6b with Figure 7a and 7b, respectively, we can observe that the interference in the superficial layer that was introduced into the ∆[HHb] is relatively smaller than that in ∆[HbO 2 ].Thus, the trend of ∆[HHb] shown in Figure 7c and 7d was visible before EMD-RLS processing.The recovered results for ∆[HHb] and its block average results are shown in Figure 7e and 7f.Although the interference in ∆[HHb] is very small, the calculations shown in Figure 7e and 7f demonstrated that processed results are still somewhat better than the unprocessed results.The quantitative comparison [Figure 7g and 7h] was conducted with the PVE compensation and approximately 91% of the ∆[HHb]

Fig. 6
Fig.6 RLS adaptive filtering to remove physiological interference.(a) The time series of HHb calculated from Swith 45mm source-detector separation and (b) its block averaged result.(c) and (d) Reference measurements from FIG. 7 EMD-RLS to remove global interference.(a) The time series of concentration changes of HHb calculated from S-D1 with 10 mm source-detector separation, and (b) its block averaged result.(c) and (d) The time series of concentration changes of HHb calculated from S-D2 with 40 mm source-detector separation.(e) and (f) EMD-RLS result for the target measurement.(g) and (h) EMD-RLS result with PVE compensation (solid line), together with the true evoked brain activity (dashed line) used for quantitative comparison and evaluation of the performance of EMD-RLS algorithm.

TABLE 1
Haemodynamic parameters, thickness, and optical properties for each layer of an adult head model.

TABLE 2
Simulation parameters for amplitude and frequency of interference oscillations and haemodynamic changes.
In our simula- have used Monte Carlo simulations of a five-layered model of the human adult head to assess the effectiveness of EMD-RLS in global interference removal from fNIRS brain activity data.Our results have shown that the method can reduce global interference within measurements of ∆[HbO 2 ] and ∆[HHb].It should be noted that the EMD-RLS method is more effective for global interference removal in ∆[HbO 2 ] than in ∆[HHb] in our study.The reason could be interpreted that the interference magnitude of ∆[HbO 2 ] is very large and the improvement of denoising for ∆[HbO 2 ] is obvious.