Janis Hofmanis

Olivier Caspary

Valérie Louis-Dorr

Member, IEEE Radu Ranta

Louis Maillard

Denoising depth EEG signals during DBS using Filtering and Subspace Decomposition

Keywords: DBS, epilepsy, denoising, subspace decomposition, filtering, EEG, EMD, GEVD, SEEG, SSA, BSS

In difficult epileptic patients, the brain structures are explored by means of depth multicontact electrodes (Stereo-ElectroEncephaloGraphy, SEEG). Recently, a novel diagnostic technique allows an accurate definition of the epileptogenic zone using Deep Brain Stimulation (DBS). The stimulation signal propagates in the brain and thus it appears on most of the other SEEG electrodes, masking the local brain electro-physiological activity. The objective of this paper is the DBS-SEEG signals detrending and denoising in order to recover the masked physiological sources. We review the main filtering methods and put forward an approach based on the combination of filtering with Generalized Eigenvalue Decomposition (GEVD). An experimental study on simulated and real SEEG shows that our approach is able to separate DBS sources from brain activity. The best results are obtained by an original SSA-GEVD approach.

I. INTRODUCTION

T He precise localization of the epileptiform discharges is a crucial diagnostic step, particularly in patients considered for resective epilepsy surgery (about 20% of the epileptic patients). Invasive EEG recordings (SEEG) prove to be necessary when scalp EEG and electrical source imagery do not provide sufficient localization information [START_REF] Gavaret | Source localization of scalp-EEG interictal spikes in posterior cortex epilepsies investigated by HR-EEG and SEEG[END_REF], [START_REF] Koessler | Source localization of ictal epileptic activity investigated by high resolution EEG and validated by SEEG[END_REF]. On the other hand, electrical deep brain stimulation (DBS) is a means of eliciting the phase clustering of local neural responses [START_REF] Parra | Gamma-band phase clustering and photosensitivity: is there an underlying mechanism common to photosensitive epilepsy and visual perception?[END_REF]. Stimulationinduced modulation of a pathological network activity represents the most likely mechanism but several questions remain open to explain the DBS effects [START_REF] Mcintyre | Uncovering the mechanism(s) of action of deep brain stimulation: activation, inhibition, or both[END_REF]. Therefore, coupling DBS with SEEG measurements is really interesting for the modeling and understanding of the (epileptic) brain, for diagnosis and, also, to access the estimation of conductivities.

The functional properties of neural connections are predominantly linear [START_REF] Friston | Learning and inference in the brain[END_REF], although not exclusively (nonlinear phenomena appear mainly on the cell scale [START_REF] Breakspear | User Research Multi-stability and non-linearity of large-scale cortical rhythms[END_REF]). SEEG electrodes have a rather large surface area (1-10 mm 2 ), thus the linearity assumption holds. Considering electrical properties of the propagation medium (brain tissue), the distances between measuring and stimulation electrodes and the characteristics of the acquisition system (sampling rate, filters), one can also consider that the DBS source generates an electric field that is mixed linearly (instantaneously) with the brain sources. Therefore, this stimulation source (perfectly deterministic in time and space) propagates through the brain and masks the electro-physiological activity recorded by other depth electrodes. It is thus interesting to separate the DBS source from the brain sources (epileptic or cognitive) in order to uncover signals that can be analyzed by neurophysiologists afterwards.

The most widely employed methods belong to the linear blind source separation (BSS) family (see, e.g. [START_REF] Hyvärinen | A fast fixed-point algorithm for independent component analysis[END_REF], [START_REF] Belouchrani | Robust whitening procedure in blind source separation context[END_REF]). In (S)EEG, the instantaneous linear mixing model comes directly from the quasi-static form of the Maxwell equations [START_REF] Hild | Source localization of EEG/MEG data by correlating columns of ICAand lead field matrices[END_REF], so BSS could be applied to separate DBS from underlying physiological activities. However, BSS is based on the key hypothesis of independent sources and time invariant mixing matrix. It is not clear whether in SEEG these assumptions are valid: DBS might trigger a synchronous brain activity and the number of neural sources is unknown and fluctuating in time. Moreover, non-linear saturation effects might appear during the DBS (see section II for further details).

Besides, in DBS setups we have prior information about the stimulation source: beginning and end time, spatial localization and time structure. Consequently, the separation can be addressed by a subspace decomposition procedure, namely Generalized Eigenvalue Decomposition (GEVD) [START_REF] Sameni | A deflation procedure for subspace decomposition[END_REF], [START_REF] Gouy-Pailler | Iterative subspace decomposition for ocular artifact removal from EEG recordings[END_REF]. In GEVD, information on a particular source can be provided by different filtering methods: a frequency band, AM/FM filter bank, nonlinear filter and energy of the signal.

The aim of this work is to propose and compare different methods for the separation of multidimensional SEEG measurements during DBS. That is why we put forward an original approach based on SSA-GEVD to detrend and denoise SEEG data simultaneously. This paper is structured as follows. The signal acquisition and the DBS are presented in section II. Section III describes the model of source mixing during DBS and several proposed methods that aim at unmixing and eliminating the DBS source from measurements. The simulations and results obtained on real SEEG signals of the filtering-GEVD approach are respectively shown in sections IV and V. Finally, section VI provides our conclusions.

II. DBS-SEEG ACQUISITION

SEEG signals are acquired by multi-electrodes implanted in the brain. The implantation sites are chosen according to noninvasive data collected during the pre-surgical phase. The number of multi-electrodes for one patient varies from 8 to 13, each consisting of 5-18 measuring contacts (electrodes) (see Fig. 1). The DBS is applied between two neighboring contacts of a multi-electrode. DBS consists of periodic biphasic impulsions with intensities between 0.2 and 3 mA, lasting from 3 to 10 Fig. 1: CT image of depth electrodes implantation scheme. seconds (Fig. 2a). The frequency spectrum (Fig. 2b) of a DBS and the amplitude A k of the k-th harmonic is:

A k = 2A τ T sinc(k τ 2T ) sin(k πτ 2T ) ( 1
)
where A is the magnitude of the stimulation, τ the width of the biphasic impulsion and T the period of the impulsion. For SEEG signals, the reference has been chosen on the surface of the head at the nasion point. The signals were digitized at a sampling rate of 4096 Hz on a 128 channels amplifier (LTM 128 Headbox; Micromed, Italy) and decimated at 512 Hz. Two anti-aliasing filters have been used: a hardware filter before digitization at 4096Hz and a digital filter before decimation. These filters are designed according to the bandwidth of the EEG signal but they are not selective enough to suppress the high-frequency harmonics of the DBS. Therefore, aliasing occurs from 300 Hz harmony onwards as shown in Fig. 3, with spectral harmonics of DBS as 50, 100, 200 and 250 Hz and aliased harmonics from 300 Hz up to 1000 Hz. The formula to calculate directly n-th DBS aliased and non-aliased harmonic is:

Fn =F hs +(-1) Fn F hs (F n mod F hs )-F hs 1+ F n F hs mod 2 (2)
where F n is n th DBS harmonic, F hs -half the sampling rate, • -integer part function (also known as floor function) and (• mod •) is modulus after division. Fig. 3 shows real DBS harmonics 300, 350, 400, 450, 500, 550 Hz and so forth, aliased at 212, 162, 112, 62, 12, 38 Hz respectively.

III. METHODS AND MODEL

Considering the dipole model of K brain sources, each measured signal vector

x(t) = [x 1 (t), • • • , xM (t)]
T from the M electrodes may be written as an instantaneous linear system:

x(t) = A s(t) (3)
where A (M ×K) is the mixing matrix, s(t) the SEEG sources and n(t) the noise. The contribution of the external stimulation signal s stim (t) can be added to the previous model, such as: where s stim (t) = a s s s (t) + s trend (t) is composed of a weighted (propagated) stimulation source s s (t) and the trend s trend (t) generated by electronic commutations in the acquisition-stimulation system. The symbol • indicates possible saturated values, again due to the acquisition system limits (in practice, this saturation appears on the electrodes close to the stimulation site, when the DBS and the commutation aterfact are strong, see figure 5 below). Depending on the context, the source separation problem can be approached by a GEVD formulation, which spatially filters mixed sources using prior information. In [START_REF] Parra | Blind source separation via generalized eigenvalue decomposition[END_REF], a unified solution summarizes the assumptions for a successful estimation. A GEVD application to separate maternal ECG from foetal ECG has shown the possibility to rank the estimated sources in order of periodicity [START_REF] Sameni | A deflation procedure for subspace decomposition[END_REF]. Different types of artefacts have been removed from EEG recordings [START_REF] Gouy-Pailler | Iterative subspace decomposition for ocular artifact removal from EEG recordings[END_REF], [START_REF] Amini | MR artifact reduction in the simultaneous acquisition of EEG and fMRI of epileptic patients[END_REF]. Finally, Tome [START_REF] Tomé | Separation of a mixture of signals using linear filtering and second order statistics[END_REF] has presented a first GEVD computation with filtered signals using a Finite Impulse Response (FIR) filter.

x(t) = x(t) + s stim (t) + n(t) (4)
In this work, we propose and evaluate three original filtering-GEVD combinations. In order to make use of prior information necessary to GEVD, that is to say the close relationship between stimulation and high frequency harmonics, which are different from the low-frequency baseline and from the frequencies of electrophysiological sources, we suggested: (a) a Savitzky-Golay filter (SGF) [START_REF] Schafer | What is a Savitzky-Golay Filter?[END_REF], (b) the Singular Spectrum Analysis (SSA) [START_REF] Golyandina | Analysis of time series structure: SSA and related techniques[END_REF] and (c) the Empirical Mode Decomposition (EMD) [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. If the SGF (such as the classical FIR) is a linear frequency filtering, SSA can be considered to be a non-parametric spectral estimation in an empirical orthogonal functions space, while EMD is an empirical decomposition in AM/FM filter banks. The three filtering/GEVD combinations will be compared with a FIR/GEVD method similar to [START_REF] Tomé | Separation of a mixture of signals using linear filtering and second order statistics[END_REF], nonlinear median filtering (suggested for artefact removal in [START_REF] Amini | MR artifact reduction in the simultaneous acquisition of EEG and fMRI of epileptic patients[END_REF]) with GEVD (that we call MED/GEVD), as well as with the classical uninformed BSS.

A. Filtering approaches 1) Savitzky-Golay Filter (SGF): the SGF is frequently applied in EEG, ECG and other biomedical fields [START_REF] Hargittai | Savitzky-golay least-squares polynomial filters in ecg signal processing[END_REF] for smoothing noisy data, and it can be demonstrated that least-squares smoothing reduces noise while maintaining the shape and height of waveform peaks. It acts as a weighted average on a sliding window of 2L+1 samples. Each sample x(t) is replaced by a linear combination x SG (t):

x SG (t) = L i=-L c i x(t + i) (5)
where the SGF coefficients c i define a FIR filter depending on the order N SG of the polynomial [START_REF] Schafer | What is a Savitzky-Golay Filter?[END_REF].

2) Singular Spectrum Analysis (SSA): the Singular Spectrum Analysis [START_REF] Elsner | Singular spectrum analysis: a new tool in time series analysis[END_REF] is a non-parametric method which uses an adaptive-data basis set to extract information assumed to be stationary in the weak sense. SSA, well-adapted for short and non-stationary time series, can be summed up in three steps. The first step builds a trajectory matrix F from x(t):

F =      x(1) x(2) • • • x(L) x(2) x(3) • • • x(L + 1) . . . . . . . . . . . . x(N -L + 1) • • • • • • x(N )      (6)
The matrix F is used to estimate the correlation matrix by:

R = F T F N -L + 1 (7)
where L defines the longest periodicity that can be estimated by SSA or a multiple of the expected period in the case of trend extraction application.

Next, an eigendecomposition is applied to R, such as:

R = L i=1 λ i U i U T i ( 8
)
where λ i and U i are respectively the nonnegative decreasing eigenvalues and the column eigenvectors of R. Selecting I ≤ L eigenvectors, usually those associated to the highest eigenvalues, yields a reconstructed trajectory matrix T:

T = I i=1 U i V T i (9)
where V = F U is the decorrelated components matrix formed by the V i decomposed sources (or temporal principal components) and the matrix U the set of U i eigenvectors. The last step is a spatial smoothing thanks to an anti-diagonal averaging of T (like hankelization) to preserve their phase. Then, the filtered signal called x SSA (t) is built from each averaging value of the anti-diagonals of smoothed T.

Let T = (T i j ) L,N -L+1 i,j=1
, each sample of x SSA (t) is given by:

x SSA (t) =                1 t t m=1 T m,t-m+1 1 ≤ t ≤ L 1 L L m=1 T m,t-m+1 L + 1 ≤ t ≤ N -L 1 N -t+1 N -L+1 m=t-L+1 T t-m+1,m L + 1 ≤ t ≤ N (10)
If SSA is used to estimate DBS measurements, one of the components should correspond to non-periodic SEEG sources and to the baseline, and all others to different periodic parts (harmonics) of stimulation. Accordingly, to automatically retain SEEG sources, we need to ensure that the eigenvalues of SEEG sources are higher than that of stimulation and thus, retaining the first component after the decomposition step, will filter off the stimulation source. This can be achieved by adding a high variating artificial trend to measurements and, after the filtration, by removing it from the resulting signal. Hence, only parameter L must be adjusted (depending on the stimulation frequency) and I can be fixed to 1. The added trend is a line increasing from -1500 to 1500 microvolts.

3) Empirical Mode Decomposition (EMD): EMD is a data-based decomposition method that processes linear/non-linear, stationary/non-stationary signals. It is employed as a denoising or detrending method according to the considered application. EMD is built on the assumption that a signal can be iteratively decomposed into a finite set of different oscillation modes called Intrinsic Mode Functions (IMFs). The algorithm ends when a monotonic IMF is found to be residual. Each i th SEEG signal x i (t) is decomposed according to:

x i (t) = J j=1 IM F i,j (t) ( 11
)
where J is the number of total IMFs. The first IMFs contain high frequency components that correspond to the noise and a large part of the stimulation. For denoising, it is enough to suppress the first 2 IMFs for each x i (t) to keep low oscillating SEEG sources that are shown in [START_REF] Hofmanis | Extraction of Deep Brain Stimulation (DBS) source in SEEG using EMD and ICA[END_REF].

B. Filtering/GEVD approach

GEVD is a decomposition that requires two matrices. The first one corresponds to the correlation matrix R X built from SEEG channels X of dimension (M × N ):

R X = 1 N X X T (12)
The second one is a correlation matrix R Y calculated likewise with a matrix Y defined by Y = H(X), where the operator H(•) represents a temporal filtering function applied to each row in X. More precisely, this operator carries out a preprocessing for each SEEG channel with one of the filtering algorithms defined in section III-A. Solving GEVD for the symmetric matrix pair (R Y , R X ) implies maximizing the Rayleigh quotient defined as:

Q(E) = E T R Y E E T R X E (13)
The optimum is found by resolving a joint diagonalization:

E T R Y E = Λ E T R X E = I M (14)
The matrix E of dimension (M × M ) contains the eigenvectors and the diagonal matrix Λ the eigenvalues of (R Y , R X ) sorted in descending order. I M is the identity matrix. With GEVD, the decomposed sources are given by V = E T X. Afterwards, to eliminate unwanted components from SEEG signals, we recompose X as follows:

X = E -T G E T X ( 15
)
where the diagonal gain matrix G = diag(1, ..., 1, 0, ..., 0) selects the desired components. The filtering-GEVD approach is summarized in Fig. 4. The idea of using GEVD is to assume that, knowing the characteristics of stimulation, after filtration we will keep only desirable SEEG signals while eliminating high frequency noise and nonlinearly saturated stimulation source. GEVD diagonalizes two covariance matrices (R X and R Y ) simultaneously (Eq. 13), giving eigenvectors related to both real and filtered versions of the signal. So, components recovered by GEVD with high energy (high eigenvalues) will represent low variating sources (spanned by both R X and R Y ). In contrast, the spatial information of the stimulation is highly reduced in R Y and thus will always be rejected by joint diagonalization.

C. Blind Source Separation (BSS)

Under the hypothesis that DBS is independent of the brain activity and linearly projected on the electrodes, one might use BSS to separate it. The goal of BSS is to recover all sources (assumed independent and linearly mixed) from the observed measures, both the mixing system A and the sources S being unknown [START_REF] Cichocki | Adaptive blind signal and image processing: learning algorithms and applications[END_REF]. For each time sample of X and S the model, in its simplest form, can be written as :

x(t) = A s(t) ⇒ s(t) = B x(t) ( 16
)
where B is a reverse linear transformation of A allowing to estimate the stimulation artefact. Among BSS algorithms, we retained one based on High Order Statistics, called FastICA (Independent Component Analysis) [START_REF] Hyvärinen | A fast fixed-point algorithm for independent component analysis[END_REF], and one based on Second Order Statistics, called SOBIRO [START_REF] Belouchrani | Robust whitening procedure in blind source separation context[END_REF].

IV. SIMULATIONS AND REAL DATASETS A. Synthetic datasets

To analyze the performance of the methods described in section III, we simulated realistic datasets. A first step was the simulation of the entire SEEG acquisition system (filtering and baseline artefact) and of the stimulation pattern (Fig. 2a) as described in section II. An important side effect of the baseline deviation provoked by the commutation between acquisition and stimulation modes is the saturation observed mainly on the electrodes close to the stimulation site, where the amplitude of the propagated pattern exceeds the range of the acquisition system 1 . An example of a simulated measured DBS is presented in Fig. 5. The stimulation lasts between seconds 1 and 6, while the baseline artefact is visible during and after the stimulation.

Multichannel depth recordings (15 intracerebral electrodes) were simulated in two steps: DBS and background activities. The propagated DBS was simulated as follows: considering that the stimulation is produced between two electrodes as a dipole, we computed the potential recorded by the electrodes using the dipole field equation [START_REF] Nunez | Electric fields of the brain: the neurophysics of EEG[END_REF]:

Φ(r) ∼ = Id cos θ 4πσr 2 (17)
In this equation, d is the distance between source I and sink -I (two stimulating poles), r is the distance from the dipole center to a measurement point (electrode), θ represents the angle between the dipole axis and the vector r to an electrode and σ is the conductivity of the propagation medium. Obviously, the head is not an infinite conductor and in general three-dimensional (3-D) spatial analysis requires accurate knowledge of the electrical properties of head tissues. Nevertheless, in SEEG, as the electrodes are placed inside the brain we consider the dipole field in an isotropic homogeneous medium as good approximation. Afterwards, we generated a deviation in the stimulation signal (to produce baseline) and padded it with 5 seconds of zero activity before and after stimulation. After applying all operations of the acquisition system previously described in this section, we obtained the multidimensional stimulation signal.

The simulated background dataset consists of N simulated non-correlated EEG signals created by the non-stationary EEG model of Rankine [START_REF] Rankine | A Nonstationary Model of Newborn EEG[END_REF], adapted for SEEG by [START_REF] Caune | EEG/SEEG Signal Modelling using Frequency and Fractal Analysis[END_REF] and mixed by some randomly generated mixing matrices A M ×N , where M is the number of simulated electrodes. The result is added to the simulated stimulation signal and to the zero-mean Gaussian noise to produce the simulated DBS-SEEG: an example is presented in figure 6. The level of noise is calculated using the mean variance of simulated SEEG sources.

B. Real datasets

To check the validation of the methods, we used 1187 selected stimulations from 38 patients who underwent the SEEG acquisition/stimulation procedure at Nancy hospital, France. To reduce the amount of data, we used the electrodes which belong to the stimulating multi-electrode, because with them the stimulation is measured with more energy. Figure 7 presents the measurements of one stimulation on a multi-electrode. Fig. 7: Real SEEG stimulation. Measurements of one multi-electrode.

V. EXPERIMENTS AND RESULTS

A. Parameter selection

The objective of the filtering process before GEVD, with the methods described in section III, consists in suppressing at the same time the noise and the stimulation artefact. Therefore, the filtering efficiency will depend on the parameter choice for each method. From the further study of optimizing parameters for filtering DBS of the 50 Hz DBS frequency at the sampling rate 512 Hz, the parameters are given: SSA: L = 10 (I being 1, as explained in section III-A2), SGF: N SG = 5 and L = 20 with a cutoff frequency at about 21 Hz [START_REF] Schafer | What is a Savitzky-Golay Filter?[END_REF], EMD: J -2 last IMFs. Indeed, these parameters depend on the dynamics of the stimulation and on the high frequency components. We complete this study with a FIR filter and a median filter in the preprocessing step in order to also test their performances: a classical FIR of order 40 with a 21 Hz cutoff frequency (as for SGF), and a nonlinear median filter of sliding window length of 6 samples (chosen to be less than one period of stimulation (∼ 10 samples) but more than the length of the stimulation impulse (∼ 4 samples)) because of its intrinsic characteristics.

B. Estimation of performance

To estimate the quality of the filtering/GEVD and BSS methods in simulated data, we used a mean signal to error ratio defined as follows:

∆ = 1 M M i=1 10 log N t=1 xi (t) 2 N t=1 (x i (t) -xi (t)) 2 ( 18
)
where M is the number of measurements (electrodes), N the number of time samples, xi the ideal mixed signals (Eq. 3) for i th electrode and xi their estimate.

To measure the performance in real data, where SEEG sources are not known, we calculated the correlation index between the decomposed stimulation signal and the pattern of stimulation of 1 s. The correlation index is calculated in several steps: for every electrode i the normalized cross-correlation function is given as:

ρ i (τ ) = E[s stim,i (t)s pattern (t -τ )] σ sstim,i σ spattern ( 19
)
and then the correlation index is calculated as:

m corr = 1 M M i=1 max τ ρ i (τ ) ( 20
)
where sstim is the estimated stimulation, s pattern the stimulation pattern (a one-second long pattern of stimulation, taken from a clear, low noise measurement where no evidence of the cerebral activity was visible) and M the number of measurements.

C. Simulated SEEG data

In figure 8 we present an example of decomposed sources by SSA-GEVD, where the first component is the baseline followed by 10 simulated brain source signals (see Section IV) and the others components consists of multiple stimulation and noise mixed sources (a similar figure can be provided from other filtering/GEVD methods). Furthermore, knowing that the baseline of commutation has a high variance, it comes out as the first component of the GEVD method and thus is easy to eliminate.

For clean signal reconstruction, we need to choose the sources of interest among the estimated sources, regardless of the estimation procedure (GEVD or BSS based). In the GEVD case, the choice is straightforward, as the estimated sources are already ordered (see figures 8 and 9): the baseline is low frequency and high variance, while the stimulation is high frequency. Consequently, the first source will be the baseline and the n GEV D last sources will represent the stimulation. However, the obtained BSS components do not have a particular order and thus some additional detection steps should be carried to keep or to eliminate a component for source reconstruction. Consequently, we have ordered them according to their mean correlation to the original sources. Then, we suppressed the component which corresponds to the baseline. Thus, the first sources will represent brain activity, while the n BSS last sources (having a low correlation with the known simulated brain sources) will represent the stimulation and will be eliminated in the reconstruction. In both cases (GEVD or BSS), this number of sources needs to be evaluated. Several situations are possible:

• if the number of brain sources is known and smaller than the number of electrodes (as in Fig. 8, where we used 10 simulated brain sources, one baseline and one stimulation source for 15 electrodes), then the reconstruction will be done afterwards by using this a priori knowledge. • if the number of sources is higher (the mixture is overdetermined or noisy), then n GEV D (respectively n BSS ) needs to be estimated. An example of estimated sources for the case of 15 brain sources, one baseline and one stimulation is given in Fig. 9. In the simulation, we have evaluated the reconstruction quality using the performance index (Eq. 18) as a function of the number of eliminated sources. The results are presented in figure 10 for different noise levels (mean values over 100 simulations). To keep the figure clear, we have only displayed the results for the best GEVD based algorithm (GEVD-SSA) and the best tested BSS algorithm (FastICA). The optimal number of sources to be eliminated varies depending on the noise level, but in all cases the GEVD-SSA method overpasses FastICA. According to the performance index (see Fig. 10), we fixed n GEV D and n BSS to 2 because of the low noise condition like real SEEG (contrary EEG). In Table I, a Gaussian additive noise of different levels allows to consider the robustness of the methods on a simulated dataset. The results indicate that SSA-GEVD and FIR-GEVD have similar performances and, with MED-GEVD, outperform the other methods. These performances are more regular with SSA-GEVD and FIR-GEVD than with MED-GEVD, which is generally less efficient with a high SNR. In the case of 15 SEEG sources, there is a small diminution of performances for all methods. However, the results still remain good enough to validate the usefulness of filtering-GEVD methods in under-determined conditions except for EMD-GEVD. The low performances of classical BSS methods could be due to the non-linearity introduced by the saturation of the stimulation source on electrodes close to the stimulation. Consequently, there is no convergence of these algorithms. Figure 11 presents one example of a reconstructed signal after the elimination of the stimulation and baseline artefacts for each method, compared to input (Fig. 11: Signal+Art.) and original (Fig. 11: Signal). All methods correctly retrieve the original sources until 5 s with the retained decomposed sources as explained previously. If we examine the results of each electrode presented in figure 12 at a 20 dB noise level, the sources are less recovered at the electrodes close to the stimulation site (in our case: 10 th and 11 th electrodes). The same can be deduced for other noise levels. Obviously, it is due to the stimulation amplitude, which is saturated in these two electrodes measurements.

Finally, if we look at the spectrogram of the 8 th electrode before and after eliminating the stimulation (Fig. 13) with SSA-GEVD and FastICA, none of the frequencies caused by the stimulation or baseline is present for SSA-GEVD, while for FastICA a small amount of the stimulation frequency and its harmonics remain on the spectrogram.

D. Real SEEG analysis

The SSA-GEVD decomposed sources of a real dataset (from Fig. 7) are presented in figure 14. The baseline appears in the first decomposed source and the stimulation in the last 4 ones. As mentioned in V-B, in real datasets, we retained only the last components such as the stimulation source and suppressed everything else. Indeed, this source can be used to study the propagation in the brain and the depth-surface relationships. The performance of the methods is calculated as indicated in [START_REF] Hofmanis | Extraction of Deep Brain Stimulation (DBS) source in SEEG using EMD and ICA[END_REF]. However, we do not know how many components representing the stimulation should be eliminated. So, this number must be calculated for each dataset separately. For this, we choose a criterion based on the eigenvalues which represent the contribution of each associated eigenvector [START_REF] Konstantinides | Statistical analysis of effective singular values in matrix rank determination[END_REF] in GEVD:

J i=1 λ i M i=1 λ i 1 2 > ρ (21)
where Λ = {λ 1 , λ 2 , . . . , λ M } is the set of decreasing eigenvalues and M is the number of SEEG channels. We calculate the number of suppressed components J (J ≤ M ) for which this criterion becomes true with J = {1, 2, . . . , M } and ρ is the threshold value. Table II shows the performance of the methods on real data when this criterion is applied. Firstly, we retain M -J components where J is estimated by the criterion (Eq. 21). The threshold value ρ is empirically chosen at 0.95. Secondly, the number n of components is 2, corresponding to the last 2 components that must be kept. From table II, the separation quality of the stimulation source is very high for all filtering-GEVD methods but SSA filtering approach is slightly better. Besides, the efficiency of the filtering-GEVD methods is considerably better compared to classical BSS algorithms. The results obtained when the number of components is determined by the criterion are even slightly less good than those with a number fixed to 2. The effect of the stimulation cannot be contained in one decomposed source but can be fixed to 2. Therefore, this approach can be considered to be a good tool for neurologists to eliminate the stimulation artefact in real SEEG recordings and retrieve the underlying sources.

VI. CONCLUSION

The recording deep brain stimulation has been analyzed from a spectral point of view. The propagated stimulation signal masks the physiological signal. To retrieve the latter, we have put forward a filtering-GEVD combination which used different a priori information deduced by the filtering step and extracted by the optimization process of the GEVD. Five methods of monodimensional filtering have been implemented: SSA and EMD are data-driven decomposition processes, a non-linear median filter (MED), whereas FIR or SGF are fixed spectral filters. The point of this method is that it processes the multidimensional signals without any biophysical model. The results obtained by the mean performance index on simulated SEEG data prove the pertinence of the filtering-GEVD approach. As a result, we validate the efficiency of the filtering-GEVD based on a 38-patient SEEG database, where it shows a non-arguable better performance than the classical BSS methods. In a low-noise condition, a realistic condition for SEEG data, results are significantly improved for simulated and real datasets. The associated process SSA-GEVD provides the best performance on real data. Additionally, in the under-determined case when brain sources are more numerous than electrodes, filtering-GEVD methods are still able to filter off the stimulation source. In the physiological signal processing and particularly in the DBS procedure, the modeling is complex because we must take into account lots of variables: the interface between electrode and matter, the physical properties of different matters and the position, thus the separation of the stimulation source can be intricate. From the neurologist's point of view, this work is of real interest: it permits the analysis of underlying signals, i.e. the physiological process generated by the stimulation. In future works, besides comparing our method with other techniques adapted from neighbouring applications (see [START_REF] Taulu | Spatiotemporal signal space separation method for rejecting nearby interference in meg measurements[END_REF] for a MEG denoising application by spatiotemporal signal space separation), further improvements can be made by combining various a priori information with an optimization process.

Fig. 2 :

 2 Fig. 2: Theoretical electrical stimulation.

Fig. 3 :

 3 Fig. 3: Frequency spectrum of a real SEEG signal during DBS.

Fig. 4 :

 4 Fig. 4: Denoising multichannel signals with Filtering-GEVD.

Fig. 5 :

 5 Fig. 5: Simulated acquisition of a DBS on an electrode close to the stimulation site.

Fig. 6 :

 6 Fig. 6: Example of 15 channels simulated DBS-SEEG with a noise level of 20 dB.

Fig. 8 :

 8 Fig. 8: Decomposed sources by SSA-GEVD using a dataset with 10 mixed SEEG sources and a noise of SNR = 20 dB.

Fig. 9 :

 9 Fig. 9: Decomposed sources by SSA-GEVD using a dataset with mixed 15 SEEG sources and a noise of SNR = 20 dB.

Fig. 10 :

 10 Fig. 10: Comparing the performance index for a different number of eliminated sources (n GEV D and n BSS ) and different noise levels, in the case of 15 brain sources.

Fig. 11 :

 11 Fig. 11: Example of denoised mono-channel signal from a synthetic dataset of 15 SEEG sources without noise. Signal+Art: simulated mono-channel SEEG signal with trend and stimulation artefacts. Signal: simulated mono-channel SEEG signal without artefacts. Following by denoised signals with the presented methods.The best extraction of artefacts is obtained with SSA-GEVD and FIR-GEVD followed by SGF-GEVD where a slight trend of baseline remains in signal around 6.5 s. The EMD-GEVD and MED-GEVD keep the residue of stimulation artefact but EMD-GEVD has the worst baseline extraction. BSS algorithms partly keep the trace of baseline and reproject more of the stimulation artefact.If we examine the results of each electrode presented in figure12at a 20 dB noise level, the sources are less recovered at the electrodes close to the stimulation site (in our case: 10 th and 11 th electrodes). The same can be deduced for other noise levels. Obviously, it is due to the stimulation amplitude, which is saturated in these two electrodes measurements.Finally, if we look at the spectrogram of the 8 th electrode before and after eliminating the stimulation (Fig.13) with SSA-GEVD and FastICA, none of the frequencies caused by the stimulation or baseline is present for SSA-GEVD, while for FastICA a small amount of the stimulation frequency and its harmonics remain on the spectrogram.

Fig. 12 :Fig. 13 :Fig. 14 :

 121314 Fig. 12: Mean performance index for each electrode at a 20 dB noise level using the presented methods.

Fig. 15 :

 15 Fig. 15: Elimination of the stimulation and baseline artefacts from SEEG data with an epileptic activity during DBS (DBS ends at 3.1 s). (a) -real SEEG as registered with the epilepsy during DBS. (b) -sources after SSA-GEVD decomposition. (c) reconstructed signal eliminating the first and the last two components.

Figure 15 -

 15 Figure15-a illustrates an epileptic crisis activity during the stimulation for a patient. After decomposition with SSA-GEVD, in figure15-b, parts of the stimulation source are clearly visible over the last two components and the baseline is present in the first. After reconstruction, the epileptic activity is recovered with a little stimulation noise (see Fig.15-c).

TABLE I :

 I Comparison of filtering-GEVD, FastICA and SOBIRO algorithms using simulated data with different levels of additive noise. Mean performance index and mean standard deviation in brackets for all 15 electrodes after 1000 trials.

			10 sources			15 sources
				SNR in dB
	Method	5	10	20	5	10	20
	SSA-GEVD	4.60 (1.90)	7.96 (2.04)	10.85 (2.79)	4.04 (1.62)	7.45 (1.73) 10.39 (1.63)
	FIR-GEVD	4.57 (2.03)	8.02 (2.61)	10.96 (1.92)	4.01 (1.73)	7.40 (2.20)	10.30 (2.33)
	SGF-GEVD	4.47 (2.61)	7.66 (1.94)	9.20 (2.08)	3.95 (2.20)	7.20 (1.65)	9.13 (1.75)
	EMD-GEVD	4.14 (2.11)	4.67 (2.80)	6.31 (2.01)	2.87 (1.83)	5.54 (1.52)	6.43 (1.66)
	MED-GEVD	4.50 (1.92)	7.61 (2.08)	11.04 (2.56)	3.98 (1.64)	7.17 (1.75)	10.25 (2.17)
	FastICA	2.72 (2.00)	4.11 (2.99)	5.00 (3.13)	2.82 (1.71)	4.23 (2.58)	4.62 (2.48)
	SOBIRO	-1.28 (2.47)	0.95 (3.21)	3.46 (2.71)	-0.26 (2.14)	2.42 (2.55)	4.96 (2.33)

TABLE II :

 II Comparison of filtering-GEVD, FastICA and SOBIRO algorithms using real datasets (1187 DBS multichannel recordings). Index of correlation and standard deviation in brackets for n fixed and calculated by the criterion.

	Method	m corr (ρ = 0.95) m corr (n = 2)
	SSA-GEVD	0.945 (0.069)	0.95 (0.065)
	FIR-GEVD	0.942 (0.075)	0.946 (0.073)
	SGF-GEVD	0.91 (0.113)	0.944 (0.06)
	EMD-GEVD	0.89 (0.146)	0.93 (0.099)
	MED-GEVD	0.918 (0.095)	0.898 (0.13)
	FastICA	X	0.81 (0.114)
	SOBIRO	X	0.822 (0.117)

Consequently, as some canals are saturated differently from others, the mixture becomes non-linear and it leads to difficulties for linear signal processing methods such as BSS.

This study was supported by the French Ministry of Health (PHRC 17-05, 2009) and CHU de Nancy.