Multi layer spectral decomposition technique for ERD estimation in EEG μ rhythms: An EEG–fMRI study

Communicated by Prof Juan Manuel Gorriz Saez

Accepted Manuscript

Multi Layer Spectral Decomposition Technique for ERD Estimation in EEG µ Rhythms: An EEG-fMRI study Saideh Ferdowsi, Vahid Abolghasemi PII: DOI: Reference:

S0925-2312(17)31667-3 10.1016/j.neucom.2017.10.016 NEUCOM 18996

To appear in:

Neurocomputing

Received date: Revised date: Accepted date:

7 April 2017 29 July 2017 17 October 2017

Please cite this article as: Saideh Ferdowsi, Vahid Abolghasemi, Multi Layer Spectral Decomposition Technique for ERD Estimation in EEG µ Rhythms: An EEG-fMRI study, Neurocomputing (2017), doi: 10.1016/j.neucom.2017.10.016

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • In this paper ERD estimation in mu rhythm when EEG acquired simultaneously with fMRI is addressed. • Singular Spectrum Analysis is used for spectral factorization.

CR IP T

• The results of ERD estimation are compared with those of conventional methods that confirm superiority of the proposed method.

AC

CE

PT

ED

M

AN US

• EEG-fMRI combination is performed as complementary experiment.

1

ACCEPTED MANUSCRIPT

Multi Layer Spectral Decomposition Technique for ERD Estimation in EEG µ Rhythms: An EEG-fMRI study Saideh Ferdowsia , Vahid Abolghasemia School of Electrical Engineering and Robotics, Shahrood University of Technology, Shahrood, IRAN

CR IP T

a

Abstract

CE

PT

ED

M

AN US

In this paper, power suppression detection in mu band known as event-related desynchronization during a simultaneous EEG-fMRI recording is addressed. Simultaneous EEG-fMRI experiments provide great opportunities to recognize relation between the active areas in the brain, discovered by fMRI, and EEG events. However, decoding EEG is challenging as a result of artifacts which MRI scanner applies on EEG. Ballistocardiogram (BCG) is one of the main destructive artifacts which deteriorates EEG rhythms specially in mu band. The proposed method, is a supervised time-frequency decomposition algorithm which estimates underlying rhythms of EEG signals recorded in MRI scanner. In the proposed method, a multi layer decomposition approach using three criteria is developed to extract brain oscillations. This goal is achieved by eliminating the artifact related components at each decomposition layer as a result of removing useless extracted components in each layer. Performance of the proposed method is evaluated using synthetic and real EEG. The achieved results from synthetic data confirm the ability of the proposed method to extract mu rhythms when different levels of noise exist. In addition, proposed algorithm is applied on a set of real EEG acquired in a simultaneous EEG-fMRI recording. Estimated ERD confirms the superiority of the developed algorithm in terms of estimating changes of mu rhythms. Moreover, an EEG-fMRI integration is performed to explore the correlation between blood oxygenation level dependent (BOLD) in fMRI and ERD of mu rhythm in EEG.

AC

Keywords: EEG–fMRI combination, Event related desynchronization in mu band, Singular spectrum analysis, Spectral factorization

1

2 3 4 5 6 7 8 9

10

1. Introduction Event-Related Desynchronization (ERD) and event-related synchronization (ERS) are two main phenomena representing changes in the power of specific frequency band of ongoing EEG activity. Neural activity acquired by EEG oscillates in various frequency bands consist of delta (up to 4 Hz), theta (4 - 7 Hz), alpha or mu (8 - 12 Hz), beta (13 - 30 Hz) and gamma (> 30 Hz) [19]. ERD and ERS refer to decrease and increase in the power of EEG brain rhythms respectively. ERD is known as suppression of the power of EEG rhythms. In contrast, ERS refers to sharp increase of the power of EEG rhythms. Various tasks may cause ERD and ERS in different EEG rhythms and brain areas. For example, movement preparation or motor Preprint submitted to Neurocomputing

October 30, 2017

ACCEPTED MANUSCRIPT

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

AC

44

CR IP T

16

AN US

15

M

14

ED

13

imagery behaviour causes ERD and ERS in both mu and beta bands[13, 7, 15]. Dujardin et al. [5] investigated event related desynchronization in alpha band due to attention. In another research, relation between event related desynchronization in beta rhythms and emotional facial expression is explored [4]. Information obtained from these changes of EEG power are utilized in various applications such as brain computer interfacing or EEG-fMRI integration [21, 6]. The most important step of these applications is reliable EEG pattern recognition (i.e. ERD and ERS). However, studies so far have mostly used simple signal processing methods mainly based on Fourier analysis [23]. The main drawback of methods working based on conventional Fourier transform is suffering from stationary assumption and poor time-frequency localization. Studying variations in amplitude, latency or scalp distribution of ERD/ERS provides useful information about cognitive and physiological states of the brain. Therefore, estimating ERD/ERS in single trial with acceptable accuracy is essential to have a successful EEG-fMRI study. Park et. al [3] developed a technique based on multivariate empirical mode decomposition (MEMD) to extract ERD and ERS in mu and beta rhythm with the aim of motor imagery classification. They developed a method based on MEMD and common spatial pattern (CSP) to detect changes in the EEG oscillations. In another research, Ting et al.[20] proposed a technique using particle filtering for spectral estimation. They used non-gaussian state noise in order to model autoregressive (AR) parameter variation computed by Monte Carlo particle filter (PF). As an extension, a non-gaussian state space formulation of time-varying autoregressive moving average (TVARMA) model also proposed by Ting et al. [20] with the aim of improving spectral estimation. ERD/ERS detection is more challenging when EEG is acquired in MRI scanner. Although simultaneous EEG-fMRI experiments provide excellent opportunity to study brain behaviour, artifacts introduced in EEG reduce the quality of signal severely. The main destructive artifact observed in EEG-fMRI experiments is Ballistocardiogram which obscures the EEG signal mainly at mu frequencies (8 − 13Hz) and below [2]. Some artifact removal techniques are proposed by researchers [8]. Independent component analysis (ICA) or optimal basis set (OBS) are two benchmarks have been suggested for BCG removal. Performance of these techniques are very sensitive to changes of their parameters. For example, in ICA or OBS selecting and removing large number of estimated components as artifact lead to removing useful EEG information. In contrast, selecting and removing small number of estimated components as artifact lead to remain BCG. The main point is that, most of the time, difference between the large number and small number is one or two components. Therefore removing less number possible sources as BCG artifact using ICA or PCA is suggested to avoid losing useful EEG information. However, using this strategy leads to remaining BCG artifact and as a result performance of the algorithms used to extract EEG features especially in mu band will be deteriorated. Based on the above discussion, an effective method is needed to extract ERD in mu band when some parts of BCG is still exist in EEG signals. In recent years, researchers developed some algorithms based on advanced signal processing techniques for estimating ERS and ERD. For example, Ferdowsi et al. [7] proposed a method based on linear prediction (LP) to extract post-movement beta rebound (PMBR) from EEG signal recorded in simultaneous EEG-fMRI experiment. In

PT

12

CE

11

45 46 47 48 49 50 51 52 53 54 55

3

ACCEPTED MANUSCRIPT

83

2. Singular Spectrum Analysis

63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81

84 85 86

SSA is a signal processing algorithm which is able to decompose a signal into a set of underlying harmonics oscillating in various frequencies [12]. Basic SSA comprises of two complementary stages: i) decomposition and ii) reconstruction. Each stage consist of two substages.

AC

87

AN US

62

M

61

ED

60

PT

59

CE

58

CR IP T

82

another work, Ritter et al. [18] proposed a BSS-based method to separate the rolandic mu and beta rhythms from EEG signal recorded while performing a motor task. They showed that extracted rolandic mu and beta rhythms are more correlated with motor task than raw EEG signals. The highest correlation was seen in the motor cortex area which are covered by channels C3, Cz and C4. In this paper, a method based on singular spectrum analysis (SSA) is proposed for EEG spectrum decomposition. SSA is a powerful technique for time series analysis. The scope of applications of SSA is too wide comprising time series decomposition, spectral factorization, filtering and parameter estimation. Applying this method on EEG signals collected over the motor cortex area leads to estimating brain rhythms oscillating in various frequency bands. Detecting changes of power of EEG in each frequency band reveals useful information about brain function. It should be noted that, EEG signals used in this study acquired inside MRI scanner and as a result their quality deteriorated by severe artifacts. The proposed method is a multi layer approach based on SSA is considered for EEG spectral decomposition to extract event related desynchronisation in mu band. In the first layer, SSA is applied on EEG recorded by C3 channel. Then, decomposition continues until pure harmonics oscillating in mu band are estimated. In other hand, the proposed algorithm tries to extract ERD by a sequential decomposition with the aim of obtaining pure harmonics. The remained parts of BCG are removed step by step in each layer when SSA applied on EEG and its underlying components. As the EEG-fMRI data used in this work consist of movement as a task, channel C3 is selected for further analysis. In the next section, the fundamental theory of singular spectrum analysis is described. In section 3, details of the proposed strategy for mu rhythm extraction are explained. Then, the simulation results for both synthetic and real data are given in section 4 which is followed by the results of EEG-fMRI integration. Finally, the paper is concluded in Section 5.

56 57

88

2.1. Decomposition Assume that x = {x1 , x2 , . . . , xN } is a nonzero time series with length N and L(1 < L < N ) is an integer value called window length. In the first substage of decomposition, embedding or the so called trajectory matrix is constructed as follows:   x 1 x2 . . . xK  x   2 x3 . . . xK+1    x x4 . . . xK+2  X=  3   ..   .  xL xL+1 4

. . . xN

ACCEPTED MANUSCRIPT

90 91 92 93

where K = N − L + 1. In order to perform embedding, the original time series x is mapped into a sequence of lagged vectors of size L by forming K lagged vectors. In this step, choosing the appropriate window length is important as it determines the dimension of the trajectory matrix. The window length is generally selected in the range of 2 < L < N2 [12]. As it is seen, rows and columns of embedding matrix are subseries of the original time series. Note that the embedding matrix is Hankel matrix because it has equal elements on anti diagonals. In the second substage of decomposition, the singular value decomposition (SVD) of X is calculated. The ith component of X are specified by ith eigenvalue, λi , and ith eigenvector ui of XXT as illustrated in equation (1). X=

d X

98 99 100 101

102 103 104 105 106 107 108 109 110 111

2.2. Reconstruction Reconstruction consists of two steps: regrouping of SVD components, and hankelization followed by converting to a time series with length N . As it was mentioned, an appropriate criterion is needed to perform accurate grouping. This criterion can be defined based on some available prior knowledge about the original time series and its components or based on the properties of eigenvalues and eigenvectors. An important index is estimated by investigating paired eigenvalues. Each periodic time series with specific frequency produces two eigen-triples with close eigenvalues. In contrast, time series related to noise or trends produce a slowly decreasing sequence of eigenvalues. Therefore, appearing two close eigenvalues after decomposition refer to existing a pure harmonic. Let J = {j1 , j2 , · · · , jm } be a subset containing indices corresponding to the m eigenˆ J is constructed by grouping elementary values of the desired components. The matrix X matrices corresponded to eigen-triples determined by J as follows:

AC

112

M

97

where d = max{i|λi > 0} = rank(X). Since the covariance matrix C = XXT is positive-definite, its eigenvalues ,(λ1 , λ2 , · · · , λd ), are also positive. √ The triple ( λi , ui , vi ) is called eigen-triple. It should be noted that the eigenvalues are sorted in decreasing order (λ1 ≥ λ2 ≥ · · · ≥ λd ), so energy are concentrated in the most dominant components corresponded to the first eigen-triples. However, selecting dominant components is not desired in all the cases and intelligent selection strategy must be applied for successful reconstruction.

ED

96

113

114

115 116 117

(1)

PT

95

i=1

p vi = XT ui / λi

CE

94

d p X λi ui viT

AN US

i=1

Xi =

CR IP T

89

ˆ J = Xj 1 + Xj 2 + · · · + Xj m X

(2)

Assume that the set of all indices, 1, 2, · · · , d, is splitted into n disjoint subsets as J1 , J2 , · · · , Jn . The trajectory matrix of the original time series x can be obtained using the following equation. ˆ J1 + X ˆ J2 + · · · + X ˆ Jn X=X (3) ˆ J is the embedding matrix corresponded to ith subset. The procedure of selecting where X i subsets J1 , J2 , · · · , Jn is known as eigen-triple grouping. In the second step of reconstruction, the obtained trajectory matrix is converted to the form of Hankel matrix and then 5

ACCEPTED MANUSCRIPT

122

3. Mu-rhythm Extraction Using SSA

127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148

AC

149

AN US

126

M

125

Consider time series x as an EEG signal acquired from a simultaneous EEG-fMRI experiment. Applying SSA to this signal results in a number of components. The main goal of this study is to extract EEG components representing changes of neural activity in mu band (i.e. components which their frequency is between 8Hz and 12Hz). Therefore, eigen-triple grouping must be performed with the aim of selecting eigenvalues that produce components oscillating in this band. In practice, a pure harmonic with a frequency ω, certain phase and amplitude after decomposition by SSA produce two eigen-triples with very close eigenvalues. Moreover, the eigen vectors are in the form of sine waves with the same frequency [12]. Therefore, heuristic criterion based on these facts facilitates visual identification of the harmonic components. Our objective here is developing an algorithm based on the above properties in order to extract components that represent activity in mu band of EEG signals. As it was mentioned, BCG artifact distorts EEG information in mu band. Therefore, an effective technique is needed to detect changes of neural activity (i.e.ERD) in this band. In order to clarify the above claims some results of applying SSA to a segment of EEG signal are presented next. Figure 1 shows plot of 50 eigenvalues obtained when SSA is applied to a single channel EEG signal. Paired eigenvalues are distinguished using stars. Each pair produces a periodic time series with specific frequency. The √ reconstructed time √ series corresponded to one of these pairs (i.e. the eigen-triples ( λ3 , u3 , v3 ) and ( λ4 , u4 , v4 ) ) and their power spectral density (PSD) are presented in Figure 2(a-b). As it can be seen, both reconstructed time series oscillate in same frequencies (11Hz) and can be considered as pure harmonics. √ √ Reconstructed components corresponding to eigen-triples, ( λ24 , u24 , v24 ) and ( λ25 , u25 , v25 ) ) are also shown in Figure 2(d-e). As it is visible in Figure 1, values of λ24 and λ25 are not close to each other and as a results the corresponding eigen vectors are not in the form of sine waves. Power spectral density of reconstructed signals show some peaks in different frequencies which imply that these components are composed of several signals with different frequencies. In the proposed method all the reconstructed signals containing multiple harmonics are decomposed using SSA to obtain the constituent parts. This procedure continues until all the harmonics oscillating between 8 − 13Hz are extracted. Component which has the highest power is selected for further analyses (i.e. EEG-fMRI combination) . Based on the above discussion about properties of the SSA results, it is clear that some quantitative criteria are needed to identify specific eigen-triples which are generated by desired harmonics. Drawing the pairwise phase-space or scatter plots of the eigenvectors or component vectors simplifies investigation of such pairs. Figure 2(c) and (f) show the phase-space plots of eigen-triples related to Figures 2(a-b) and 2(d-e) respectively. The phase-space of pure harmonics such as sine and cosine waves with equal frequencies, amplitudes, and phases lead to the scatterplot which its points placed on a circle (Figure

ED

124

PT

123

CE

120

CR IP T

121

transformed to a time series. Assume that xpq be an element of matrix X, hankelization is performed by averaging along the elements with indices p + q = const. When hankelization is completed the obtained matrix can be converted to a time series with length N using the procedure followed by embedding step but in reverse direction.

118 119

150 151 152 153 154 155 156 157 158 159 160 161

6

ACCEPTED MANUSCRIPT

3.2 3

2.6 2.4 2.2 2 1.8 0

5

10

15

20

25

30

35

40

45

AN US

Index of eigenvalues

CR IP T

Amplitude

2.8

50

1 0 -1 0

1

2

Time (second)

1 0 -1

3

0

PSD 0

5

10

20

25

Frequency (HZ)

30

35

40

0

0

100

200

300

400

(a)

40

2

0.1

10

15

20

25

30

Frequency (HZ)

35

0.08 0.06 0.04 0.02

0

0 -0.02

40

0.04 -0.06 -0.08

0

-0.1

-0.1 -0.1

100

200

300

400

-0.05

0

0.05

0.1

u3

500

(b)

(c)

0

0

5

10

15

20

Frequncy (HZ)

25

30

35

2

3

Time (second)

0.1 0.05

0

5

10

15

20

Frequency (HZ)

25

30

35

40

−0.15

-0.1 0

-0.15

−0.15 200

(d)

300

0 -0.05

0.15

Amplitude

AC

1

20

0

40

0

100

0

40

0.15

0

0

−1

3

Time (second)

CE

PSD

5

u25

1

PT

0

80

Amplitude

Amplitude

0

Amplitude

0

0.1

3

1

1

−1

500

PSD

-0.1

0

Amplitude

0.1

15

200

M

0

Amplitude

400

200

2

Time (second)

ED

PSD

400

1

u4

Amplitude

Amplitude

Figure 1: Logarithms of 50 eigenvalues

400

500

0

100

200

(e)

300

400

500

-0.1

-0.05

0

0.05

0.1

0.15

u24

(f)

Figure 2: EEG components and their PSD; (a) From top to bottom: reconstructed time series, PSD and eiegn vector corresponded to λ3 , (b) From top to bottom: reconstructed time series, PSD and eigen vector corresponded to λ4 , (c)Phase-space of eigen vectors; u4 versus u3 , (d) From top to bottom: reconstructed time series, PSD and eigen vector corresponded to λ24 , (e) From top to bottom: reconstructed time series, PSD and eigen vector corresponded to λ25 and (f) Phase-space of eigen vectors; u25 versus u24 .

162 163

2(c)). In contrast, quasi periodic or non-periodic signals produce irregular phase-space plot (Figure 2(f)). Moreover, comparing the power spectral density of pure harmonics 7

ACCEPTED MANUSCRIPT

164 165

CR IP T

166

and quasi periodic components imply that the PSD of pure harmonics are more spiky. Therefore, kurtosis is used here as an index to measure the sharpness of the components’s PSD which are extracted by SSA. In this study, three heuristic criteria are defined to recognize eigen-triples producing pure harmonics. Equation 4 presents first criterion which is used to looking for close eigenvalues: λi (4) ∀ i, j < L : |1 − | < T r1 λj

AN US

where T r1 is threshold defined based on the application. This criterion has also been used by [10, 14] with T r1 = 0.05. Moreover, another complementary criterion based on the phase-space of eigen vectors is proposed here. Phase-space quantification with the aim of instantaneous phase analysis is the main purpose of this criterion. On the other hand, phase measurement is used here as a tool to quantifying and recognizing the pattern of phase-space map. This is a feature that has not been considered much in previous researchs about EEG. As it was mentioned, phase-space map of eigen vectors produced by pure harmonics is in the form of circle. Phase of each point in the phase-space map, θ, is computed using the following equation: ∀ i, j < L : θij = arctan (

ui ) uj

(5)

ED

M

Consider θij as a vector with length L containing calculated values of phase for all the points. Inspecting the obtained values for phase shows that θ changes periodically when the pattern of phase-space is circle wise. Therefore autocorrelation of θ is a good indication to recognize pattern of phase-space map: X |Aut(θij )| > T r2 (6)

PT

L

where P¯ is the mean of PSD, σ is the standard deviation and N is the number of data points. The PSDs of pure harmonics have higher kurtosis than the PSD of quasi periodic time series. Extracted components which their kurtosis is higher than T r3 are selected as pure harmonics, otherwise they should be decomposed into their underlying components. The optimal values of T r2 and T r3 are selected empirically. In what follows a summary of the proposed method is presented:

AC

167

CE

The third suggested criterion for discriminating pure harmonics from quasi periodic components is Kurtosis of each PSD which is computed as follows [22]: PN Pi −P¯ (7) kurtosis = i=1 4 N − 3 > T r3 σ 168

169 170

171 172

173 174 175 176 177 178

1. SSA is applied to a segment of one EEG channel with length N acquired over motor cortex area. 2. Three heuristic criteria are defined to classify principal components into two main groups including pure harmonics and components oscillating in multiple frequencies. 3. SSA is applied on all the components categorized in the group of non-pure harmonic signals. 8

ACCEPTED MANUSCRIPT

185

186

187 188 189 190 191 192 193 194 195 196

197 198 199 200 201 202 203 204 205 206 207 208

4. Experimental Results

In this section performance of the proposed method is evaluated using synthetic and real EEG data. Moreover, the results obtained by the proposed method are compared with the results of original SSA, Empirical Mode Decomposition (EMD) [17] and CLP-BSE [7]. EMD is a well known algorithm which decomposes a signal into a set of components modulated in amplitude and frequency named intrinsic mode function (IMF). Each IMF presents signal activity in a specific frequency band. However, EMD has some drawbacks such as lacking a solid mathematical formulation which leads to restriction on the performance of this method. In this study, IMFs are categorized as source of interest or noise based on their power spectrum. IMFs which their spectrum show oscillation around mu rhythms are selected as source of interest. 4.1. Synthetic data Seven signals including simulated brain rhythm in five frequency bands, ballistocardiogram artifact and noise are generated to be used as synthetic EEG sources. Sampling frequency of these signals is selected as 256Hz. Simulated EEG sources shown in Figure 3. Two first sources oscillate in 8Hz and 10Hz presenting EEG oscillation in mu band. 3rd and 6th sources presenting EEG rhythms in beta band. 4th source refers to theta band, 5th source presenting BCG artifact and the last simulated source representing the random gaussian noise. In order to generate synthetic single channel EEG signal, a linear combination of simulated sources is considered. A random 1 × 7 row vector is used for this purpose. The mixing vector has gaussian distribution and its number are in the range of [-1 1]. Figure 4 presents a sample mixture of simulated sources and its corresponding PSD. All the simulations have been implemented in MATLAB software. The proposed method is applied on synthetic EEG signals (one example is given in Figure 4(a)) with the aim of extracting mu rhythm. In order to evaluate performance of the proposed method root mean square error (RMSE) between the extracted component and the original mu rhythm is computed when different noise levels are added to EEG. Mathematical formulation of RMSE is presented as: v u N u1 X t (si − sˆi )2 (8) RM SE = 10log N 1

AC

209

CR IP T

184

AN US

183

M

182

ED

181

4. The results obtained in the 3rd step are examined in terms of qualifying criteria defined in the second step. 5. Third step is repeated until no component is remained in the category of non-pure harmonic signal. 6. Pure harmonics which their frequencies are laid in the mu band and has the highest energy among the other extracted harmonics are selected for ERD estimation and EEG-fMRI combination for further processing.

PT

180

CE

179

210 211 212

where si presents sample of the actual mu rhythm and sˆi presents sample of the extracted component as mu activity (i.e. components oscillating in 8Hz and 10Hz). In this experiment window length is chosen L = 512 which is half of signal length. Figure 5(a) presents 9

0

1

2 Time(Second)

CR IP T

AN US

Amplitude

Amplitude

ACCEPTED MANUSCRIPT

3

4

(a)

0

10

20

30 40 Frequency (HZ)

50

60

(b)

Figure 3: (a)Synthetic EEG sources (b) PSD of Synthetic EEG sources. The constituting frequency components of each source can be seen corresponding PSD.

−5

0

300 Amplitude

M

0

100 0

PT

0

1

ED

Amplitude

5

10

2 Time (Second)

3

4

(a)

20

30 Frequency (HZ)

40

50

60

(b)

213

first 10 components obtained by applying the proposed method on synthetic EEG signal shown in Figure 4. As it is seen, first four components from the top can be considered as pure harmonics. This is also confirmed by investigating numerical values of criteria defined in equations 4,6 and 7. So, the first four components are categorized in the group of pure harmonics and remaining components are categorized in the group of non-pure harmonics. As it was mentioned before, SSA applied on all the components placed in the category of non-pure harmonics. In order to decrease the number of steps of applying SSA in the third step, the proposed technique is applied on the sum of components which have similar spectrum. This can be simply done by adding components corresponded to the same singular values. Figure 2(c) presents the first 10 components obtained after performing third step of the proposed algorithm for sum of 7th and 8th components shown in Figure 5(a). As it is seen from Figure 5(c), the component presenting EEG oscillation in 8Hz is extracted successfully in the second run of the algorithm.

AC

214

CE

Figure 4: (a)An Example of synthetic EEG signal, (b) PSD of synthetic EEG signal

215 216 217 218 219 220 221 222

223 224 225

10

0

1

2 Time(Second)

3

4

CR IP T 0

10

20

30 40 frequency (HZ)

50

60

50

60

(b)

AC

Amplitude

CE

PT

Amplitude

ED

M

(a)

AN US

Amplitude

Amplitude

ACCEPTED MANUSCRIPT

0

1

2 Time(Second)

3

4

(c)

0

10

20

30 frequency (HZ)

40

(d)

Figure 5: From top to button: (a-b) extracted components by the proposed method in the first layer and their corresponding PSDs, (c-d) components extracted from signals 7th and 8th in (a) in the second layer and their corresponding PSDs.

11

ACCEPTED MANUSCRIPT

229 230 231 232 233 234 235

-30

-30

-45

-50

-40

-45

0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

-55

0

Tr1

0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

M

(a)

-40

Tr2 = 0.7

-45

-50

-50

0

Tr2 = 0.8

Tr2 = 0.07

RMSE (dB)

RMSE (dB)

RMSE(dB)

-40

Tr2 = 0.9

Tr2 = 0.8

-35

Tr2 = 0.7

-35

-35

Tr2 = 0.9

Tr2 = 0.9 Tr2 = 0.8

CR IP T

228

In order to find the optimal threshold values the following experiment is designed. Different values are selected for thresholds and performance of the proposed method is evaluated. Based on the literature and pervious experiment these threshold values are selected: 0.02 < T r1 < 0.1, 0.7 < T r2 < 0.9 and 420 < T r3 < 450. Increasing steps for Tr1, Tr2 and Tr3 are 0.01, 0.1 and 10 respectively. Figure 6 presents the obtained results when different values are selected for thresholds. Figure 6 (a), (b) and (c) show the results when T r3 is 420, 440 and 450. As it is seen, acceptable performance is obtained when T r1 < 0.06, T r2 > 0.8 and T r3 > 440. However, in order to avoid high computational cost, 0.05, 0.8 and 440 are selected for T r1 , T r2 and T r3 respectively. It should be noted that, window size is selected 512 in this experiment.

AN US

226 227

-55

0

0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

Tr1

Tr1

(b)

(c)

ED

Figure 6: Calculated RM SE for different threshold values and no noise. (a) T r3 = 420, (b) T r3 = 440 and (c) T r3 = 450 236

239 240 241 242 243 244

AC

245

PT

238

In another experiment, effect of changing window length on the performance of the proposed method is investigated. Obtained results show that selected window length is optimal to extract sources of interest. As it is seen from Figure 7, choosing small window length results in decreasing RMSE value specially for higher level of noise. Since the aim of this study is to extract mu rhythms which are low frequency harmonics, small window is not a suitable choice for extracting these rhythms. Although increasing window length leads to improving RMSE, processing time of calculating singular values and reconstruction increase. Therefor, in order to avoid increasing computational complexity and loosing low frequency oscillations half of EEG signal length is selected as optimal window size. Table 1 presents the RMSE averaged over 50 trials when level of noise is changed. In each trial, a new random vector is used to generate synthetic EEG signal. It is observed that the proposed method achieves smaller RMSE compared to other techniques.

CE

237

246 247 248 249 250

251 252 253 254 255 256

4.2. Real EEG Data The real EEG data used in this research were recorded from four healthy subjects in a simultaneous EEG-fMRI experiment. All the subjects are right handed and healthy. The study was approved by the local ethics committee of Kings College London. The EEG was acquired using the Neuroscan Maglink RT system (impedance was kept within 10-20Kohms), providing 64-channel comprising 62 scalp electrodes, one 12

ACCEPTED MANUSCRIPT

0

20

No Noise

Level of Noise (dB)

5

10

0

-20 -30

L=768 L=512 L=256 L=128 L=64

-40 -50 -60

CR IP T

RMSE (dB)

-10

Figure 7: Calculated RM SE for different window length and noise condition. Table 1: RM SE when different noise conditions exists. XCalculated XXRMSE(dB) XXX Proposed method SSA EMD EEG rhythm X X

261 262 263 264 265 266 267 268 269 270 271

AC

272

M

260

ED

259

273

274 275 276 277 278 279

280 281 282 283

-35.554 -28.076 -21.415 -16.640 -10.035

-38.147 -28.108 -20.957 -15.086 -9.064

ECG electrode and one EOG electrode. The sampling rate of raw EEG data was set at 10KHz. The experimental design used in this study consists of three blocks including cued movement, free movement, and visual control. During the cued and free movement, subjects are asked to move a joystick and during the visual control, participants should not move the joystick. In all 3 blocks, there were 9 trials per block and each trial lasted for 2.4 seconds. Algorithm proposed by Niazi et al. [16] is used to remove pulse artifact from recorded EEG. Then, in order to remove high frequency noise and baseline, a Butterworth filter with cut off frequencies 0.5Hz and 45Hz is applied to EEG signals. In the third step, EEG signals are down sampled to 250Hz and reference free recordings are obtained using the common average referencing. Eye blink and eye movement can be corrected by ICA or the proposed method. Finally, BCG artifact is removed using the method developed by Allen et al. [1]. As it was mentioned before, the aim of the proposed algorithm is to extract mu rhythm activity from EEG signals when BCG artifact partially exists. In spite of removing BCG artifact using conventional techniques, still some amounts of this artifact remained in EEG and distort meaningful information of signal, especially in mu band. Figure 8 shows a segment of EEG before and after BCG removal. Figures 8(b)(c) represent the result of applying OBS when 3 and 5 principal components are removed respectively for BCG correction. As it is seen, Ballistocardiogram is still visible in Figure 8(b). In contrast, BCG is completely removed in Figure 8(c). Although selecting and removing more principal components as BCG lead to have more clean data, useful information may also be lost along with these components. The proposed method is applied to a segment of real EEG after BCG removal using OBS algorithm. The window length and threshold values are selected as same as values used for synthetic data. Figure 9 presents the power of EEG signal in mu band when three algorithms including conventional filtering, empirical mode decomposition (EMD) and proposed method are

PT

258

CE

257

-51.835 -47.694 -41.052 -32.156 -23.739

AN US

no noise 20 10 5 0

13

ACCEPTED MANUSCRIPT

200 0 -200

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15 16

9

10 11 12 13 14 15 16

9

10 11 12 13 14 15 16

(a)

0 150

0

1

2

3

4

5

6

7

8

(b)

150 0

0

1

2

3

4

5

6

7

8

AN US

-150

CR IP T

Amplitude

150

(c)

Time (second)

Figure 8: A segment of EEG (a) before, after BCG removal using OBS (b) by removing 3 components and (c) by removing 5 components.

288 289 290 291 292 293 294 295 296

AC

297

M

287

ED

286

applied on EEG segment (Figure 8(b)). The red lines show movement triggers recorded during the experiment. Figure 9(a) shows PSD of an extracted rhythm in mu band when a band pass filter with cut off frequencies 9Hz and 11Hz is applied to extract EEG brain rhythms in mu band. Figure 9(b) shows the result obtained by EMD which is a well known technique for signal processing purposes. Figure 9(c) presents the power of components oscillating in 10Hz extracted by the proposed method when all added together. These components (i.e. 10Hz waves) are selected to investigate ERD because they have the highest energy among other extracted components in mu band. As it is clear in the Figure 9 ERD is more visible in the result obtained by the proposed method. The reason of more fluctuations in Figures 9(a-b) can be due to the fact that filters are not ideal and leakage happens around cut off frequencies. EMD alone is not able to decompose component oscillating in very close frequencies. In other word, since there is frequency overlap between the BCG artifact and brain mu rhythms, a time-frequency decomposition technique with ability to separate very close harmonics is needed. It is important to note that most of the fluctuations seen in Figures 9(a-b) are due to BCG artifact and not due to mu rhythm changes. In addition, for better demonstration, the proposed method is applied to EEG segment shown in Figure 8(c) to illustrate how EEG information are lost as a results of removing more principle components by OBS. The result is shown in Figure 9(d). As it is seen, ERD is weakened due to removing more component during BCG removal. Numerical results of this observation is given in Figure 10. As it is seen, when 5 components are eliminated for BCG removal, the ERD percentage has been decreased significantly. Table 2 shows the numerical results calculated for ERD in mu band when the above techniques are applied on EEG signal of four subjects. Moreover, CLP-BSE [7] which

PT

285

CE

284

298 299 300 301 302 303 304 305 306 307

14

ACCEPTED MANUSCRIPT

0

0

2

4

6

0

2

4

6

0

2

4

6

0

2

4

8

10

12

8

10

8

10

CR IP T

500

(a)

14

16

Amplitude

1500

0

(b)

200

0

(c)

0

14

16

12

14

16

14

16

AN US

150

12

6

8

(d)

10

12

Time (Second)

Figure 9: Power of EEG signal in mu band after applying: (a) filtering, (b) EMD, (c) proposed method on Figure 8(b); (d) proposed method on Figure 8(c).

500

OBS - 3 components OBS - 5 components

M

1000

ED

ERD %

1500

Subject1

Subject2

Subject3

Subject4

PT

0

CE

Figure 10: ERD percentage for two different BCG removal strategies.

AC

has been designed for post movement beta rebound is also applied to this data for ERD estimation. Comparing the results shows the proposed method outperforms CLP-BSE. The reason is that this method uses a reference signal which is generated by filtering the EEG signal which is noisy here and mislead the decomposition step. ERD is calculated using the following equation: ERD =

308 309 310 311 312 313

Pref erence − Pactivation × 100% Pref erence

(9)

where Pref erence is the signal power during a reference interval and Pactivation is the signal power when EEG desynchronization occurs. The obtained results presented in Table 2 show that the proposed method is able to extract ERD more effectively. Signal to noise ratio is improved as a results of eliminating the remained parts of the BCG artifact and as a consequence the calculated values for desynchronization in mu-band activity is increased. Next, we present the results of EEG15

ACCEPTED MANUSCRIPT

Table 2: Percentage of calculated value for ERD using different methods when 3 principal components are removed by OBS. PP method PP Proposed method EMD Filtering CLP-BSE Subject PP

319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337

AC

338

4.3. EEG-fMRI Combination In this section the extracted component from EEG signal demonstrating ERD is utilized to generate a time course for fMRI analysis. The most important aim is to investigate the correlation between mu rhythm desynchronization and fMRI voxels. In this research, all the fMRI analysis are performed using general linear model (GLM) which is known as a model-based technique [9]. As the experiment includes hand movements, calculated power time-course represents the instantaneous interaction between the EEG activity in mu band and motor task. In order to generate the regressor for fMRI analysis, the power of EEG time series is convolved with hemodynamic response function (HRF)[11]. GLM employes the estimated regressor to predict the BOLD response in fMRI data which are collected simultaneously with EEG. In the first stage, the collected fMRI data are preprocessed using SPM toolbox [9]. Then estimated regressor using proposed method and EMD are used to predict BOLD response in fMRI images. Figure 11 shows the estimated BOLD after fMRI group analysis. The contrast images are estimated for cued and free movements versus control. The highlighted voxels represent the areas detected as the significant variance in the fMRI images. The activation map is thresholded at p < 0.05 (Bonferroni-corrected). Figures 11(a-b) show the result of using regressor obtained by the proposed method to predict BOLD related to cued and free movement. Estimated T-scores for these movements are 3.12 and 2.96 respectively. As it is seen from Figure 11, active voxels are correlated with ERD located in motor cortex area. Figures 11(c) presents BOLD correlated with regressor obtained by EMD for cued movement with T-score 2.41. For free movement no correlated BOLD was found using EMD regressor.

AN US

318

fMRI correlation in alpha band using EEG rhythm extracted in this section.

M

317

339

340 341 342 343 344 345 346 347 348

716 584 459 721

ED

316

173 183 165 214

PT

315

201 191 184 235

CE

314

1173 1082 994 1175

CR IP T

1 2 3 4

5. Conclusion In this paper a novel method for extracting event related desynchronization in mu band is proposed. Moreover, EEG-fMRI combination is performed to find the MRI active voxels that are responsible for suppression of EEG power in mu band. EEG recorded simultaneously with fMRI are severely defected by noise. Although several algorithms have been proposed for gradient and BCG removal, some parts of these destructive artifacts are remained. In this study a method based on singular spectrum analysis is proposed to extract EEG harmonics oscillating in mu band when some parts of BCG artifact exist. The obtained results from both EEG analysis and EEG-fMRI combination show that the extracted harmonics clearly represents power suppression in mu brain rhythms. 16

ACCEPTED MANUSCRIPT

AN US

(b)

CR IP T

(a)

(c)

351 352

353 354 355

[2] G. Bonmassar, P. L. Purdon, I. P. Jskelinen, K. Chiappa, V. Solo, E. N. Brown, and J. W. Belliveau, Motion and ballistocardiogram artifact removal for interleaved recording of EEG and EPs during MRI, NeuroImage, 16 (2002), pp. 1127 – 1141.

AC

356

[1] P. J. Allen, G. Polizzi, K. Krakow, D. R. Fish, and L. Lemieux, Identification of EEG events in the MR scanner: The problem of pulse artifact and a method for its subtraction, NeuroImage, 8 (1998), pp. 229 – 239.

PT

350

References

CE

349

ED

M

Figure 11: Group analysis of EEG-fMRI combination to investigate the correlation of BOLD and ERD in mu band: (a) Cued movement (b) Free movement for proposed method, and (c) Cued movement for EMD.

357 358 359 360

361 362 363 364

365 366

[3] P. Cheolsoo, D. Looney, N. ur Rehman, A. Ahrabian, and D. Mandic, Classification of motor imagery bci using multivariate empirical mode decomposition, IEEE Transactions on Neural Systems and Rehabilitation Engineering, 21 (2013), pp. 10–22. [4] N. Cooper, A. Simpson, A. Till, K. Simmons, and I. Puzzo, Beta eventrelated desynchronization as an index of individual differences in processing human facial expression: further investigations of autistic traits in typically developing adults, Frontiers in Human Neuroscience, 7 (2013), p. 159. [5] K. Dujardin, P. Derambure, L. Defebvre, J. Bourriez, J. Jacquesson, and J. Guieu, Evaluation of event-related desynchronization (erd) during a recogni17

ACCEPTED MANUSCRIPT

373 374

375 376 377

378 379 380

381 382 383

384 385

386 387 388

389 390 391

392 393

[8] S. Ferdowsi, S. Sanei, V. Abolghasemi, J. Nottage, and O. O’Daly, Removing ballistocardiogram artifact from eeg using short- and long-term linear predictor, IEEE Transactions on Biomedical Engineering, 60 (2013), pp. 1900–1911. [9] K. J. Friston, A. P. Holmes, K. J. Worsley, J. P. Poline, C. D. Frith, and R. S. J. Frackowiak, Statistical parametric maps in functional imaging: A general linear approach, Hum. Brain Mapp., 2 (1994), pp. 189–210. [10] F. Ghaderi, H. Mohseni, and S. Sanei, Localizing heart sounds in respiratory signals using singular spectrum analysis, IEEE Transactions on Biomedical Engineering, 58 (2011), pp. 3360–3367. [11] G. H. Glover, Deconvolution of Impulse Response in Event-Related BOLD fMRI, NeuroImage, 9 (1999), pp. 416 – 429. [12] N. Golyandina, V. Nekrutkin, and A. Zhigljavsky, Analysis of Time Series Structure: SSA and Related Techniques, Chapman & Hall/CRC Monographs on Statistics & Applied Probability, CRC Press, 2001. [13] D. Krusienski, G. Schalk, D. McFarland, and J. Wolpaw, A µ-rhythm matched filter for continuous control of a brain-computer interface, IEEE Transactions on Biomedical Engineering, 54 (2007), pp. 273–280. [14] J. Mamou and E. Feleppa, Singular spectrum analysis applied to ultrasonic detection and imaging of brachytherapy seeds, The Journal of the Acoustical Society of America, 121 (2007), pp. 1790–1801.

AC

394

[7] S. Ferdowsi, S. Sanei, and V. Abolghasemi, A Predictive Modeling Approach for Joint EEG-fMRI Analysis, International Journal of Neural Systems, 25 (2015), p. 1440008.

CR IP T

372

AN US

371

M

370

[6] S. Ferdowsi, V. Abolghasemi, and S. Sanei, EEG-FMRI integration using a partially constrained tensor factorization, in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 2013, pp. 6191–6195.

ED

369

tion task: effect of attention, Electroencephalography and Clinical Neurophysiology, 86 (1993), pp. 353 – 356.

PT

368

CE

367

395

396 397

398 399 400

401 402 403

[15] Y. Meirovitch, H. Harris, E. Dayan, A. Arieli, and T. Flash, Alpha and beta band event-related desynchronization reflects kinematic regularities, Journal of Neuroscience, 35 (2015), pp. 1627–1637. [16] R. K. Niazy, C. F. Beckmann, G. D. Iannetti, J. M. Brady, and S. M. Smith, Removal of fMRI environment artifacts from EEG data using optimal basis sets, NeuroImage, 28 (2005), pp. 720 – 737. [17] J. Nunes and E. Delechelle, Empirical mode decomposition: Applications on signal and image processing, Advances in Adaptive Data Analysis, 01 (2009), pp. 125– 175. 18

ACCEPTED MANUSCRIPT

406

[18] P. Ritter, M. Moosmann, and A. Villringer, Rolandic alpha and beta eeg rhythms’ strengths are inversely related to fmri-bold signal in primary somatosensory and motor cortex, Human Brain Mapping, 30 (2009), pp. 1168–1187.

407

[19] S. Sanei and J. Chambers, EEG Signal Processing, Wiley, 2007.

408 409 410 411

412 413 414 415

416 417

[20] C. Ting, S. Salleh, Z. M. Zainuddin, and A. Bahar, Spectral Estimation of Nonstationary EEG Using Particle Filtering With Application to Event-Related Desynchronization (ERD), IEEE Transactions on Biomedical Engineering, 58 (2011), pp. 321–331.

CR IP T

405

[21] P. Wei, W. He, Y. Zhou, and L. Wang, Performance of motor imagery braincomputer interface based on anodal transcranial direct current stimulation modulation, IEEE Transactions on Neural Systems and Rehabilitation Engineering, 21 (2013), pp. 404–415.

AN US

404

[22] P. H. WESTFALL, Kurtosis as Peakedness, American statistician, 68 (2014), p. 191195.

422

Biography

423 424 425

Saideh Ferdowsi received the Ph.D. degree from the University of Surrey, Guildford, UK, in 2012. She is currently an assistant professor at the Department of Electrical Engineering, Shahrood University of Technology, Shahoord, Iran. Her main research interests include blind source separation, EEG-fMRI integration techniques, and biomedical signal and image processing.

AC

426

CE

PT

420

ED

419

M

421

[23] H. Yuan, A. Doud, A. Gururajan, and B. He, Cortical imaging of event-related desynchronization during online control of brain-computer interface using minimumnorm estimates in frequency domain, IEEE Transactions on Neural Systems and Rehabilitation Engineering, 16 (2008), pp. 425–431.

418

427 428 429

430 431 432 433 434

Vahid Abolghasemi received the Ph.D. degree from the University of Surrey, Guildford, UK, in 2011. He is currently an assistant professor at the Department of Electrical Engineering, Shahrood University of Technology, Shahoord, Iran. Before January 2013, he has worked as a Postdoctoral researcher at Brunel University, UK. His main research 19

ACCEPTED MANUSCRIPT

435 436

interests include sparse signal and image analysis, dictionary learning, compressed sensing, and blind source separation, mainly in biomedical applications.

AC

CE

PT

ED

M

AN US

CR IP T

437

20