International Journal of Scientific & Engineering Research, Volume 4, Issue 7, July-2013 2280

ISSN 2229-5518

Motion Artifact Reduction in fNIRS Signals by

ARMA Modeling Based Kalman Filtering

Mehdi Amian, S. Kamaledin Setarehdan

Abstract— Functional Near Infrared Spectroscopy (fNIRS) is a technique that is used for noninvasive measurement of the oxyhemoglobin (HbO2) and deoxyhemoglobin (Hb) concentrations in the brain tissue. Since the ratio of the concentration of these two agents is correlated with the neuronal activity, fNIRS can be used for the monitoring and quantifying the cortical activity. The portability of fNIRS makes it a

good candidate for studies involving subject’s movement. The fNIRS measurements, however, are sensitive to artifacts generated by

subject’s head motion. This makes fNIRS signals less effective in such applications. In this paper, the autoregressive moving average (ARMA) modeling of the fNIRS signal is used for state-space representation of the signal which is then fed to the Kalman filter for estimating the motionless signal. Results were compared to the autoregressive model (AR) based approach and showed that the ARMA models outperform AR models because of better modeling of the fNIRS signals. We showed that the signal to noise ratio (SNR) is about 3 dB higher for ARMA based method.

Index Terms— Brain, Gaussian noise, linear model, state estimation.

—————————— ——————————

1 INTRODUCTION

HIS paper’s focus is on Functional near infrared spectros- copy (fNIRS) as a relatively recent technique for noninva- sive measurement of oxygenated hemoglobin (oxy-Hb)
and deoxygenated hemoglobin (deoxy-Hb) concentrations in the human brain [1]. Other applications include quality con- trol, pharmacology and medical diagnoses [2]. Basically, a typ- ical fNIRS system is composed of one or a number of light sources in near infrared (NIR) range (700-900 nm), and several detectors that collect the reflected photons from the brain tis- sue. In this spectrum, the blood and vital tissues are relatively transparent and by analysing the collected light intensities (fNIRS signals), the properties of the medium through which the light has passed can be extracted [2],[3].
This technique is affordable, portable, and capable of being used in real field applications such as monitoring pilots dur- ing flight [2],[4],[5]. Also, fNIRS is safe compared to other im- aging techniques such as X ray imaging, positron emission tomography (PET), nuclear medicine and computed tomogra- phy (CT). As such, it has been widely used for studies with vulnerable populations such as neonates. However, such ap- plications entail inevitable head movements. As the head mo- tion can increase the blood flow through the scalp and rarely causes an increase in brain’s blood pressure [2] , therefore, motion artifact may change original signal and lead to incor- rect result or misguided diagnose.
Reducing head motion artifacts is therefore a key aspect in signal processing area of fNIRS studies and there have been a

————————————————

Mehdi Amian has received an M.Sc. degree in 2011 in biomedical eni- gineering from University of Tehran, Tehran, Iran. E-mail: Mehdi.amian@alumni.ut.ac.ir

S. Kamaledin Setarehdan is associate professor of biomedical engineering at

University of Tehran. E-mail: ksetareh@ut.ac.ir
number of attempts. Adaptive filtering is used as the main artifact removal method [6]. The Wiener filter was also applied to the FNIRS signals [7]. A wavelet based approach was used by Molavi et al. [8]. Although these and other methods could reduce the motion artifact from fNIRS signals but each method has some specific requirements and limitations.
For example, in the adaptive filter based method [6], the algo- rithm needs additional hardware and sensors which makes it more complex and costly. The Wiener filter based method does not need extra instruments [7]; however, it needs to have the whole data simultaneously and therefore is not applicable in real time applications.
In contrast, the Kalman filter based approaches, such as the one presented by Izzetoglu et al. do not need extra sensors and also can perform motion artifact reduction in real time [1]. In the Kalman filter based methods it is necessary to model the input signal by a linear model first. For example, an auto- regressive (AR) model is used in [1] to model the fNIRS sys- tem first, and then the AR model is transformed into state space representation. Next the Kalman filter is applied to es- timate motionless data from motion corrupted data.
In this paper, we have employed the Autoregressive Moving Average (ARMA) model instead of the commonly used AR model for motion artifact removing from the NIRS signals for the first time. Our results showed improved motion reduction over previous studies. The ARMA model is more comprehen- sive than AR. Selecting an appropriate transform, as will be presented later, also leads to better results. Our results show an improvement of 3 dB in the signal to noise ration of the fNIRS signals.
The rest of the paper is organized as follows. In Section 4 the proposed algorithm is described in details. The data set used in this work is explained in that section as well. The results of the application of the proposed method to this data set are brought in section 5. Finally Section 6 concludes the paper
.

IJSER © 2013 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 4, Issue 7, July-2013 2281

ISSN 2229-5518

2 MATERIALS AND METHODS

2.1 fNIRS Data

Detailed The fNIRS data that used in this work were recorded by a three channel fNIRS probe attached to the subject’s fore- head at CONQUER CollabOrative at Drexel University (3508
Market Street, Philadelphia, PA 1904, Drexel University). fNIRS data were collected using a continuous wave fNIRS system. The fNIRS system is composed of three subsystems: 1) fNIRS sensors that consist of one light source and three photo detectors. The light source is a multi-wavelength light emit- ting diode (LED) manufactured by Epitex Inc. type L4*730/4*850 – 40Q96-I. The LED comes in a STEM TO- 5 package at 730 nm and 850 nm wavelengths with an output power of 5 to 15 mW. The photo detectors are manufactured by Burr-Brown Corporation type OPT101 and come in an 8- pin DIP package. 2) A control box for operating the LEDs and photo detectors. 3) A desktop computer running the COBI Studio software developed in the laboratory for data acquisi- tion and real-time data visualization. Three channels are used to record fNIRS signals. Source-detector distances for the channels are 2.8, 2.8, and 1 cm. Sampling frequency is 2 Hz. Six healthy, right handed individuals (3 males) with no history of neurological, psychological, or psychiatric disorders who were analgesic-free were recruited from the Drexel University community. All participants signed the informed consent form approved by the Institutional Review Board (IRB) at Drexel University. Examples of the recorded signals are shown in Fig.
1.

Fig. 1. Example of fNIR signals including three stages, stage 1 is rest, and then subject was moving his/her head during stage 2, and stage 3 is post motion, solid line was used for deoxy-Hb and the dotted line for oxy-Hb.

2.2 Protocol

For In this study we have concentrated on the motion artifact reduction in the NIRS signals. This artifact is generated due to the head motions of the subject. Fig. 1 shows the motion cor- rupted oxy-Hb and deoxy-Hb in solid and dotted lines respec- tively. The following protocol was used.
The first stage includes one minute baseline signal in which the subject was instructed to sit still and relax (rest stage). This stage is considered as motionless signal. Then the subject was instructed to move his/her head in a steady frequency around
0.3 Hz for 30 seconds (three movements in each ten seconds). This motion artifact is considered as slow head motion. The experiment ended with two minutes post recording with no motion (post motion stage).
As a result each of the two oxy-Hb and deoxy-Hb signals was contaminated by the motion artifacts.
Fig. 2 shows a typical fNIRS signal including two rest and head motion stages. The post motion stage is removed.

Fig. 2. Example of a typical fNIR signal with rest and head mo- tion stages (stage 1 and stage 2 respectively).

In our method we first model the fNIRS signal using the commonly used AR model as well as ARMA model. The linear model is then transformed into state space representations. Finally Kalman filter is applied using the information in the state space matrices and, the motionless signal is estimated. The results are compared by calculating ∆SNR for both AR and ARMA based methods.

2.3 Linear Model

In most cases, it is likely that not the entire recorded data is contaminated by motion artifact. Usually, motion artifact cor- rupted parts of the signal are easily distinguishable from mo- tionless parts. We will use motionless parts of NIRS data as well as motion corrupted parts to calculate required parame- ters.
At first, the fNIRS motionless signal (stage 1) is modeled as an autoregressive (AR) model. The optimal AR order for the sys- tem (fNIRS motionless signal) was obtained by Akaike Infor- mation Criterion as N=4. We can see form (1) that in an AR model of order of 4, how a sample of signal in the present time
point 𝑘 is linearly relating to its previous values and to the noise. In (1), 𝑤𝑘 is a white Gaussian noise that its variance will
be calculated later.

IJSER © 2013 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 4, Issue 7, July-2013 2282

ISSN 2229-5518

𝑥𝑘 = 𝑎1 𝑥𝑘−1 + ⋯ + 𝑎4 𝑥𝑘−4 + 𝑤𝑘 (1)
where 𝑎𝑖 , i=1,...,4 are AR coefficients. We did not obtain the
coefficients manually, but we used MATLAB to obtain them
directly and automatically. The model which software as-
sumes for an AR model is as following.
M(q) y(n) = e(n) (2) where, 𝑦(𝑛) is a discrete function of 𝑛, e(n) is a Gaussian white noise, and 𝑀(𝑞) includes AR coefficients and is as fol-
lows;
𝑀(𝑞) = 1 + 𝑚1𝑞−1 + 𝑚2 𝑞−2 + 𝑚3 𝑞−3 + 𝑚4𝑞−4 (3)
artifact we assume:
𝒗�𝒌+𝟏 = 𝒛𝒌+𝟏 − 𝐻𝒙𝒌 (11)
where, 𝐻 is the same as that it was brought in (8).
In (11) index 𝑘 includes those time points during which
head motions of the subject exist.
We assume the statistics of the system and the measure-
ment noise to be wk and vk respectively, and also they are in-
dependent of each other, with white Gaussian distributions
𝒘𝒌~ 𝑁 (0, 𝑄), 𝒗𝒌 ~ 𝑁(0, 𝑅). Kalman filtering equations are as
followings:
1. Time update equations:
𝒙�𝒌 = 𝐴𝒙�𝒌−𝟏
In (3) operator 𝑞 is delay tap in fact, that after operating on
𝑃𝑘 = 𝐴𝑃𝑘−1 𝐴
+ 𝑄
the function 𝑦(𝑛), creates function lags.
2. Measurement update equations:
MATLAB uses the forward-backward algorithm to estimate
𝐾𝑘 = 𝑃𝑘 𝐻 (𝐻𝑃𝑘 𝐻
+ 𝑅)
the parameters of the linear model. In the forward model sum
of a least-square criterion is minimized and analogous criteri-
on for a time-reversed model.
At the next stage, the linear model is transformed into the
In which 𝑃
𝒙�𝒌 = 𝒙�𝒌 + 𝐾𝑘 (𝒛𝒌 − 𝐻𝒙�𝒌 )
𝑃𝑘 = (𝐼 − 𝐾𝑘 𝐻)𝑃𝑘
is a priori covariance matrix (error covariance
state space representation as below;
matrix before update) as follows:
𝒆𝒌 − = 𝒙𝒌 − 𝒙�𝒌
, 𝑃𝑘 = 𝐸[𝒆𝒌
𝒆𝒌 ]
Where,
where, 𝐾𝑘 is Kalman gain matrix [9].
Then output signal of the Kalman filter and the contami-

xk

wk

nated input signal are depicted in the plots below. One can see

   

strength and success of the method in eliminating motion arti-

xk −1

,  0 
(6)
facts by visually comparison of the two signals. A quantitative

  

 

xk N +1 N ×1

  

 

 0  N ×1

criterion as ∆𝑆𝑁𝑅 is as below.
∆𝑆𝑁𝑅 = 𝑆𝑁𝑅𝑒 − 𝑆𝑁𝑅𝑖 (19)
where, the estimation of the signal to noise ratio, SNRe, is
computed as
𝒛𝒌 is the motion corrupted signal vector, 𝒙𝒌 is the motion- less signal vector, 𝒘𝑘 is the measurement noise vector and 𝒗𝒌
is the motion artifact vector. In general, 𝐴 is an N×N matrix and 𝐻 is an1×N row vector as followings:

SNRe

σ 2

= ( x )

e

(20)

2 2

− mN

mN −1

mN −2

 − m2

m1

where, 𝜎𝑥 is variance of the motionless signal, and 𝜎𝑒 is the
variance of estimation error that is difference between the NIR
signal and the modified signal after applying filtering method.

 1 0 0

 0 1 0

 0

 0 0

 0 0

 1 0

N ×N

(7)
𝑒(𝑛) = 𝑥(𝑛) − 𝑥�(𝑛) (21) The input signal to noise ratio, SNRi, is computed as;

2

x

and
𝐻 = [1 0 … 0 0]1×𝑁 (8)

SNRe = ( 2 )

e

2

(22)
As it was mentioned before, N=4 here.
To estimate the covariance of system's noise, we assume
𝒘� 𝒌+𝟏 = 𝒛𝒌+𝟏 − 𝐺𝒙𝒌 (9)
where
𝐺 = [−𝑚1 − 𝑚2 − 𝑚3 − 𝑚4 ] (10)
𝒙𝒌 is the AR sample vector and 𝒛𝒌 is the measured signal
vector. Index k contains those time points in which we have
where, 𝜎𝑣 is the variance of the motion artifact.
Now we develop an ARMA model for the NIR data that is
a richer model than AR. Optimal order of model is obtained as
N=4. The parameters of model directly and automatically
were estimated using MATLAB. The model that MATLAB
assumes for ARMA is as follows:
𝑁(𝑞)𝑦(𝑛) = 𝐶(𝑞)𝑒(𝑛) (23)
where
motionless signal. Eventually, variance of w� k+1 is calculated
and considered as an estimation of the system's noise. Similar-
and
𝑁(𝑞) = 1 + 𝑛1 𝑞−1
+ 𝑛2 𝑞−2
+ 𝑛3 𝑞−3
+ 𝑛4 𝑞−4
(24)
ly, to estimate the variance of the measuring noise or motion
𝐶(𝑞) = 1 + 𝑐1 𝑞−1 + 𝑐2 𝑞−2 + 𝑐3 𝑞−3 + 𝑎4 𝑞−4 (25)

IJSER © 2013 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 4, Issue 7, July-2013 2283

ISSN 2229-5518

For estimating ni and ci the software uses a recursive ap-
proach through which minimizes estimation error. The cost
function is the determinant of input covariance matrix.
Then this liner model is converted into a space state repre-
sentation. There are numerous state space representations for
a linear system, however we chose the following representa-
tion that resulted in better outcome:
Where

TABLE 1

RESULTS OF APPLYING AR BASED METHOD TO SIMULATED DATA

Signal index

1

2

3

4

5

6

∆𝑆𝑁𝑅

10.39

4.12

9.08

7.70

10

8.94

In Fig. 3, (a) an example of a white and Gaussian noise is shown. Such noise added to an fNIRS signal (stage 1) and the resulted signal is shown in (b), and in (c) the cleaned signal by AR based method is printed in solid line and the original sig- nal from stage 1 in dashed- line.

2 2

2

0.6

0.4

0.2

0

-0.2

-0.4

-0.6

-0.8

0 10 20 30 40 50 60

0.5

0

-0.5

-1

-1.5

-2

0 10 20 30 40 50 60

At this step the variance of the system's noise, wk, should
be estimated. ARMA model can be written as:
(a)

-0.2

(b)

Xhat

time (sec)

𝑦(𝑘) = −𝑛1 𝑦(𝑘 − 1) − ⋯ − 𝑛4 𝑦(𝑘 − 4) + 𝑒(𝑘)
+ 𝑐1 𝑒(𝑘 − 1)+. . . +𝑐4 𝑒(𝑘 − 4)
We can rewrite (30) as
𝑒(𝑘) = 𝑦(𝑘) + 𝑛1 𝑦(𝑘 − 1) + ⋯ + 𝑛4 𝑦(𝑘 − 4) − 𝑐1 𝑒(𝑘 − 1)
− ⋯ − 𝑐4 𝑒(𝑘 − 4)
(30)
(31)

-0.4

-0.6

-0.8

-1

data

One can obtain value of noise iteratively from (31) by
knowing the initial value of the noise. We assume zero value
for initial values of noise. After finding noise in this way, the
variance of the noise can be easily calculated.
To estimate the variance of the measurement noise, vk, that
we assume it to be the motion artifact, the procedure is as fol-
lowing. An estimation of measurement noise can be in form of
(c)

-1.2

-1.4

-1.6

0 10 20 30 40 50 60

time(sec)

(32).
𝒗�𝒌+𝟏 = 𝒛𝒌+𝟏 − 𝐻𝒙𝒌 (32)

Fig. 3. (a). A typical pattern of an artificial noise (we assume as white and Gaussian) (b). The noise in (a) added to an fNIRS signal (c). The contami- nated signal in (b) is cleaned by AR based method (dashed line) and the original motionless signal (solid line).

After finding 𝒗�𝒌 ,variance of measurement noise can be
easily calculated.

3 RESULTS

3.1 Simulated Data

To evaluate the accuracy of the method, the proposed algo- rithms, both AR and ARMA based ones, are applied to the simulated data. Simulated data is in fact the same clear NIR signal that is contaminated by a known artificial noise. Some white noises with definite covariance added to one of the clear NIR signals, and then the variance of such noise is estimated through the algorithms explained before. If algorithm works well the variance will be estimated truly.
Results of applying AR based algorithm to simulated data shown in Table 1.
The result of testing algorithm relating to ARMA based meth- od is shown in Table2.

TABLE 2

RESULTS OF APPLYING ARMA BASED METHOD TO SIMULATED DATA

covariance of

Gaussian noise

0.10

0.07

0.08

0.09

0.07

0.12

0.10

estimated covariance of noise

0.12

0.08

0.10

0.10

0.08

0.14

0.14

∆𝑆𝑁𝑅

3.21

1.63

2.63

2.45

1.82

2.83

3.70

In Fig. 4 (a) we demonstrated another white and Gaussian noise. Such noise added to a motionless fNIRS signal (stage 1)

IJSER © 2013 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 4, Issue 7, July-2013 2284

ISSN 2229-5518

and the corrupted is shown in (b). In (c) the cleaned signal by the ARMA model and the original motionless signal are shown.

2.5

Xhat data

2

0.8

0.6

0.4

0.2

0

-0.2

-0.4

-0.6

-0.8

0 10 20 30 40 50 60

0.5

0

-0.5

-1

-1.5

-2

0 10 20 30 40 50 60 time (sec)

1.5

1

0.5

0

(a) (b)

-0.2

0 5 10 15 20 25 30 time (sec)

-0.4

-0.6

-0.8

Xhat data

Fig. 5. An actual experimental contaminated signal (solid line) from stage 2 where the subject was moving his/her head and denoised signal (dashed line) after using AR based method.

-1

-1.2

-1.4

-1.6

0 10 20 30 40 50 60 time (sec)

(c)
Notice that, unlike to Fig. s 3, 4, the motion corrupted part of signal (stage) 2 is platted in Fig. 5. However, in Fig. s 3, 4, the motionless part of signal (stage 1) was adulterated with artificial noise. Therefore, the range of horizontal axis in Fig- ure 5 (and also Fig. 6) is different from Fig. s 3, 4.
The quantitative results due to computing through ARMA

Fig. 4. (a). A typical pattern of an artificial noise (we assume as white and Gaussian) (b). The noise in (a) added to an fNIRS signal (c). The contami- nated signal in (b) is cleaned by ARMA based method (dashed line) and the original motionless signal (solid line).

based method brought in Table 4.

TABLE 4

RESULTS OF APPLYING ARMA BASED METHOD TO REAL DATA

3.2 Real (actual) Data

In this part we use the proposed method to estimate clear data from motion corrupted real NIR signals. Table 3 shows the
calculated values of ∆SNR for several signals in AR method.

TABLE 3

RESULTS OF APPLYING AR BASED METHOD TO REAL DATA

An example of the motion artifact reduction of a contami- nated real signal in ARMA method, is shown in Figure 6 where the contaminated signal drawn as solid line and the cleared on as dashed line.

Signal index

1

2

3

4

5

6

∆𝑆𝑁𝑅

10.39

4.12

9.08

7.70

10

8.94

Fig. 5 demonstrates one of the actual motion corrupted fNIRS signals (stage 2) and its clear estimated one by pro- posed AR based method.

2.5

2

1.5

Xhat data

1

0.5

0

0 5 10 15 20 25 30

time (sec)

Fig. 6. Real contaminated signal (solid line) and denoised signal

(dashed line) after using ARMA based method.

IJSER © 2013 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 4, Issue 7, July-2013 2285

ISSN 2229-5518

4 CONCLUSION

In the previous section, the result of applying AR based method as well as ARMA based method elucidated. The cleared signals of both simulated and real motion artifact cor- rupted, in ARMA based approach have higher quality in com- parison to AR. In essence, ARMA model includes more terms rather than AR. The additional terms in ARMA lead to a better
fitted linear model of NIRS system. Furthermore, ∆𝑆𝑁𝑅 was
introduced as a quantitative measure that establishes a numer-
ical criterion for comparing alternative methods. As we saw
from Tables 1 and 2, ∆𝑆𝑁𝑅 grows more than 2 dB in ARMA
results. From Tables 3 and 4, one can see an improvement of
about 3 dB in ∆𝑆𝑁𝑅 in ARMA results in comparison to AR.

ACKNOWLEDGMENT

We highly appreciate Professor Kambiz Pourrezaei and Dr. Zeinab Barati, PhD, from Drexel University, Philadelphia, PA, who provided us with the data used in this work. Our grateful thanks go to them for their useful and helpful guide to com- plete our work.

REFERENCES

[1] M. Izzetoglu, P. Chitrapu, S. Bunce and B. Onaral, “Motion artifact cancella- tion in NIR spectroscopy using discrete Kalman filtering,” BioMedialEngineer- ing Online, vol. 9, pp. 1-10, 2010.

[2] F. Mattehews, B.A. Pearlmutter, T.E. Wrad, C. Soraghan, C. Markham, “He- modynamics for Brain-Computer Interfaces,” IEEE Signal Processing Magazine, pp. 87-94, 2008.

[3] J. Gervain, J. Mehler , J.F. Werker, C.A. Nelson, G. Csibra, S. Lioyd-fox, M.

Shukla, R.N. Aslin, ”Near-infrared spectroscopy: A report from the McDon- nell infant methodology consortium,” Developmental Cognitive Neuroscience, vol. 1, pp. 22-46, 2011.

[4] M. Cope, “The application of near infrared spectroscopy to noninvasive

monitoring of cerebral oxygenation in the newborn infant,” Ph.D. dissertion, Department of Medical Physics Bioengineering, College London, London, England, 1991

[5] F.F. Jobsis, “Noninvasive infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science, vol. 198, pp. 1264-

1267, 1977.

[6] M. Izzetoglu, A. Devaraj, K. Izzetoglu, S. Bunce, B. Onaral, “Motion artifact removal in fNIR signals using adaptive filtering,” Proceeding of the EMBS,

2003.

[7] M. Izzetoglu, A. Devaraj, S. Bunce, B. Onaral, “Motion artifact cancellation in NIR spectroscopy using wiener filtering,” IEEE Transaction on Biomedical Engi- neering, vol. 52, 2005.

[8] B. Molavi, G. Dumont, B. shadgan, “Motion artifact removal from muscle NIR spectroscopy measurements,” IEEE Canadian Conferrence on Electrical and Computer Engineering , 2010.

[9] M. S. Grewal, A. P. Andrews, “Kalman filtering: Theory and Practice Using

MATLAB,” Second Edition, John Wiley & Sons Inc. , 2001.

IJSER © 2013 http://www.ijser.org