EEMD and bidimensional RLS to suppress physiological interference for heterogeneous distribution in fNIRS study

Yan Zhang*||, Dan Liu*, Qisong Wang*, Xin Liu, Chunling Yang*¶||, Jinwei Sun*, Jingyang Lu and Peter Rolfe* *School of Electrical Engineering and Automaton Harbin Institute of Technology, Harbin, P. R. China School of Transportation Science and Engineering Harbin Institute of Technology, Harbin, P. R. China Intelligent Fusion Technology, Inc., Germantown, MD, USA zyhit@hit.edu.cn ¶yangcl1@hit.edu.cn


Introduction
Continuous-wave Near-Infrared spectroscopy (NIRS) which utilizes a constant-frequency or low-frequency modulated diode as the light source has been widely used in obtaining hemodynamic information of oxyhemoglobin, deoxyhemoglobin, etc. [1][2][3] The photodiode is often used as the detector.Compared to other cerebral function testing technologies such as electroencephalograph (EEG), magnetoencephalography (MEG), positron emission tomography (PET) and functional magnetic resonance imaging (fMRI), continuous-wave NIRS has shown great capability in various cerebral activity signal detection with advantages of simple construction, lowcost and easily realized miniaturization. 4However, the detection accuracy of NIRS cerebral function can be signi¯cantly in°uenced by the physiological activities. 5n recent years, many international research groups have devoted a great deal of e®ort to studying physiological interference issues. 6,7Physiological disturbances include cardiac cycle, respiration, spontaneous low frequency oscillation, and ultra-low frequency oscillation. 8The sources of its e®ects are two-fold.One is the physiological disturbance in the outer cerebral tissue such as the scalp, skull and cerebrospinal °uid.The other is the physiological disturbance in deep cerebral tissue such as cerebral gray matter and white matter. 9he physiological interferences are therefore also referred to as global disturbances or systemic physiological disturbances.
To deal with this type of interference e®ectively, the following four solutions have drawn wide attention in the international community: First, based on adaptive ¯ltering approach, the noise effect is suppressed with the help of auxiliary test equipment.For example, the signal that characterizes heartbeat and respiration is collected through the instruments such as blood pressure monitor, respirometer, 10 spirometer, carbon dioxide analyzer and pulse oximeter. 11The received signal can be used as reference signal to eliminate the disturbance using least mean square algorithm and Kalman ¯lter algorithm.Second, using a prior frequency (cardiac cycle, $1 Hz; respiration, $ 0:25 Hz; low frequency oscillation, $0:1 Hz), quasi-sinusoidal signal is regarded as the noise component, and the interference is tracked and then removed through the linear adaptive model and Kalman ¯lter algorithm. 12Third, considering the cerebral function activity has regional characteristics and the physiological interference has global characteristics, subtraction between the measured signals from the active and inactive regions can separate the cerebral function signals from the interference signals.The typical theory includes principal component analysis (PCA) 13 and independent component analysis (ICA). 14Fourth, the multi-distance measurement method consisting of pairs of light source-detectors is used: the close-separation source-detector pair is used as the reference signal, the far source-detector pair is used to acquire the synthetic signal containing a target signal and the interference, an adaptive ¯lter model 15,16 is constructed for the inhibition of interference.The above methods inhibit the physiological interference and obtain a good noise suppression e®ect to a certain extent but pay insu±cient attention to the heterogeneity of cerebral tissue.
In 2008, Saager et al. studied the heterogeneity of cerebral tissue based on the multi-distance measurement methods. 17Research shows that 90% of subjects have the cerebral tissue without serious heterogeneity, while the other 10% of subjects show signi¯cant heterogeneity.At this point, there are some di®erences in the distribution of capillaries, arteries and veins in di®erent parts of the cerebral tissue, and the traditional strati¯ed cerebral tissue model is ine®ective.In the near-infrared cerebral function tests, physiological disturbances come from the di®erent physiological activities, which may have multiple components.When the nonuniformity of cerebral tissue is serious, di®erent physiological activities will lead to di®erent physiological disturbances in di®erent positions.For this situation, a more viable solution is to separately estimate the di®erent types of interferences.One way is to obtain a reference signal for each physiological disturbance via a blood pressure monitor, a respirometer, etc., and then to track di®erent physiological disturbances using Kalman ¯ltering, but this method requires additional equipments.Another method is to use multiple prior sin/cos signals as the reference signals of physiological interference to estimate the physiological interference by Kalman ¯ltering, but this needs the prior knowledge of the frequency information of the physiological interference of the subjects.However, it is often hard for implementation in practice due to individual di®erences.
Considering multi-component characteristics of physiological disturbance and the uneven distribution of interfering components, a new method based on multi-distance measurement is proposed to improve the cerebral function detection accuracy.In this study, the near-spacing detector is used to acquire the extra cerebral tissue synthetic signals that are mainly caused by the physiological interference.The far-spacing signals acquired from far-detector are su±ciently sensitive to brain activity.The near channel signal is treated as the reference signal and is decomposed into a series of simple intrinsic mode functions using ensemble empirical mode decomposition (EEMD).The intrinsic mode functions have good instantaneous frequency characteristics and are suitable for nonlinear nonstationary signal analysis.The temporal correlation of neighboring sampling values is studied, and the estimated points are used as the center for smoothing processing.The physiological interference is further estimated using the recursive least square (RLS) algorithm.Monte Carlo simulation results show that this algorithm can e®ectively inhibit physiological interference and improve the detection accuracy of cerebral activity signal.

Theory 2.1. EEMD
Empirical mode decomposition (EMD), as an adaptive time-frequency analysis method, is widely used in biomedical signal processing and is especially suitable for analyzing nonlinear and nonstationary signals. 18However, there are some shortcomings in the application of EMD.One of the main drawbacks is the existence of serious mode mixing phenomenon in EMD, which leads to receiving di®erent time-scale or similar time-scale information in the intrinsic mode function components. 19The mode mixing in EMD makes the physical meaning of the individual intrinsic mode function unclear.In order to overcome the mode mixing problems, Huang et al. 20 propose the EEMD by adding a white noise to the original signal.As the white noise distributes uniformly in the whole timefrequency space, the problem of mode mixing can be avoided by adding the white noise signal to make the frequency distribution of the synthetic signal uniform.
Denote the measurement of oxyhemoglobin, deoxyhemoglobin concentration change xðtÞ, white noise series are added with the zero mean and constant variance.The signal is then decomposed into IMFs using EMD.Repeat the same procedure with di®erent white noise series each time and adopt the means of the ensemble corresponding IMFs to avoid the white noise e®ect to achieve the ¯nal IMF components.
where c j ðtÞ is the jth IMF component achieved using EEMD.c i;j ðtÞ is the jth IMF component by adding white noise series at the ith time, 1 i M.
The ¯nal results using EEMD are as follows: where rðtÞ is the ¯nal residue component describing the signal average trend.N is a positive integer describing the total number of the IMF components using EEMD.

EEMD-BRLS
The signal xðkÞ is the Á½HbO 2 ] or Á[HHb] acquired from the short-distance measurement S-D1 at instant k and is used to be decomposed for multiple components.The decomposed IMF components are achieved by EEMD.The ¯rst dimension is the order i of the EEMD decomposition, which corresponds to time-scales.The larger the value of i is, the smaller the time-scale of the IMF component is.The other dimension is the time series k, which indicates the signal changes with time.In the previous study, the EMD-RLS algorithm 21 considers only the relationship between IMF components of di®erent order at the current moment k and physiological interference and does not take into account that the signals of the adjacent sampling points are correlated.
In order to capture the correlation information in physiological disturbance estimation and achieve better estimation results, we explore the method of EEMD and bidimensional RLS (EEMD-BRLS) by investigating and analyzing the tapped delay line ¯lter.The block diagram of EEMD-BRLS is shown in Fig. 1.For a speci¯c time k, we select the window of length 2M þ 1 with the time k as the center.The two-dimensional weight w i;j ðkÞ as described in Eq. ( 3) describes the relationship between IMF components and physiological disturbances.This process is similar to smoothing operations in signal processing and is a posterior form of estimation since observed data after the point of interest are used for estimation.As more associated data accumulates around the estimated point, smoothing operation tends to give better estimation than ¯ltering.Here, the estimate of physiological disturbance can be expressed as where M is the half-window size, N is the decomposition order for EEMD.Equation (3) takes into account the correlation between the IMF order and time, so it is called EEMD-BRLS.The bidirectional weight can be interpreted as the i th IMF component c i ðkÞ passing through an FIR ¯lter of length 2M þ 1 with coe±cient w i;j ðkÞ.Therefore, it is equivalent to a ¯lter bank consisting of N ¯lters, each of which processes one IMF component and the ¯nal result is the sum of the outputs of all the ¯lters.Compared with EEMD-RLS, EEMD-BRLS considers the correlation between adjacent sampling points at the expense of increased computational complexity.Similar to EEMD-RLS, the cost function is the weighted square sum as follows: The parameter 0 < 1 is the forgetting factor that controls the memory span of the algorithm.Di®erentiate the equation with respect to w i;j ðkÞ, set it to zero, we can get, w l;m ðkÞR l;i; m;j ðkÞ ¼ p i;j ðkÞ; where Rðl; i; m; jÞ and p i;j ðkÞ are using Eq. ( 8), if the RðkÞ is not singular, wðkÞ can be calculated as The dimension of RðkÞ is ð2M þ 1ÞN Â ð2M þ 1ÞN.Given wðkÞ, the cerebral function signal can be extracted based on the matrix inversion lemma.
In the EEMD-BRLS algorithm, a larger M will improve the performance of the estimation because more time-related signals are considered.But once M reaches a certain value, the system performance cannot be improved signi¯cantly, even though the amount of computing has increased dramatically.Therefore, in practical application, the half-window length M is selected according to speci¯c system requirement.

Monte Carlo Simulation
Our study used multi-distance probe arrangement and two wavelengths of the light source (750 nm and 830 nm).A ¯ve-layer slab model consisting of scalp (sc), skull (sk), cerebrospinal °uid (csf), gray matter (gm) and white matter (wm) was used as a human adult head model (Fig. 2).We adopt particular source-detector separations to distinguish between processes occurring at di®erent tissue depths, as indicated by the `banana-shaped' regions that encompass statistically determined photon paths.The short separation with source-detector1 is used to probe the super¯cial tissue layer, and the long distance with source-detector2 to probe the deeper tissue layer.The short source-detector spacing was set to 10 mm and the long source-detector spacing was set to 40 mm.The generation of the NIRS time series was conducted with Monte Carlo simulations and standard analysis of functional hemodynamic changes.
The simulated cerebral function activity process is based on the block design experiment paradigm. 5he experiment is divided into the Stimu Phase and the Rest Phase (Fig. 3).The Stimu and Rest phases both lasted 20 s and alternated.The entire test time is 200 s.As shown in Fig. 3, sampling rate is 10 Hz, the entire measurement process collects 2000 points of data.For each layer of cerebral tissue, both cardiac cycle and respiratory-induced physiological disturbances are simulated, and the cerebral Stimu Phase is achieved through an increase in HbO 2 concentration and a decrease in HHb concentration.
In the process of the simulation, the scattering coe±cient of cerebral tissue in each layer is kept constant, and the change of light intensity is determined by the change of absorption coe±cient.Therefore, the change of HbO 2 and HHb concentration (Á[HbO 2 ] and Á[HHb]) can be calculated with the modi¯ed Lambert-Beer law as raw measurements of hemodynamic parameters over time during cerebral functional activity.In order to eliminate the physiological interference, we can

Results
Since the consideration of both di®erent orders of IMF component and temporal correlation guarantees the EEMD-BRLS algorithm with high accuracy, we ¯rst analyzed the e®ects of the algorithm in testing cerebral function.Then, targeting the issue of di®erent in°uences from various physiological activities on reference signal and synthetic signal caused by heterogeneity in the brain, the performances of RLS, EEMD-RLS and EEMD-BRLS algorithms are compared and analyzed using the mean squared errors (MSEs) from statistical results.
Simulation experiment is based on the ¯ve-layer cerebral tissue model suggested in earlier studies, and the Monte Carlo simulation program.We utilize the block design experimental paradigm of cerebral function, and take physiological interferences caused by cardiac cycle (frequency of 1.1 Hz), breathing (frequency of 0.18 Hz), low frequency oscillation (frequency of 0.1 Hz), as well as ultra-low frequency oscillation (frequency of 0.04 Hz) into consideration.Physiological interferences are appended to changes in hemodynamics of cerebral tissues, then, white noise with mean value of 0 and 1/100 standard deviation is generated to simulate random noise from actual measurements.Finally, the information of optical density changes is gathered through Monte Carlo simulation.The value of the length of half-window is set as 5 considering the complexity of the algorithm in order to get better simulation results.
Concentration changes of oxyhemoglobin reconstructed with EEMD-BRLS algorithm is shown in Fig. 4. Results are represented in 200 s time series on the left, and the block averages are on the right.Figures 4(a

Figs. 4(g) and 4(h), as well as Figs. 5(g) and 5(h).
From these results, we knew that utilizing the EEMD-BRLS algorithm based on multi-distance measurement could suppress physiological interferences and e®ectively extract cerebral function signals.
Under normal conditions, heterogeneity of hemodynamic parameters in extracerebral tissues and deep cerebral cortex is not severe, thus adaptive ¯ltering based on multi-distance measurement can reduce physiological interference very e®ectively.Saager et al. use both linear and circular light source-detector layouts to study the heterogeneity of local cerebral tissues, and observed obvious heterogeneity in the cerebral tissues in 10% of testers. 17onsidering the fact that degree of disturbance on di®erent cerebral regions is distinct from diverse physiological activities because of the heterogeneity in cerebral tissues, percentage of various physiological activities in reference signal and synthetic signal is changed during simulation.Since the e®ects from di®erent interferences are similar, we will only study the example of heartbeat interference.Although heartbeat is global physiological interference, heterogeneity of cerebral tissues resulted from di®erences in the volume, distribution and shape of blood vessels could lead to nonuniform in°uences by heartbeat on di®erent regions.Focusing on control parameters in previous studies, we assume that disturbance by heartbeat in the region detected by near-separation optical detector is sig-ni¯cantly larger than that in any other region, and the exceeding amount is measured in percentage.Zero percent represents the same amount of heartbeat interference compared to other tissues, which means homogeneity, while 100% means the interference from heartbeat in this region is one-fold higher than other tissues, and the heterogeneity is severe.
During Monte Carlo simulation, control parameters of physiological disturbance amplitudes among di®erent tissues (amplitude control parameters) are valued with 20% intervals, and the MSE is obtained by taking the average of 100 independent experiments.The simulation results are shown in Fig. 6, where Fig. 6(a) contains the MSE of concentration changes of oxyhemoglobin, and Fig. 6(b) contains the MSE of concentration changes of deoxyhemoglobin.From Fig. 6(a), we conclude from the MSE value that the RLS algorithm could most accurately perform cerebral function detection while the amplitude control parameter remain unchanged (cerebral tissues have ideal homogeneity).When the amplitude control parameter increases (tissues begin to show severe heterogeneity), the e®ectiveness of RLS algorithm deteriorates; however, it does not a®ect EEMD-RLS or EEMD-BRLS algorithm signi¯cantly, showing the great robustness of the latter two algorithms.Decomposing reference signals, assigning di®erent values to various IMF components, and adjusting weights self-adaptively, where unrelated IMF components are assigned with small weights and large weights are assigned to strongly related components, these characteristics of the EEMD earned better interference estimations for the algorithm.
During optimization, EEMD-BRLS took the relevance of adjacent time of each IMF component into consideration, and used smoothing to process algorithms, thus the MSE calculated with EEMD-BRLS is signi¯cantly smaller than that with EEMD-RLS algorithm.Considering Fig. 6(b), the trend of RLS algorithm is similar to that in Fig. 6(a), and EEMD-BRLS algorithm also performed better than the EEMD-RLS algorithm.Compared to the results of changes in concentration of oxyhemoglobin, RLS algorithm performed very well compared to the other two EEMD methods when the change in amplitude control parameters is less than around 45% in the deoxyhemoglobin result.There are two reasons for this phenomenon, the ¯rst is that physiological interference has smaller amplitudes for the detection of oxyhemoglobin, which means heterogeneity has far less e®ect on the detection of change in concentration of oxyhemoglobin than that on deoxyhemoglobin; secondly, although EEMD decomposition can provide intrinsic mode functions of di®erent simple components, °uctuations may occur during decomposition process, and could bring errors to the ultimate optimization results.Thus, when the tissue heterogeneity is not severe, advantages of EEMD-RLS and EEMD-BRLS algorithms could not o®set the °uctuations during EEMD decomposition and could lead to results worse than the RLS algorithm.
To analyze the performance of EEMD-BRLS algorithm with di®erent half-window lengths, statistical analyses are performed on MSEs when halfwindow length (M) ranged from 1 to 5, and the result is shown in Fig. 7. Corresponding MSE values with alternating concentrations of oxygenated and deoxyhemoglobin are plotted in Figs.7(a) and 7(b), respectively.Similar conclusions are drawn from both ¯gures, the larger the M is, the smaller the MSE value is, which means the better results are obtained when half-window length is large.However, large M values do not necessarily mean good results.Our simulation results showed that when M is larger than 5, the trend for smaller MSE value slows down, and when M is larger than 8, MSE does not signi¯cantly decrease.Overall, considering the complexity of algorithms, when studying the extraction of cerebral function signals, half-window length of 5 is the most appropriate, and MSE value increases along with the increase in amplitude control parameters; however, EEMD-BRLS algorithm has better robustness compared to RLS algorithm.

Discussion
Physiological interference can greatly degrade the performance of fNIRS measurement of evoked brain activity response.Many di®erent methodologies have been explored to resolve this problem.One of the approaches is to identify and separate the interference components in fNIRS study including using auxiliary physiological measurements and to analyze the signal with an adaptive ¯ltering algorithm.The instruments such as electrocardiogram (ECG), the pulse oximeter, chest band respirometer, spirometer and capnograph can be used to achieve auxiliary physiological measurements.This method has the drawback of the indispensability of additional equipment.Another approach adopted adaptive ¯ltering methods to separate the brain activity response and the physiological interference using the priori frequencies.The priori frequencies (e.g., the cardiac, $ 1 Hz; respiratory, $ 0:25 Hz; and Mayer wave frequency, $ 0:1 Hz) have also been assumed to be the sine/cosine terms. 12However, this methodology is not appropriate due to individual di®erence.Thus, recently multi-distance measurement and adaptive ¯ltering have been proposed to suppress the physiological interference and the advantage of this method is convenient and e®ective to extract the signal for real-time application.
It should be noted that in the measurement method, the performance of adaptive ¯ltering declines due to heterogeneous problem.In this current paper, we employ the time-frequency analyzing method to decompose the reference signal and then using RLS algorithm to estimate the brain activity response.The accuracy of the EEMD-RLS algorithm and EEMD-BRLS algorithm depends on the reference signal decomposition.Therefore, this paper uses the improved EMD algorithm EEMD to improve the accuracy of the decomposition.In addition, the EEMD-RLS and EEMD-BRLS algorithms have no obvious di®erences when the tissue heterogeneity is not serious, while the RLS algorithm is easy to implement and can be utilized for real-time signal processing without EEMD processing procedure.The EEMD-RLS and EEMD-BRLS algorithms can signi¯cantly improve system robustness when the distributions of disruptions caused by di®erent physiological activities in different regions are not consistent.
An important matter that should be considered here is the EEMD-RLS and EEMD-BRLS algorithms can not only conduct the cerebral function signals detection, but also provide the extracerebral hemodynamic changes for cerebral function study as they can provide the time-frequency information of the outer layer of the brain via EEMD.During the cerebral stimulation, the hemodynamic parameters of extracerebral layers are not always independent of the response to the cerebral cortex.For example, the cardiac cycle increase during ¯nger motility experiments increases hemodynamic correlations between the outer layer of the cerebral and the cerebral cortex, which may lead to over-¯tting in cerebral function testing via multiple-range measurements.The cerebral function signals can be extracted by EEMD-RLS and EEMD-BRLS approaches.At the same time, the information of frequency components regarding the outer layers of cerebral tissue can be achieved via EEMD method to further analyze the hemodynamics of the outer layer of cerebral tissue.
A truly rigorous evaluation of this method requires the simulation of the heterogeneous head model, which unfortunately, is not easy to realize for the in vivo experiments as such experiments need the selection of the heterogeneous subjects.Meanwhile, the PVE e®ect in vivo cannot be exactly compensated and the quantitative comparison of the recovery response and the true response of brain activity is di±cult.Thus, the control parameters are altered with ¯xed steps to imitate the heterogeneous head tissue.Monte Carlo simulation is further implemented for quantitative analysis here as the preliminary study.As for the deep study of this method in vivo with practical NIRS measurement, this will be the subject of our on-going research and will be published in due course.

Conclusions
In this paper, multi-distance measurement method has been investigated regarding eliminating the

Fig. 2 .
Fig. 2. Schematic of a multi-distance fNIRS probe arrangement together with the ¯ve layer model of the head.
) and 4(b) are acquired with nearseparation detector, while Figs.4(c) and 4(d) are acquired with far-separation detector.Based on the multi-distance measurement and EEMD-BRLS algorithm, changes in oxyhemoglobin concentration as well as the coherent averages are represented in Figs.4(e) and 4(f).After partial volume e®ect (PVE) compensation of changes in oxyhemoglobin concentration collected with the EEMD-BRLS algorithm, results and block averages are shown in Figs.4(g) and 4(h), in which the full line represents extracted value, and the dashed line represents the actual value.Concentration changes of deoxyhemoglobin reconstructed with EEMD-BRLS algorithm are shown in Fig. 5. Results are represented in 200 s time series on the left, and coherent averages are on the right.Figures 5(a) and 5(b) are acquired with near-end detector, while Figs.5(c) and 5(d) are acquired with remote detector.Based on multidistance measurement and EEMD-BRLS algorithm, changes in deoxyhemoglobin concentration as well as the coherent averages are represented in Figs.5(e) and 5(f).After PVE compensation of changes in deoxyhemoglobin concentration are collected with the EEMD-BRLS algorithm, results and coherent averages are shown in Figs.5(g) and 5(h), in which the full line represents extracted value, and the dashed line represents the actual value.Observing Figs. 4 and 5, we noticed that there is obvious physiological interference of synthetic signals in Figs.4(a) and 4(c), as well as in Figs.5(a) and 5(c), along with their block averages.Extracted changes in hemodynamic parameters and actual changes in hemodynamics are quite close from

Fig. 6 .Fig. 7 .
Fig. 6.MSE versus the variation of amplitude control parameters for three di®erent algorithms.(a) The MSE of concentration changes of oxyhemoglobin and (b) the MSE of concentration changes of deoxyhemoglobin.