Fabien Lotte
email: fabien.lotte@gmail.com

Anatole L Écuyer
email: anatole.lecuyer@irisa.fr.

Bruno Arnaldi
email: arnaldi@irisa.fr.

A Lécuyer

FuRIA: An Inverse Solution based Feature Extraction Algorithm using Fuzzy Set Theory for Brain-Computer Interfaces

Keywords: brain-computer interface (BCI), feature extraction, inverse solution, fuzzy sets, electroencephalography (EEG) EDICS Category: Biomedical Signal Processing (BIO) I

This paper presents FuRIA, a trainable feature extraction algorithm for non-invasive Brain-Computer Interfaces (BCI). FuRIA is based on inverse solutions and on the new concepts of fuzzy Region Of Interest (ROI) and fuzzy frequency band. FuRIA can automatically identify the relevant ROI and frequency bands for the discrimination of mental states, even for multiclass BCI. Once identified, the activity in these ROI and frequency bands can be used as features for any classifier. The evaluations of FuRIA showed that the extracted features were interpretable and can lead to high classification accuracies.

 [START_REF] Baillet | Electromagnetic brain mapping[END_REF]. Hence, they may localize the neural sources of activity within the brain, thus recovering a physiologically relevant information. Recently, inverse solutions were shown to be promising feature extraction methods for EEG-based BCI [START_REF] Lotte | Furia: A novel feature extraction algorithm for brain-computer interfaces using inverse models and fuzzy region of interest[END_REF].

In this paper we propose a trainable feature extraction algorithm for BCI which relies on inverse solutions. This algorithm is called FuRIA which stands for Fuzzy Region of Interest Activity. FuRIA aims at automatically identifying what are, for a given subject, the relevant Regions Of Interest (ROI) and their associated frequency bands for the discrimination of mental states. The activity in these ROI and associated frequency bands can be used as features for any classifier.

The paper is organized as follows: Section II details the concept of inverse solution and reviews its use for BCI design. Section III describes in details the FuRIA algorithm we propose. Section IV presents the results of the FuRIA evaluations. Finally, Section V discusses the results and concludes.

II. BACKGROUND: INVERSE SOLUTIONS AND BCI

A. Inverse solutions

When using EEG, the signals m(t) (m ∈ ℜ Ne, 1 with N e being the number of electrodes used) recorded at time t on the scalp can be modeled by a linear combination of brain dipole activity c(t) (c ∈ ℜ 3×Nv,1 with N v being the number of dipoles considered). This is called the forward problem [START_REF] Baillet | Electromagnetic brain mapping[END_REF], which can be modeled by:

m(t) = Kc(t) (1)
where K is an N e × (3 × N v ) matrix called the leadfield matrix which represents the physical properties (conduction) of the head. This matrix is a head model in which each dipole is modeled by a volume element known as a voxel. Typical head models are composed of thousands of voxels [START_REF] Baillet | Electromagnetic brain mapping[END_REF]. The c(t) vector holds the orientation and amplitude of each dipole, according to the three dimensions of the head model space.

Inverse solutions aim at estimating the brain dipole activity ĉ(t) by using only the scalp measurements m(t) and the leadfield matrix (head model) K, which can be modeled by:

ĉ(t) = T m(t) (2)
where T is the generalized inverse of K. As N v >> N e , this problem has no unique solution and additional constraints must be added to solve it. Depending on the constraints used, different inverse solutions are obtained which leads to different T matrices [START_REF] Michel | Eeg source imaging[END_REF]. There are two main kinds of inverse solutions:

distributed solutions and equivalent dipole solutions [5] [3]. Distributed solutions estimate the amplitudes and orientations of a large number of voxels distributed in all the cortex or in all the brain whereas equivalent dipole solutions estimate the position, amplitude and orientation of few sources (typically one or two), each one modeled by an equivalent dipole.

Congedo has shown that with any linear and discrete inverse solution the amplitude γ v (also known as current density) in voxel v could be formulated as a quadratic form [START_REF] Congedo | Subspace projection filters for real-time brain electromagnetic imaging[END_REF]:

γ v (t) = m(t) T Q v m(t) (3)
Here, Q v is an N e ×N e matrix denoted as the inverse operator for voxel v [START_REF] Congedo | Subspace projection filters for real-time brain electromagnetic imaging[END_REF]. The superscript T denotes transpose. Typically Q v = T T v T v with T v being the v th row of T . Such a γ v measure has been used as a feature for BCI (see next section). The current density in a given Region Of Interest (ROI) Ω, i.e., in a set of voxels, can be computed as follows:

γ Ω (t) = v∈Ω γ v (t) = m(t) T Q Ω m(t) with Q Ω = v∈Ω Q v (4)
This notation is very convenient as it allows very fast computations, whatever the number of voxels in the ROI Ω. Actually the Q Ω matrix can be computed offline, and the size of this matrix depends only on the number of electrodes used.

B. A review of inverse solution-based BCI

Concerning BCI, inverse solutions extract physiologically relevant information which appears as an attractive type of feature. Recently, a few studies have focused on inverse solutions for feature extraction and have obtained promising results [START_REF] Lotte | Furia: A novel feature extraction algorithm for brain-computer interfaces using inverse models and fuzzy region of interest[END_REF] [7] [START_REF] Kamousi | Classification of motor imagery tasks for brain-computer interface applications by means of two equivalent dipoles analysis[END_REF]. In order to design BCI, inverse solutions are generally used in two different ways:

• As a direct feature extraction technique: in this case, either the brain current density computed in a number of ROI [START_REF] Congedo | Classsification of movement intention by spatially filtered electromagnetic inverse solutions[END_REF] or the neural source positions [10] [8] are used as features to identify a mental state. • As a preprocessing method preceding feature extraction: in this case, the inverse solution is used to estimate c(t) from which the features are extracted [START_REF] De Peralta Menendez | Non-invasive estimation of local field potentials for neuroprosthesis control[END_REF] [12] [START_REF] Noirhomme | Single trial eeg source reconstruction for brain-computer interface[END_REF]. These methods have all obtained very satisfying results. Moreover, it has been observed that extracting features from c(t) (the source domain) is more efficient than extracting them directly from m(t) (the sensor domain) [START_REF] De Peralta Menendez | Non-invasive estimation of local field potentials for neuroprosthesis control[END_REF] [7] [START_REF] Besserve | Non-invasive classification of cortical activities for brain-computer interfaces: a variable selection approach[END_REF].

In spite of these promising results, some limitations remain. Indeed, current methods are either general-purpose, i.e., they have the ability to deal with any kind of mental state, or they generate a consise set of discriminative features but rarely both at the same time. Several methods require strong a priori knowledge on the neurophysiological mechanisms involved by the mental states used, and hence, are not generalpurpose at all [START_REF] Qin | Motor imagery classification by means of source analysis for brain computer interface applications[END_REF] [8] [START_REF] Babiloni | The estimation of cortical activity for brain-computer interface: Applications in a domotic context[END_REF]. With these methods, the ROI to be used must be defined beforehand and by hand. These methods are currently limited to the use of mental states that involve the motor and sensorimotor areas of the brain. Moreover, a non-general-purpose method, exclusively based on a priori knowledge, will not be able to adapt to each subject's specificities, and will have, most probably, nonoptimal performances.

A few general-purpose methods, based on distributed inverse solutions, have been proposed [START_REF] De Peralta Menendez | Non-invasive estimation of local field potentials for neuroprosthesis control[END_REF] [7] [START_REF] Besserve | Non-invasive classification of cortical activities for brain-computer interfaces: a variable selection approach[END_REF]. Their main drawback is that they must extract one or several features per voxel, which generates a very large number of features. Thus, this requires the use of feature selection techniques [START_REF] De Peralta Menendez | Non-invasive estimation of local field potentials for neuroprosthesis control[END_REF] [7] [START_REF] Besserve | Non-invasive classification of cortical activities for brain-computer interfaces: a variable selection approach[END_REF]. Even though this solution gives good results, the number of features used remains generally relatively high, particularly in comparison with the number of features extracted by nongeneral-purpose methods [START_REF] Qin | Motor imagery classification by means of source analysis for brain computer interface applications[END_REF] [START_REF] Kamousi | Classification of motor imagery tasks for brain-computer interface applications by means of two equivalent dipoles analysis[END_REF]. Moreover, in these methods, all voxels are processed independently, whereas the current densities in neighboring voxels are generally correlated. As such, these voxels should be gathered in brain regions.

Congedo et al have proposed a method which is generalpurpose and generates a small set of discriminative features, which are determined after gathering the voxels into a few non-overlapping ROI [START_REF] Congedo | Classsification of movement intention by spatially filtered electromagnetic inverse solutions[END_REF]. However, this method still needs improvements as it is not completely automatic and currently limited to the use of two ROI [START_REF] Congedo | Classsification of movement intention by spatially filtered electromagnetic inverse solutions[END_REF].

III. THE FURIA FEATURE EXTRACTION ALGORITHM

FuRIA is an inverse solution-based algorithm which can learn and use subject-specific features for mental state classification. A feature extracted with FuRIA corresponds to the activity in a given ROI and its associated frequency band. Contrary to existing methods, FuRIA can automatically identify relevant ROI, as well as the frequency bands in which these ROI current densities are discriminant. Finally, FuRIA also introduces the concepts of fuzzy ROI and fuzzy frequency bands to obtain increased classification performances.

FuRIA aims at being modular in the sense that various kinds of inverse solutions could be used within it. This section briefly describes the inverse solutions that could be used within FuRIA and the specific one that we used in our implementation. It then describes in details the FuRIA feature extraction algorithm.

A. Inverse solutions for FuRIA

Equivalent dipole solutions have been popularly employed for BCI as they can obtain satisfactory results with only two features and without using a trained classifier [START_REF] Qin | Motor imagery classification by means of source analysis for brain computer interface applications[END_REF] [START_REF] Kamousi | Classification of motor imagery tasks for brain-computer interface applications by means of two equivalent dipoles analysis[END_REF]. However, FuRIA is based on linear and distributed inverse solutions. Actually, distributed solutions enable the use of a large number of dipoles rather than a few equivalent dipoles. As such they provide more information and are more flexible.

On the other hand, the use of linear inverse solutions appears as essential for BCI applications. Indeed, the computational demand of non-linear inverse solutions makes them inappropriate for online BCI operation. Several linear and distributed inverse solutions have been used for BCI, such as ELECTRA [START_REF] De Peralta Menendez | Non-invasive estimation of local field potentials for neuroprosthesis control[END_REF], sLORETA (standardized low resolution electromagnetic tomography) [START_REF] Congedo | Classsification of movement intention by spatially filtered electromagnetic inverse solutions[END_REF] or the depth-weighted minimum norm technique [START_REF] Babiloni | The estimation of cortical activity for brain-computer interface: Applications in a domotic context[END_REF]. FuRIA can be used with any of them.

In our implementation of FuRIA, we used sLORETA which is an instantaneous, discrete and linear inverse solution [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sloreta): technical details[END_REF]. sLORETA is known to have very good localization properties [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sloreta): technical details[END_REF] including no localization bias in the presence of measurement and biological noise [START_REF]Discrete, 3d distributed, linear imaging methods of electric neuronal activity. part 1: exact, zero error localization[END_REF]. Moreover, it has been proven experimentally that sLORETA is suitable for the design of EEG-based BCI [START_REF] Congedo | Classsification of movement intention by spatially filtered electromagnetic inverse solutions[END_REF]. With sLORETA, the inverse operator

Q v is computed by Q v = T T v S -1 v T v with S v = T v K v ,
T being obtained with the minimum-norm method [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sloreta): technical details[END_REF] [START_REF] Congedo | Subspace projection filters for real-time brain electromagnetic imaging[END_REF]. More details on sLORETA can be found in [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sloreta): technical details[END_REF] [START_REF]Discrete, 3d distributed, linear imaging methods of electric neuronal activity. part 1: exact, zero error localization[END_REF]. Even if we used sLORETA, it should be stressed that the principle of the FuRIA algorithm is not dependent on the inverse solution chosen. Consequently, any other distributed and linear inverse solution could be used instead of sLORETA.

B. Overview of the FuRIA algorithm 1) Training of FuRIA:

To be used, FuRIA has first to be trained using a set Θ = {(m(t), C) 1..Nt } of labelled training data. C is the class of m(t), i.e., the mental state of the subject while m(t) was recorded. This training phase aims at finding subject-specific ROI Ω l and frequency bands Φ l that contain the most relevant information for mental state discrimination. This phase is accomplished offline, in three steps:

• Identification of statistically discriminant voxels and frequencies: The goal of this step is to identify the ordered pairs

w k = (f i , v j ) (k ∈ [1, N k ]
), with f i being a frequency and v j being a voxel, with the largest discriminative power. In order to do so, we rely on a statistical test for comparing, between the different classes, the mean current density in each frequency f i (i ∈ [1, N f ]) and in each voxel v j (j ∈ [1, N v ]).

• Creation of ROI and frequency bands: This step aims at clustering the ordered pairs w k selected in the previous step into a smaller number of ordered pairs W l (l ∈ [1, N w ]). Here W l = (Φ l , Ω l ), with Φ l being a frequency band, i.e., a set of frequencies, and Ω l being a ROI, i.e., a set of voxels. In order to do so, we rely on a clustering algorithm for finding clusters of voxels and frequencies and for transforming these clusters into ROI associated to frequency bands. Thus, by the end of this step, we have created a set of N w ordered pairs W l (one pair per cluster) in which each ROI is associated to a single frequency band. This frequency band should gather the frequencies in which the activity of the ROI voxels is discriminant. • Fuzzification of ROI and frequency bands: The found ROI Ω l are turned into fuzzy ROI Ωl and the frequency bands Φ l into fuzzy frequency bands Φl . This aims at giving more importance to the more discriminant voxels and frequencies, while still using the information contained in the less discriminant ones. The overall objective is to increase the discriminative power of the W l pairs.

2) Use of FuRIA for feature extraction: Once the fuzzy pairs Wl = ( Φl , Ωl ) have been identified, FuRIA can be used for feature extraction.

The l th feature extracted is the average current density in the fuzzy ROI Ωl after band-pass filtering EEG signals in the associated fuzzy frequency band Φl . All these steps are detailled hereafter.

In this paper, we evaluate a specific implementation of FuRIA. However, it is worth noting that other implementations could be used (e.g., different clustering algorithms or different statistical tests) as long as they are consistent with the algorithm proposed.

C. First training step: identification of statistically discriminant voxels and frequencies 1) Algorithm: The first training step of FuRIA aims at identifying the pairs w k of voxels v j and frequencies f i which are the most discriminant, i.e., the pairs of voxels and frequencies whose current density is the most different between classes. To do so, we perform a statistical test that compares the mean current density between classes for each pair w k = (f i , v j ). To this end, each training EEG record m(t) passes through the following procedure (see Fig. 1):

1) m(t) is decomposed into frequency bands by using a set of filters h i . Each filter h i is a 2-Hz wide band-pass filter centered on frequency f i . We denote as m i (t) the signal resulting from the filtering of m(t) by h i . 2) the current density γ i,j (t) in voxel v j , for frequency f i is computed using the inverse solution:

γ i,j (t) = m i (t) T Q vj m i (t) (5)
3) finally, γ i,j (t) is averaged over a time window of interest, which starts at sample t 0 and is N s samples long. The obtained value is then log-transformed:

< γ i,j >= log 1 N s t0+Ns t=t0 γ i,j (t) (6)
Fig. 1. Computation of the average current densities < γ i,j > in all frequencies f i and voxels v j , from a training data m(t).

All < γ i,j > are then gathered into statistical samples according to the label of their corresponding training record m(t). These samples are compared using the statistical test. In other words, this statistical test compares the mean < γ i,j > between classes and thus gives the discriminative power of each pair w k . Pairs w k which obtained a p-value higher than a given (user-selected) threshold α are not considered anymore in the remaining of the training process. The other pairs are denoted as "significant". This procedure should remove numerous voxels and frequencies and should only keep the ones that are related to the mental states considered. It should be noted that the w k selected can be different from one subject to the other.

2) Implementation: In our implementation of FuRIA, the h i filters were either Finite Impulse Response (FIR) or Infinite Impulse Response (IIR) filters. We used windowed sinc filters as FIR filters [START_REF] Smith | The Scientist & Engineer's Guide to Digital Signal Processing[END_REF] and Yule-Walker filters as IIR filters [START_REF] Friedlander | The modified yule-walker method of arma spectral estimation[END_REF]. These filters were chosen as they enable the design of custom filters, which is needed for the fuzzification step (see section III-E). Concerning the statistical test used, we employed Multiple Comparisons Randomization (MCR) tests as described by Holmes [START_REF] Holmes | Nonparametric analysis of statistic images from functional mapping experiments[END_REF]. More precisely, for BCI with only 2 classes, we used MCR t-tests whereas for multiclass BCI (with 3 or more classes) we used MCR ANOVA (ANalysis Of VAriance), as ANOVA can compare multiple conditions, and as such, can deal with more than 2 classes.

D. Second training step: creation of ROI and frequency bands 1) Algorithm: This step aims at gathering pairs of significant voxels and frequencies into pairs of ROI

Ω l (l ∈ [1, N w ])
and frequency bands Φ l , each ROI being associated to a single frequency band. Basically, a given ROI would gather significant voxels, and the frequency band associated to this ROI would gather the frequencies in which these voxels activity is discriminant. Creating such ROI and frequency bands aims at obtaining a compact feature set. Indeed, using the activity in a few ROI and frequency bands as features should lead to much fewer features than when considering voxels and frequencies alone, as generally done [START_REF] De Peralta Menendez | Non-invasive estimation of local field potentials for neuroprosthesis control[END_REF] [START_REF] Noirhomme | Single trial eeg source reconstruction for brain-computer interface[END_REF]. Moreover, the activities in neighboring voxels and frequencies tend to be statistically correlated [START_REF] Michel | Eeg source imaging[END_REF]. Thus, it may be appropriate to use these voxels and frequencies together rather than independently.

We would like to gather voxels belonging to the same neural source into the same ROI. We would also like to gather into a single frequency band the frequencies in which each ROI is discriminant. Thus, it is desirable to find clusters gathering both voxels and frequencies. To do so, we associate to each significant pair w k = (f i , v j ) the feature vector [x j , y j , z j , f i ] in which x j , y j , z j are the spatial coordinates of v j . Each one of these 4 coordinates was normalized to zero mean and unit variance over all j in order to deal with the different ranges between space and frequencies. Then, we apply a given clustering algorithm to all these vectors. Finally, for each one of the N w clusters obtained automatically, we gather into the same ROI Ω l all the voxels whose spatial coordinates correspond to the x j , y j , z j coordinates of one of the vectors from this cluster. We also associate to this ROI Ω l the frequency band Φ l = [f min , f max ]. Here, f min and f max are respectively the minimal and maximal value of the coordinate f i among all the vectors belonging to this cluster. This clustering gives a set of ordered pairs W l = (Φ l , Ω l ) which are expected to be discriminant.

2) Implementation: When using sLORETA, the neural sources tend to appear as local maximums of the current density [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sloreta): technical details[END_REF]. Hence, we assume they would also appear as local maximums of the statistics obtained in output of the statistical test. Even though it is only an assumption, this seems to be often verified in practice (see results in Section IV-C for instance). We used Mean Shift as the clustering algorithm since it gathers vectors attracted by the same local maximum of the underlying density function [START_REF] Comaniciu | Mean shift -a robust approach toward feature space analysis[END_REF]. However, as the voxels coordinates and frequencies considered are regularly spaced within their numerical domain, the underlying density function for the vectors [x j , y j , z j , f i ] will be relatively flat, and thus will prevent a proper use of Mean Shift. To cope with this problem, we used a slightly modified version of the Mean Shift algorithm for clustering. This slight modification simply consists in replacing D, the standard density estimate at point P used in Mean Shift:

D(P ) = 1 N k H d N k k=1 χ 1 H (P -P k ) (7)
by D, a weighted density estimate at point P :

D(P ) = 1 N k H d N k k=1 s k • χ 1 H (P -P k ) (8)
with H being the smoothing parameter, χ a kernel function, here the Epanechnikov kernel (which was found to be the optimal kernel for Mean Shift [START_REF] Comaniciu | Mean shift -a robust approach toward feature space analysis[END_REF]), P the current vector, P k , the k th vector from the data set, and d the dimensionality.

Finally s k = 1-p k , with p k being the p-value obtained by the w k pair, during the statistical test performed in the previous step. The vector P k is the vector corresponding to the w k pair. This leads to the following form for the sample mean shift vector :

M H (P ) = 1 n P Pi∈SH (P ) s k • (P k -P ) (9)
Here, S H (P ) is the d-dimensional sphere of radius H centered at P , with n P vectors inside. This weighted version of Mean Shift will gather into the same cluster all the vectors attracted by the same local maximum of the statistics. As such, this version of Mean Shift aims at gathering altogether the vectors corresponding to the same neural source, as these sources are assumed to be local maximums of the statistics (see above).

E. Third training step: fuzzification of ROI and frequency bands

The last training step of FuRIA consists in fuzzifying the previously obtained ROI Ω l and frequency bands Φ l . Actually, a ROI can be seen as a conventional (or "crisp") sets of voxels whereas a frequency band can be seen as a crisp set of frequencies. However, in a ROI or in a frequency band, all the voxels and frequencies do not have the same discriminative power. Nevertheless, all these elements still carry information that could be used, making it hard to choose which of them should be kept. Consequently, we believe that all significant voxels and frequencies should be used, but the voxels and frequencies that are less discriminant should belong "less" to their ROI and frequency band than the others. Thus, we propose to consider ROI and frequency bands as fuzzy sets [START_REF] Zadeh | Fuzzy sets," Fuzzy sets, fuzzy logic, and fuzzy systems: selected papers by[END_REF] of voxels and frequencies, in which all voxels and frequencies are given a degree of membership into the ROI or frequency band to which they belong. We denote such kinds of ROI and frequency bands as fuzzy ROI and fuzzy frequency bands. Thus, a fuzzy membership function µ is associated to each ROI and each frequency band.

1) Design of fuzzy ROI and fuzzy frequency bands from a given fuzzy membership function:

a) Algorithm: As mentionned above, a crisp ROI Ω l is a set of voxels which current density is computed according to Eq. 4. A fuzzy ROI Ωl is not defined by a set of voxels anymore but by a fuzzy membership function µ Sl . This function provides the degree of membership, in [0, 1], of any existing voxel to the fuzzy ROI Ωl . Contrary to crisp ROI for which the activity in all its voxels is used equally (see Eq. 4), we define the current density γ Ωl (t) in a fuzzy ROI as follows:

γ Ωl (t) = Nv j=1 µ Sl (v j )γ vj (t) (10)
This leads to:

γ Ωl (t) = m(t) T Q Ωl m(t) with Q Ωl = Nv j=1 µ Sl (v j )Q vj (11)
Similarly, we associate a fuzzy membership function µ F l to each frequency band Φ l . The function µ F l provides the degree of membership, in [0, 1], of any existing frequency to the fuzzy frequency band Φ l . We can note that this function has exactly the same form as the magnitude response of a digital filter. This means that to band-pass filter a signal in a given fuzzy frequency band we have to design a custom digital filter that has the desired fuzzy membership function µ F l as magnitude response.

b) Implementation: We used the window technique to automatically design FIR filters from the desired magnitude response, i.e., from a fuzzy membership function to a fuzzy frequency band [START_REF] Smith | The Scientist & Engineer's Guide to Digital Signal Processing[END_REF]. To design the IIR filters, we used the Yule-Walker method [START_REF] Friedlander | The modified yule-walker method of arma spectral estimation[END_REF]. This explains why we used a windowed-sinc FIR filter or a Yule-Walker IIR filter for the first training step of FuRIA (see section III-C).

2) Set up of the fuzzy membership functions:

a) Algorithm: To determine the kind of fuzzy membership functions to be used as well as their parameters, we first compute the discrimination scores d vj and d fi of each voxel v j and frequency f i respectively, for each pair W l = (Φ l , Ω l ):

d vj = 1 N f fi∈Φ l s vj ,fi and d fi = 1 N v vj ∈Ω l s vj ,fi (12)
where s vj ,fi = 1p j,i with p j,i the p-value obtained during the first training step, for voxel v j at frequency f i . In order to highly emphasize the contribution of the most discriminative voxels and frequencies, which are generally much less numerous than others, we chose exponential fuzzy membership functions that are Gaussian:

µ S l (v j ) = exp(-1 2 ( dv j -dvmax σv ) 2 ) v j ∈ Ω l 0 otherwise ( 13
)
µ F l (f i ) = exp(-1 2 ( d f i -d f max σ f ) 2 ) f i ∈ Φ l 0 otherwise ( 14
)
where d vmax and d f max are the maximal scores among voxels and frequencies respectively. The values of the σ should be optimized in order to maximize the discriminative power of the pairs Wl = ( Φl , Ωl ). Consequently, we use

σ v = 1 card(Ω l )
and σ f = 1 card(Φ l ) as initial values, and then optimize these values using the adaptive gradient ascent procedure described in Algorithm 1. Here λ f and λ v are positive learning rates,

Algorithm 1 Adaptive gradient ascent algorithm (ǫ f ,ǫ v ,λ f ,λ v ) 1: newF ← F (σ v , σ f ) 2: repeat 3: oldF ← newF 4: ∆F f ← F (σ v , σ f ) -F (σ v , σ f + ǫ f ) 5: σ f ← σ f -λ f ∆F f ǫ f {σ f update} 6: newF ← F (σ v , σ f ) 7: if newF < oldF then 8: λ f ← λ f 2 {λ f Adaptation} 9: else 10: λ f ← λ f + 0.1λ f {λ f Adaptation} 11: end if 12: oldF ← newF 13: ∆F v ← F (σ v , σ f ) -F (σ v + ǫ v , σ f ) 14: σ v ← σ v -λ v ∆Fv ǫv {σ v update} 15: newF ← F (σ v , σ f ) 16:
if newF < oldF then 17:

λ v ← λv 2 {λ v Adaptation} 18:
else 19:

λ v ← λ v + 0.1λ v {λ v Adaptation} 20: end if 21: until | ∆F f ǫ f | ≤ 0.001 and | ∆Fv ǫv | ≤ 0.001
and ǫ v and ǫ f are small positive increments used to estimate the derivatives of the function F . This function F is the fitness function that we want to maximize and that evaluates the discriminative power of a given pair Wl . This fitness function is equal to the statistics obtained with a statistical test that compares the current density in Ωl and Φl between the different classes. More precisely, for each training record m(t), this record is first band-pass filtered in the Φl frequency band by using the corresponding IIR or FIR filter. Then, the current density in Ωl is computed using Eq. 10, and averaged over a given time window and log-transformed as in Eq. 6.

One should note that the obtained values < γ Ωl , Φl > depend on σ v and σ f which are used to compute the band-pass filter and the ROI current density. These < γ Ωl , Φl > are then arranged by class label. The statistical test finally compares the mean value of these < γ Ωl , Φl > between the different classes, the null hypothesis H 0 being "the mean value of the < γ Ωl , Φl > is not different between the classes". The obtained statistics is used as the value of the fitness function F . Thus, algorithm 1 selects the values of σ f and σ v that maximize the discriminative power of a given Wl . Naturally, this procedure is performed for each pair W l . A gradient ascent optimization seemed appropriate here as we experimentally observed that the fitness function F was not monotonic. Rather, this fitness function generally had an optimum for small values of the σ parameter. b) Implementation: Concerning the statistical test used to compute the fitness function F , we used a t-test for binary BCI and an ANOVA for multiclass BCI. These tests are the same ones as those used during the first training step (see III-C). Concerning the values of ǫ f , ǫ v , λ f and λ v , we performed extensive experimental tests and found that ǫ f = ǫ v = 0.0001 and λ f = λ v = 10 -5 were appropriate values. We used these values in all our experiments.

At the end of this offline training, a set of fuzzy ROI Ωl associated to fuzzy frequency bands Φl has been identified. They can now be used for feature extraction, possibly online.

F. Feature Extraction with FuRIA

Once the training is achieved, feature extraction with FuRIA consists in computing the current density in each fuzzy ROI and fuzzy frequency band and in using these current density values as features. More formally, it consists in filtering the EEG signals m(t), once for each one of the N w fuzzy ROI Ωl obtained, using the FIR or IIR filter corresponding to Φl .

Then, < γ Ωl , Φl >, the current density in Ωl , is computed using Eq. 10 and averaged over a given time window and logtransformed as in Eq. 6. The N w current densities < γ Ωl , Φl > are then concatenated into a N w dimensional feature vector [< γ Ω1 , Φ1 >, < γ Ω2 , Φ2 >, . . . , < γ ΩNw , ΦNw >]. Such a feature vector can then be used as an input of any classifier, e.g., a Support Vector Machine, this classifier being in charge of estimating the class of m(t) from the features. Figure 2 summarizes the principle of feature extraction using FuRIA.

G. Model selection

In our implementation of FuRIA, two hyperparameters should be defined by the user: 1) α, the significance threshold used in the statistical test of the first training step (see Section III-C) and 2) H, the smoothing parameters used in the Mean Shift clustering during the second training step (see Section III-D). The values of these hyperparameters have an impact on the number of features N w and on the extension and shape of the ROI and frequency bands. As this impacts the performance of the recognition algorithm, we must select the most appropriate hyperparameters. A solution could be to test several values for H and α and select the couple that enables the best classification on a training set, estimated using a given classifier and Cross Validation (CV). However, we noticed that this method favors models with numerous features, which is not desirable. Furthermore, we observed experimentally that there were generally models with a classification accuracy only slightly lower than the best one, but with much less features. Indeed, if we plot the number of features versus the classification accuracy for several models, the resulting curve tends to be relatively flat for large numbers of features, and suddendly decreases for a smaller number of features (see Figure 3). Ideally, we would like to use the model corresponding to the point of the curve situated just before this sudden decrease of classification accuracy (point A, in blue, on Figure 3), as it would be the best tradeoff between a high classification accuracy and a small number of features. Indeed, we believe that such a model with few features should be prefered as it is probably more robust, less computationally demanding, more easily interpretable and it should ease the training of the classifier. Consequently, we propose a simple model selection criterion C p :

C p = 2 • A -N w ( 15
)
where A is the accuracy (in percent) obtained using CV on a training set. The model with the highest C p is the one that should be prefered. Thus, this criterion is still based on CV but it penalizes models with many features. We also consider that models with a number of features lower than the number of mental states should be avoided. Actually, we consider that a mental state is generated by at least one brain region. In pratice, we observed that this proposed criterion gave satisfactory results. This is shown in the next section which is devoted to the evaluation of FuRIA.

IV. EVALUATIONS OF FURIA

These evaluations have two objectives. First we want to assess the impact of the different hyper-parameters as well as the contribution of the fuzzification processes on the performances. Then we also want to globally assess the efficiency of FuRIA, by comparing a BCI based on FuRIA with other state-of-the-art BCI systems used during BCI competitions.

In order to assess FuRIA, we evaluated it on four different subjects, available from two data sets of BCI competitions 2003 and 2005. For all these evaluations, we worked with a standard head model, composed of three concentric spheres and containing 2394 voxels. This head model has been obtained using LORETA-Key, a software dedicated to inverse solutions (see [START_REF] Pascual-Marqui | Loreta/sloreta website[END_REF] for more details on this software). It could have been better to use subject-specific head models, however, the information needed to build such models was not available in the BCI competition data. For the training of FuRIA, we considered the frequencies located in the 3-45 Hz frequency band with a step of 1 Hz between two consecutive frequencies. Concerning the MCR tests, we used a value of 1000 for the number of random permutations. The FuRIA features were classified using a Gaussian Support Vector Machine (SVM), as SVM is one of the most popular and efficient classifiers used for BCI [START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces[END_REF]. When dealing with multiclass problems, several SVM were combined using the One-Versus-the-Rest (OVR) scheme in order to design a multiclass classifier. The optimal SVM hyperparameters were selected using 10*10 fold stratified cross validation. The descriptions of the data sets and the results obtained are presented in the following sections.

A. EEG data sets 1) BCI competition 2003 -data set IV:

The first data set used was the EEG data set IV of the BCI competition 2003 [START_REF] Blankertz | The bci competition 2003: Progress and perspectives in detection and discrimination of eeg single trials[END_REF], provided by the Berlin group. These data contain EEG signals recorded while a subject was performing self-paced left and right finger tapping tasks. EEG signals were sampled at 100 Hz, recorded using 28 electrodes and comprised the 500 ms before the actual movement. 314 trials were available for training and 100 for testing. The goal of the competitors was to forecast, for each trial, the hand that was used. For this data set, we used FuRIA to learn and extract features on the last 250 ms time window of the data, i.e., we used t 0 = 25 and N s = 25 in Eq. 6. According to several studies, this time window should be the most informative [START_REF] Wang | Bci competition 2003-data set iv: an algorithm based on cssd and fda for classifying single-trial eeg[END_REF] [START_REF] Congedo | Classsification of movement intention by spatially filtered electromagnetic inverse solutions[END_REF]. Here, we used FIR and IIR filters with 24 points and an order 8 respectively.

2) BCI competition 2005 -data set IIIa:

The second data set used was the EEG data set IIIa of the BCI competition 2005 [START_REF] Blankertz | The bci competition iii: Validating alternative approaches to actual bci problems[END_REF], provided by the Graz group [START_REF] Schlögl | Characterization of four-class motor imagery eeg data for the bci-competition 2005[END_REF]. These data were recorded while three subjects S1, S2 and S3, were performing a 4-class motor imagery task. They were instructed to imagine left hand, right hand, foot or tongue movements. For both training and testing, 60 trials were available per class. Trials were sampled a 256 Hz and were recorded using 60 electrodes. Each trial lasted 7 seconds without taking into account the inter trial periods of random lengths. The subjects were instructed to perform the motor imagery tasks during the last 3 seconds of each trial. For all subjects, we subsampled the data at 128 Hz, and used as time window for FuRIA these last 3 seconds, i.e., we used t 0 = 512 and N s = 384 in Eq. 6. For training FuRIA, we ignored all trials contamined by artifacts, leaving approximatly between 25 and 45 trials per class, depending on the subject. Here, we used FIR and IIR filters with 50 points and an order 10 respectively.

B. Evaluation of the influence of hyperparameters and fuzzification processes

In this section, we evaluated FuRIA for different values of the hyperparameter α (see Section III-C), among {0.01, 0.05, 0.1, 0.25, 0.5}, and different values of the hyperparameter H (see Section III-D), among {0.75, 1, 1.25, 1.5, 1.75, 2}. The goal was to assess the impact of these hyperparameters on the classification performances and on the number of generated features. For all these data sets we also compared the results obtained when using FuRIA without the fuzzification process (i.e., using only crisp ROI and frequency bands), with only the spatial fuzzification, with only the frequential fuzzification and with the full (spatial and frequential) fuzzification. In the following, these four conditions are denoted as "Raw", "Freq", "Space" and "All" respectively. We also computed the results for both FIR and IIR filters. However, due to space limitations, only the results for the filters that gave the best results are presented.

1) BCI competition 2003 -data set IV: Tables I and II display the mean classification accuracies obtained on the test set when using FIR filters, for different values of α and H respectively. The mean number of features, denoted as "#Features", is also displayed. In each table, the best result for each condition and each subject is displayed in bold figures. The differences between the Raw condition and the fuzzy conditions, as revealed by a paired t-test, showed that all the fuzzy conditions performed better than the Raw condition on average. However, this difference is only significant for the "All" (p < 0.05) and "Space" (p << 0.01) conditions. 2) BCI competition 2005 -data set IIIa: Tables III and IV display the mean classification accuracies obtained for each of the three subjects on the test set, for different values of α and H respectively. Please note that, for subject S2, no results are presented for α = 0.01 as no pair w k was found significant with this threshold. Here again, only the results obtained with the best filter are displayed, that is IIR filters for S1 and S2 and FIR filters for S3. 3) Discussion: We performed a paired t-test to investigate the overall statistical differences between all the conditions, across all subjects, filter kinds and hyperparameter values. This revealed that globally, all fuzzy conditions performed better than the "Raw" one, and that this difference was significant (p < 0.001). This suggests that for a given set of ROI and frequency bands, fuzzifying them is likely to increase their classification performances. Both the "Space" and "All" conditions performed significantly better than the "Freq" one (p < 0.05). There was no significant difference between these two conditions (p > 0.05), even though the "Space" condition gave slightly better results than the "All" one, on average. We then suggest to use, by default, the spatial fuzzification only, or to try the different fuzzifications and select the most appropriate one for a given subject. The superiority of the spatial fuzzification over the frequential fuzzification could be due to the spatial resolution (thousands of voxels) being higher than the frequential resolution (tens of frequencies). As such, the spatial fuzzification could be more deeply exploited.

Concerning the effects of the threshold α, it can be noticed that the "Raw" condition reached its best performances using the traditional values α = 0.01 or α = 0.05 whereas the fuzzy conditions generally reached their best results for the values α = 0.05 or α = 0.1. This suggests that the fuzzification process enables to use efficiently the information contained in less discriminant voxels and frequencies in order to improve the performances. Regarding the results for the hyperparameter H, the best accuracy is almost always obtained for the lowest value of H. However, accuracies very close to these ones can be obtained with higher value of H and hence, much less features. This supports the use of the previously proposed model selection criterion which penalizes models with a large number of features (see section III-G).

C. Comparison with BCI competition results

In this section, we assess the global efficiency of the FuRIA features by comparing a BCI based on FuRIA features with BCI used by the BCI competition participants. To perform this comparison, we needed to select some parameters without the knowledge of the test set. Consequently, we only relied on cross validation scores on the available training sets. Thus, we selected the kind of filter (FIR or IIR) and the kind of fuzzification (spatial, frequential or both) according to their average cross validation score on the training set. To select the optimal hyperparameters H and α, we relied on the model selection criterion proposed in section III-G, Equation 15.

1) BCI competition 2003 -data set IV: Based only on the training set, the parameter selection procedure described above found that H = 1.75, α = 0.05, FIR filters and only the spatial fuzzification was the most appropriate configuration. This resulted in only 2 features. These parameters led to an accuracy of 84% on the test set, that is, exactly the same score as the winner of the competition [START_REF] Wang | Bci competition 2003-data set iv: an algorithm based on cssd and fda for classifying single-trial eeg[END_REF]. This suggests that the method is efficient, especially when considering the fact that only two features were used. Fig. 4 displays the two ROI and frequency bands learnt by FuRIA. Interestingly, these ROI lie in the left and right motor areas, and the frequency bands lie in the β (≃ 13-30 Hz) band, which is consistent with the literature on the subject [START_REF] Wang | Bci competition 2003-data set iv: an algorithm based on cssd and fda for classifying single-trial eeg[END_REF] [9] [START_REF] Pfurtscheller | Event-related eeg/meg synchronization and desynchronization: basic principles[END_REF]. This suggests that the FuRIA features are interpretable features, which can be used to check what has been learnt or even to extract knowledge about the brain dynamics.

2) BCI competition 2005 -data set IIIa: Table V sums up the parameters used for each subject and selected using only the training set. The resulting number of features is also displayed. Interestingly, different parameters are needed for each subject. This highlights the high subject-specificity of EEG signals and the need to use trainable feature extraction methods. Concerning this data set, the goal of the participants was to provide a continuous classification, i.e., a class label for each time point, here from second 4 to second 7 of each trial. It was indeed during this period that the motor imagery tasks were performed. Hence, the FuRIA features were first learnt on the time window from second 4 to second 7 of each 14-31 Hz 14-28 Hz Fig. 4. The fuzzy ROI (in red) and their corresponding frequency bands that were automatically obtained by using FuRIA for data set IV of BCI competition 2003. The brighter the red color of the voxel, the higher the voxel degree of membership µ Sl (v j ). The brain is seen from the top, nose up. These pictures were obtained with the LORETA-Key software [START_REF] Pascual-Marqui | Loreta/sloreta website[END_REF]. trial. We classified each point by using the FuRIA features extracted from the last 1 second window. Classification outputs were also aggregated across time, which is known to improve accuracy [START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces[END_REF]. In order to do so we used a different multiclass Gaussian SVM (made of several binary SVM combined using the OVR scheme) for each time point of the trial period from second 4 to second 7. These SVM were trained on the FuRIA features extracted on the 1 second window preceding their corresponding time point. Then, to classify each trial, the final output O f (t) at time t was estimated by using the individual outputs O i (t) of the SVM corresponding to the previous time points situated from second 4 to the present time point:

O f (t) = t k=0 w k O i (k) (16)
where the outputs O i are vectors containing the output of each SVM used in the OVR scheme, i.e., one output per class. Here, t = 0 corresponds to second 4 of the trial. We defined w k = A k -25 with A k being the 10*10 fold CV accuracy (in percent) of the k th SVM on the training set. Hence, this method is a weighted combination of classification outputs across time. These weights were chosen so as to ignore contributions from randomly performing classifiers (with a CV error of 25 %, as there are 4 classes) and emphasize contributions of well performing classifiers. The final class C f attributed to a given point t was the one for which

C f = argmax(O f (t)).
The performance measure used was the maximal classification accuracy. Performances obtained using our methods as well as performances obtained by BCI competition 2005 participants on data set IIIa [START_REF] Schlögl | Results of the bci-competition 2005 for datasets iiia and iiib[END_REF] are reported for comparisons in table VI. All these participants also used SVM as classifiers (combined with other classifiers for the 2nd participant). For feature extraction, they all used a frequential information (band-pass filters or amplitude spectra) combined with spatial filters such as common spatial patterns, independant component analysis, and/or surface Laplacian [START_REF] Blankertz | The bci competition iii: Validating alternative approaches to actual bci problems[END_REF]. As shown by table VI, our method outperformed the one of the overall winner of the competition on 2 subjects out of 3 (subjects 1 and 3) and even reached the best score among all participants on subject 3. However, the overall winner reached a really impressive score on subject 2 (the "worst" subject according to the general performances), leading him to the best overall results. Globally, our method reached the second position, with a score only slightly lower than the one of the winner. This shows the efficiency of FuRIA, especially when considering the few features used. Fig. 5 displays the fuzzy ROI and corresponding fuzzy frequency bands automatically learnt by FuRIA for subject 3. Interestingly, the fuzzy ROI identified as relevant are located in the left and right motor areas, for frequency bands located within the µ (9-12 Hz) and β (13-30 Hz) rhythms. This is consistent with the literature on motor imagery [START_REF] Pfurtscheller | Event-related eeg/meg synchronization and desynchronization: basic principles[END_REF] which, again, enhances the interpretability of the extracted features. Fig. 5. The fuzzy ROI (in red) and their corresponding fuzzy frequency band (which is equivalent to the filter magnitude response) that were automatically obtained by using FuRIA, on subject 3 from data set IIIa of BCI competition 2005. The brighter the voxel red color, the higher the voxel degree of membership µ Sl (v j ). The brain is seen from the top, nose up. These pictures were obtained with the LORETA-Key software [START_REF] Pascual-Marqui | Loreta/sloreta website[END_REF].

V. CONCLUSION This paper has presented FuRIA, a trainable feature extraction algorithm for Brain-Computer Interfaces which is based on inverse solutions. This algorithm can be trained to automatically identify relevant regions of interest and their associated frequency bands for the discrimination of mental states, in binary as well as in multiclass BCI. This paper also introduced the concepts of fuzzy ROI and fuzzy frequency bands that enable to use efficiently the available information and, thus, to increase the classification performances.

The evaluation of the proposed method, using sLORETA as the inverse solution and an SVM as classifier, showed its efficiency in terms of classification accuracy. Actually, the obtained results were comparable with those of BCI competition winners. Indeed, it seems that the inverse solution combined with the FuRIA training, acts as a spatial filter that removes the background activity and the noise not correlated with the targeted mental states. As such it focuses on relevant, subjectspecific, brain activity features. An additionnal advantage of FuRIA is the interpretability of the learnt and extracted features. The main drawback of FuRIA is its long training process. Indeed, the discriminative power of several voxels and frequencies should be investigated, which can be long if the number of voxels, electrodes or training data is large, and if time consuming statistical tests such as the ones based on permutation tests are used. However, as this training is performed offline, this point not seem critical.

In future works, it would be interesting to study the influence of both the spatial resolution (number of electrodes used, number of voxels in the chosen head model) and the frequential resolution (number of frequencies investigated) on the performances, in order to possibly reduce the training times and/or improve the quality of the feature extraction. It would also be interesting to take into account the temporal information in FuRIA. To this end, an attractive possibility would be to replace the classical band-pass filters by wavelets, which have proven to be particularly adpated for neuroelectric signals and especially EEG [28]. Finally, it could be interesting to use the dipole orientation in inverse solution-based BCI as current methods only consider the dipole amplitude.

Fig. 2 .

 2 Fig. 2. Feature extraction with FuRIA.

Fig. 3 .

 3 Fig. 3. Example of a plot of the number of features versus the 10*10 fold cross validation accuracy of several FuRIA models. The point A corresponds to the model we would like to automatically select. These data come from subject S1 from the BCI competition 2005 (see Section IV-A2).

TABLE I DATA

 I SET IV, BCI COMPETITION 2003, TEST SET: CLASSIFICATION ACCURACY (%) AND NUMBER OF FEATURES FOR DIFFERENT VALUES OF α, AVERAGED OVER THE DIFFERENT VALUES OF H .

	α	#Features	Raw	Freq	Space	All
	0.01	10.17	80.5	81	82	80.33
	0.05	10.17	79.5	80.67	82.67	82.5
	0.1	9.17	77	78.17	84.67	81.83
	0.25	9.67	76.67 77.17	83.17	76.17
	0.5	8.83	72.17 74.83	78	77.33
			TABLE II
	DATA SET IV, BCI COMPETITION 2003, TEST SET: CLASSIFICATION
	ACCURACY (%) AND NUMBER OF FEATURES FOR DIFFERENT VALUES OF
	H , AVERAGED OVER THE DIFFERENT VALUES OF α.
	H	#Features	Raw Freq	Space	All
	0.75	31.8	83.6	82	83.4 82.6
	1	12.6	80.2	79.6	82.6 80.8
	1.25	6.4	78.2	79	82.6 79.2
	1.5	2.8	74.4	76.8	80 78.4
	1.75	2	73.2	76.6	82 78.2
	2	2	73.4	76.2	82 78.6

TABLE III DATA

 III SET IIIA, BCI COMPETITION 2005, TEST SET: CLASSIFICATION ACCURACY (%) AND NUMBER OF FEATURES FOR DIFFERENT VALUES OF α, AVERAGED OVER THE DIFFERENT VALUES OF H .

	Subject	α	#Features	Raw	Freq	Space	All
		0.01	9.17	75.28 76.02	74.26	73.15
		0.05	9.33	77.96 79.63	76.76	78.46
	S1	0.1	10.67	75.93 77.22	79.17	77.78
		0.25	7.83	60.93 67.96	72.31	68.98
		0.5	9	62.13 66.76	75.19	71.39
		0.05	7.83	54.58 44.31	58.19	55.278
		0.1	11.17	53.33 50.97	57.78	56.11
	S2	0.25	11.67	50.97 43.89	47.22	45.69
		0.5	14.67	50.83 47.92	50.69	51.39
		0.01	6.17	71.67 69.72	71.94	68.75
		0.05	9.83	73.06 76.94	70.42	77.78
	S3	0.1	9.17	70.42 72.08	72.92	72.5
		0.25	11.83	66.39	75	67.92	68.75
		0.5	12.5	64.17 71.39	68.47	67.08

TABLE IV DATA

 IV SET IIIA, BCI COMPETITION 2005, TEST SET: CLASSIFICATION ACCURACY (%) AND NUMBER OF FEATURES FOR DIFFERENT VALUES OF H , AVERAGED OVER THE DIFFERENT VALUES OF α.

	Subject	H	#Features	Raw	Freq	Space	All
		0.75	25.4	86.22 88.67	87.44	89.33
		1	12.2	84.22 85.67	86.56	85.89
	S1	1.25	7.2	77.44 79.33	79.11	77.66
		1.5	5	70.89 73.67	76.11	72.66
		1.75	3.2	58.56 62.33	68.11	63.56
		2	2.2	45.33 51.44	55.89	54.56
		0.75	23.6	61.39 49.44	60.83	56.94
		1	13	58.96 51.46	58.33	56.25
	S2	1.25	8.6	57.50 55.28	56.67	55.83
		1.5	5.2	54.44 53.06	55.56	55
		1.75	2.6	55.56 53.89	56.67	55
		2	1.4	56.11 45.83	54.44	50
		0.75	26	81.67 80.17	83.5	81
		1	13.4	78.33	80	80	80.17
	S3	1.25	7.6	76	78.5	78.17	79.5
		1.5	5.6	74.33	78.5	77.83	79.83
		1.75	4	55.5	63.83	51.83	58.33
		2	2.8	49	57.17	50.67	47

TABLE V PARAMETERS

 V USED FOR DATA SET IIIA FROM BCI COMPETITION 2005

		H	α	filter kind fuzzification	feature number
	S1	1	0.5	IIR	fuzzy space	11
	S2	1	0.05	IIR	fuzzy space	12
	S3 1.25	0.1	FIR	all fuzzy	5

TABLE VI MAXIMAL

 VI CLASSIFICATION ACCURACY FOR THE TEST SETS OF DATA SET IIIA FROM BCI COMPETITION 2005 (%)

		S1	S2	S3	Mean
	Winner	86.67 81.67 85.00	84.44
	2 nd	92.78 57.50 78.33	76.20
	3 rd	96.11 55.83 64.17	72.04
	FuRIA	90.56 69.17 88.33	82.68

ACKNOWLEDGEMENTS

This work was supported by the French National Research Agency within the Open-ViBE project and grant ANR05RNTL01601. The authors also thank Dr. Schlögl, Pr. Yanez-Suarez, Dr. Maby, Dr. Congedo, Dr. Gribonval, Ms. Rosendale and anonymous reviewers for their help and support.


Abstract-This paper presents a general method for studying the mechanical properties of a ferrofluid seal by using a three-dimensional analytical approach based on the coulombian model of a magnet. The fundamental Maxwell's equations lead us to define the concept of magnetic energy of the ferrofluid seal by using only the three-dimensional equations of the magnetic field created by ring permanent magnets radially magnetized.

Our study corresponds to the specific case when the ferrofluid is submitted to a very high magnetic field. Under these conditions, the mechanical properties of the ferrofluid depend only on the magnetic field created by the permanent magnets. Throughout this paper, the remanence polarization J of the magnets used is higher than 1T. The magnetic field we use in order to align the magnetic particles is very intense, greater than 400 kAm -1 . Consequently, the magnetic particles are assumed to be saturated and the magnetic field they create can be omitted.

In this paper, a cylindrical structure consisting of two outer ring permanent magnets and an inner non-magnetic cylinder is considered. In addition, a ferrofluid seal is placed between them. The calculation of the magnetic pressure of the ferrofluid seal has been analytically established in three dimensions in order to determine its shape. Moreover, the geometrical evolution of the ferrofluid seal shape is presented when the inner non-magnetic cylinder crushes the ferrofluid seal. The radial stiffness of the ferrofluid seal is determined in three dimensions when the inner cylinder is decentered. Furthermore, a way of obtaining the ferrofluid seal static capacity is discussed.

It must be emphasized here that the three-dimensional method presented in this paper differs from the ones commonly used in the literature because the only term we consider in the equation of equilibrium describing the state of the ferrofluid seal is determined accurately. Such a method has been used for the design of ironless loudspeakers that have demonstrated some promising acoustical tests.

INTRODUCTION

Analytical calculations of the magnetic field produced by ring permanent magnets are useful for the design of devices using permanent magnets and ferrofluid seals. In this paper, we show how these analytical calculations can be used for studying and optimizing ferrofluid seals in ironless structures. To our knowledge, this way of studying a ferrofluid seal appears for the first time in the literature and we think that it can be useful for the design of magnetic bearings using ferrofluid seals or ironless loudspeakers. Indeed, some recent papers dealing with ironless loudspeakers have clearly taken the same assumptions as the ones used in this paper for studying audio speakers containing ferrofluid [START_REF] Rosensweig | Study of audio speakers containing ferrofluid[END_REF]. The first studies dealing with ferrofluids were certainly carried out by Rosensweig [START_REF] Rosensweig | Ferrohydrodynamics[END_REF] and many authors have been interested in the promising applications using ferrofluids [START_REF] Raj | Advances in ferrofluid in ferrofluid technology[END_REF].

On the other hand, some progress in the analytical calculations of the magnetic field produced by ring permanent magnets [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF]- [START_REF] Selvaggi | Calculating the external magnetic field from permanent magnets in permanent-magnet motors -an alternative method[END_REF] have allowed us to determine accurately the magnetic energy of the ferrofluid by using the fundamental's Maxwell equations.

We can say that studies on ferrofluids deal generally with either their chemical and physical properties or their mechanical behavior [START_REF] Choi | Simulation of hydrostatical equilibrium of ferrofluid subject to magneto-static field[END_REF]. Ferrofluids have various engineering applications [START_REF] Shah | Ferrofluid squeeze film in a long bearing[END_REF], such as damping [START_REF] Bajkowski | A model for a magnetorheological damper[END_REF], for medical applications [START_REF] Park | New design of the magnetic fluid linear pump to reduce the discontinuities of the pumping forces[END_REF] or for ironless loudspeakers [START_REF] Lemarquand | Ironless loudspeakers[END_REF]- [START_REF] Berkouk | Analytical calculation of ironless loudspeaker motors[END_REF]. These last applications are the ones we are currently studying.

Moreover, ferrofluids are also often used as squeeze films [START_REF] Tarapov | Movement of a magnetizable fluid in lubricating layer of a cylindrical bearing[END_REF] in bearings or as rotating shaft seals. Since then, numerous studies in the field of ferrofluid dynamic bearings have been carried out. Both static and dynamic characteristics of these bearings have been studied theoretically [START_REF] Walker | Ferrohydrodynamics thrust bearings[END_REF]- [START_REF] Cunha | A new boundary integral formulation to describe three-dimensional motions of interfaces between magnetic fluids[END_REF]. Moreover, recent trends in the ferrofluid lubrication applications are described and discussed, taking into account various phenomena including cavitation - [START_REF] Meng | An analysis on the magnetic fluid seal capacity[END_REF]. This paper presents a general method for studying the mechanical properties of a ferrofluid seal by using a three-dimensional analytical approach based on the coulombian model of a magnet. The structure considered is composed of two parts that are either made out of permanent magnets (the stator) or non magnetic materials (the moving part). Consequently, no reluctance effect is observed. The moving part has a translation movement, but this study does not take its effects into account. Thus, we only deal with the static properties of the ferrofluid seal. The magnetic field created by the ring permanent magnets is very intense. Therefore, the ferrofluid is largely saturated. Although numerical approaches are commonly used to study ferrofluid seals, analytical approaches with exact formulations allow a great flexibility in optimizing the dimensions of the configurations studied. Consequently, in our specific use of the ferrofluid seal, it is more interesting to use an exact analytical term rather than a lot of simplified terms in the Bernoulli equation for describing the state of the ferrolfuid seal.

It is noted that this study is entirely based on the magnetic pressure calculation which is analytical throughout this paper. Indeed, such an analytical approach allows a precise study of the ferrofluid seal shape when it is crushed, which is often an occurring state during the functioning. The seal capacity is also determined and discussed and a analytical method for obtaining the radial stiffness of the ferrofluid seal is also presented.

We think that this method paves the way for new and promosing applications using both permanent magnets and ferrofluid seals. All the analytical expressions used in this paper have been published in previous papers and are a good alternative to the finite element method for the design of such structures.

STRUCTURE DESCRIPTION

Geometry

The studied structure consists of two outer ring permanent magnets and an inner non-magnetic cylinder. Such a structure is considered in this paper because it allows us to have important gradient values. Moreover, this structure is not very complicated to build since it requires only two ring permanent magnets. As a result, the magnets are stacked and their magnetizations are both radial, but in opposed directions (Fig. 1). The non-magnetic cylinder is centered with the two ring permanent magnets. Moreover, both outer ring permanent

h Z 0 R r in r out Figure 1
. Geometry : two outer ring permanent magnets and an inner non-magnetic cylinder with a ferrofluid seal between them; the ring inner radius is r in , the ring outer radius is r out , the height of a ring permanent magnet is h. magnets have the same dimensions and the same magnetization in magnitude. The ring inner radius is r in , the ring outer radius is r out and the ring permanent magnet height is h. The z axis is a symmetry axis.

Basic equations

The magnetic field produced by the ring permanent magnets can be determined by using a fully analytical approach. The fundamental equations of Maxwell lead us to determine the magnetic energy in the ferrofluid. Indeed, let us consider the four fundamental Maxwell's equations:

∇. B = 0 (1)

∇ ∧ H = j (2)
∇. D = ρ (3)
∇ ∧ E = - ∂ B ∂t ( 4
)
where B is the magnetic induction field, H is the magnetic field, j is the volume current density, D is the electric flux density, E is the electrostatic field and ρ is the electrical charge. In our structure, no currents are considered as the magnetic field is created only by the permanent magnets. The vector fields B and H are defined for all points in space with the following relation:

B = µ 0 H + J (5)
where µ 0 is the permeability of the vacuum and J is the polarization vector of the magnet. When the magnetic field is determined outside of the magnet, we have J = 0. The analogy between the Mawxell's equations leads us to write that:

∇. H = - ∇. J µ 0 = σ * µ 0 (6)
where σ * corresponds to a fictitious magnetic pole density. On the other hand, the magnetic field H verifies the following equation:

∇ ∧ H = 0 (7)
Thus, H can de deducted from a scalar potential φ( r) by

H = -∇ (φ( r)) (8)
Let us consider a structure owing several ring permanent magnets, that is, an assembly of magnetic pole surface densities and magnetic pole volume densities; with ( 6) and ( 7), we can write that

φ( r) = 1 4πµ 0   i S i J k .d S i r -r ′ + j V j -∇. J k r -r ′ dV j   (9)
where J k is the magnetic polarization vector owing to the ring permanent magnet k and rr ′ is the distance between the observation point and a magnetic charge contribution. We deduct immedialtely that the magnetic field created by the ring permanent magnets is determined as follows: As there is a cylindrical symmetry in our configuration, the magnetic field depends only on two variables r and z. This relation is used in the next section for showing the magnetic field in our configuration.

H = -∇.   1 4πµ 0   i S i J k .d S i r -r ′ + j V j -∇. J k r -r ′ dV j     (10)

Magnetic field created by the two ring permanent magnets

As shown in Fig 1, our geometry is comoosed of two ring permanent magnets. The magnetic field created by the two outer ring permanent magnets depends on both the radius, r, and the altitude, z. The radial component of the magnetic field created by the two outer ring permanent magnets is denoted H r (r, z) and the axial component is denoted H z (r, z). It is noted that the expressions of these components have been established in previous papers [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF]- [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF]. The magnetic field H is represented in Fig. 2. The ring inner radius r in equals 25 mm, the ring outer radius equals 28 mm and the permanent magnet height, h, equals 3 mm. The magnetic pole surface density +σ * equals 1 T. Figure 2 shows that the magnetic field H is the strongest near the permanent magnets, especially where the magnetic field gradient is the most important.

MAGNETIC PRESSURE ON THE FERROFLUID SEAL

Assumptions

It is emphasized that the structure considered in this paper contain no iron. Moreover, the further calculations are presented for magnets with 1 T remanent magnetization in order to normalize the results.

In fact, the magnets used in the prototypes are Neodymium Iron Boron ones for which the remanent magnetization can reach 1.5 T. Furthermore, the interesting regions of space are the ones where the ferrofluid goes. Besides, in the field gradient due to the external magnetic field created by the two outer ring permanent magnets, the magnetic particles are attracted towards the higher intensity regions of the magnetic field. With these 1 T normalized magnets and for the proposed configurations, the magnetic field intensity there is greater than 400 kA/m. Of course, all the field values are proportional to remanent magnetization value. Then, we use commercial ferrofluids either from the company Ferrotec or Ferrolabs. Such ferrofluids have a saturation magnetization, M s , smaller than 32 kA/m and a particle concentration below 5.5 %. It is to be noted that for bearing or loudspeaker applications, a great bearing effect is seeked which requires high saturation magnetizations. Therefore, the magnetic field, H, created by the permanent magnets is far higher than the ferrofluid critical field [START_REF] Matthies | Numerical simulation of normalfield instability in the static and dynamic case[END_REF]. So, the ferrofluid is totally saturated and its magnetization is denoted M s . Another consequence is that the magnetic relative permeability of the ferrofluid equals one. Thereby, the presence of ferrofluid does not modify the magnetic field created by the permanent magnets and the field created by the ferrofluid itself is omitted. Moreover, the ferromagnetic particles are assumed to be small saturated spheres which can be freely oriented in all the directions of space. Besides, this study is a static one. Thus, all the particles of the saturated ferrofluid are aligned with the permanent magnet orienting field. Consequently, the ferrofluid magnetization has the same direction as the orienting field. Furthermore, the aggregation in chains of the ferrofluid particles is omitted [START_REF] Ivanov | Ferrofluid aggregation in chains under the influence of a magnetic field[END_REF]. It is noted that, when the device is at laboratory temperature and at rest, the aggregation phenomenon is observed.

Moreover, the thermal energy E T (E T = kT where k is Boltzmann's constant and T is the absolute temperature in degrees Kelvin) and the gravitational energy E G (E G = ∆ρV gL where V is the volume for a spherical particle, L is the elevation in the gravitational field, g is the standard gravity, ∆ρ is the difference between the ferrofluid density and the outer fluid) are neglected.

Eventually, this paper deals with the ferrofluid free boundary surface. Its shape depends on the result of the force competition at this boundary. On one hand, the magnetic pressure, p m , is exerted. It is given by:

p m (r, z) = µ 0 M s . H(r, z) = µ 0 M s H r (r, z) 2 + H z (r, z) 2 (11)
where both magnetic field components H r (r, z) and H z (r, z) are analytically calculated [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF]. It is to be noted that this expression is given below in a fundamental form. However, this expression can be written in a form more suitable for numerical treatments [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF]. The magnetic pressure p m (r, z) is given by [START_REF] Selvaggi | Calculating the external magnetic field from permanent magnets in permanent-magnet motors -an alternative method[END_REF].

p m (r, z) = µ 0 M s H r (r, z) 2 + H z (r, z) 2 (12)
where M s is the intensity of magnetization of a magnetic particle. The radial component H r (r, z) of the magnetic field created by the two permanent magnets is defined by [START_REF] Choi | Simulation of hydrostatical equilibrium of ferrofluid subject to magneto-static field[END_REF].

H r (r, z) = σ * πµ 0 i(1 + u) (ς(u 1 ) -ς(u 2 )) (13)
with

ς(u) = ξ 1 (-(a 1 d + b 1 (c + e 1 )))F * i sinh -1 [ √ -c+d-e 1 √ c+e 1 +du ], c+d+e 1 c-d+e 1 d √ -c + d -e 1 e 1 d(1+u) c+e 1 +du √ 1 -u 2 + ξ 1 (b 1 c -a 1 d)Π * e 1 c-d+e 1 , i sinh -1 [ √ -c+d+e 1 c+e 1 +du ], c+d+e 1 c-d+e 1 d √ -c + d -e 1 e 1 d(1+u) c+e 1 +du √ 1 -u 2 + ξ 2 (-(a 2 d + b 2 (c + e 2 ))) F * i sinh -1 [ √ -c+d-e 2 √ c+e 2 +du ], c+d+e 2 c-d+e 2 d √ -c + d -e 2 e 2 d(1+u) c+e 2 +du √ 1 -u 2 + ξ 2 (b 2 c -a 2 d)Π * e 2 c-d+e 2 , i sinh -1 [ √ -c+d+e 2 c+e 2 +du ], c+d+e 2 c-d+e 2 d √ -c + d -e 2 e 2 d(1+u) c+e 2 +du √ 1 -u 2 - η 3 ((a 3 d -b 3 e 3 )) F * i sinh -1 [ √ -d-e 3 √ e 3 +du ], -d-e 3 d+e 3 d √ -d -e 3 (-c + e 3 ) d(1+u) e 3 +du √ 1 -u 2 - η 3 (b 3 c -a 3 d)Π * -c+e 3 d+e 3 , i sinh -1 [ √ -d+e 3 e 3 +du ], -d+e 3 d+e 3 d √ -d -e 3 (-c + e 3 ) d(1+u) e 3 +du √ 1 -u 2 - η 4 (a 4 d -b 4 e 4 ) F * i sinh -1 [ √ -d-
The axial component of the magnetic field created by the two ring permanent magnets is given by [START_REF] Bajkowski | A model for a magnetorheological damper[END_REF].

H z (r, z) = σ * πµ 0   -r in K * -4rr in (r-r in ) 2 +z 2 (r -r in ) 2 + z 2   + σ * πµ 0   r in K * -4rr in (r-r in ) 2 +(z-h) 2 (r -r in ) 2 + (z -h) 2   - σ * πµ 0   r in K * -4rr in (r-r in ) 2 +z 2 (r -r in ) 2 + z 2   + σ * πµ 0   r in K * -4rr in (r-r in ) 2 +(z+h) 2 (r -r in ) 2 + (z + h) 2   (15)
ξ i = d(-1 + u) c + e i + du (16)
η i = d(-1 + u) e i + du (17)
where K * [m] is given in terms of the incomplete elliptic integral of the first kind by ( 18)

K * [m] = F * [ π 2 , m] (18)
F * [φ, m] is given in terms of the elliptic integral of the first kind by [START_REF] Berkouk | Analytical calculation of ironless loudspeaker motors[END_REF]:

F * [φ, m] = θ=φ θ=0 1 1 -m sin(θ) 2 dθ ( 19
)
Parameters a 1 r in rz b 1 -r 2 in z c r 2 + r 2 in d -2rr in e 1 z 2 a 2 -r in r(z -h) b 2 r 2 in (z -h) e 2 (z -h) 2 a 3 r in rz b 3 -r 2 in z e 3 r 2 + r 2 in + z 2 a 4 r in r(-z -h) b 4 -r 2 in (-z -h) e 4 r 2 + r 2 in + (z + h) 2
Table 1. Definition of the parameters used in ( 14)

Π * [n, φ, m] is given in terms of the incomplete elliptic integral of the third kind by ( 20)

Π * [n, φ, m] = φ 0 1 (1 -n sin(θ) 2 ) 1 -m sin(θ) 2 dθ (20)
The parameters used in ( 14) are defined in Table 1. Moreover, when we input [START_REF] Choi | Simulation of hydrostatical equilibrium of ferrofluid subject to magneto-static field[END_REF] in Mathematica, we have to take the real part of H r (r, z) because of the noise calculus. What's more, the parameter i used in ( 14) is the imaginary number (i 2 = -1).

As the ferrofluid considered in this paper has a saturation magnetization of 32 kA/m and the magnetic field is greater than 400 kA/m, the value of this pressure is higher than 12800 N/m 2 . On the other hand, the surface tension exists. But when the values of both the surface tension coefficient, A, (A equals 0.0256 kg/s 2 for the used ferrofluids) and the radius of curvature are considered, the effect of the surface tension can be omitted: this latter does not deform the free boundary surface.

Finally, for hydrodynamic pressures which equal zero or have low values, the seal free boundary surface is a magnetic iso-pressure surface. Then, its mechanical properties, such as its capacity or its stiffness can be studied. Further dynamic studies could also be carried out with a perturbation method from this rest position [START_REF] Kuzhir | Free boundary of lubricant film in ferrofluid journal bearings[END_REF]. The magnetic pressure p m (r, z) has been determined with [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF]. Figure 3 shows that the magnetic pressure is greater near the two ring permanent magnets, especially where the magnetic field gradient is the strongest. It can be noted that such a figure is useful to see quickly where the magnetic pressure is the most important. Moreover, Fig.

(3) gives also an indication about how the magnetic pressure varies according to r and z; we see that this variation is very important either to the permanent magnets. Furthermore, this representation is interesting since it shows that the potential energy is concentrated in a very small ferrofluid volume. As a consequence, it gives an indication about what quantity of ferrofluid should be used to design a ferrofluid seal. Indeed, some viscous effects can appear when the inner cylinder moves between the two outer ring permanent magnets. If a large quantity of ferrofluid is used, the ferrofluid seal is thick, the potential energy increases but the viscous effects become an actual drawback according to the dynamic movement of the inner cylinder. If a too small amount of ferrofluid is used, the viscous effects disappear but all the interesting properties of the ferrofluid seal (damping, stability, linearity,...) disappear as well. In short, to a given geometry ( here two ring permanent magnets with an inner non-magnetic cylinder) corresponds to an adequate quantity of ferrofluid which has some interesting physical properties with very little viscous effects.

Potential energy of the ferrofluid seal

The potential energy of the ferrofluid seal is defined by ( 21):

E m = - (Ω) p m (r, z)dV (21)
where (Ω) is the ferrofluid seal volume. Indeed, this energy allows the calculation of the seal mechanical properties. As a remark, the magnetic pressure is given in N/m 2 and the potential energy in J.

STUDY OF THE FERROFLUID SEAL SHAPE

Ferrofluid seal shape under a magnetic field

We first analyze the ferrofluid shape without taking into account the inner non-magnetic cylinder. As we assume each magnetic particle to be a small sphere which can be freely oriented in all the directions of space, each magnetic particle is oriented with the direction of the orienting field in order to minimize its potential magnetic energy. More exactly, the magnetic particles go where the magnetic field is the most important by always being aligned with it. Figure 4 is a two dimensional contour representation of the potential energy of the ferrofluid seal. Such a figure is useful to see the points in space which have the same potential energy. This is interesting because this element of information allows us to study the deformation of the seal when the inner cylinder is radially decentered. The aim of this figure is to show where the ferrofluid goes at first: the warmest areas are the first occupied ones. The shape of the ferrofluid seal between r = 24.6 mm and r = 25 mm can be approximated as half an ellipse. Consequently, when the ferrofluid seal thickness is smaller than 0.4 mm, its contour can be written in terms of an equation of an ellipse [START_REF] Tiperi | Overall characteristics of bearings lubricated ferrofluids[END_REF].

(r -r i ) 2 a 2 i + z 2 b 2 i = 1 (22)
To illustrate [START_REF] Tiperi | Overall characteristics of bearings lubricated ferrofluids[END_REF], Table [START_REF] Rosensweig | Ferrohydrodynamics[END_REF] gives the parameters of the equations of ellipse describing the ferrofluid seal in Fig. 4 for r between 24.6 mm and 25 mm. The approximation error between the equations of ellipse and the real contour shape of the ferrofluid seal is also given. When some ferrofluid is added, the shape of the ferrofluid seal changes and for Table 2. Parameters (in mm) of the equations of an ellipse describing small ferrofluid seals r < 24.6 mm, as represented in Fig. 5, the ferrofluid seal contour cannot be represented in terms of an equation of ellipse. Furthermore, we can point out that if a large amount of ferrofluid is added, the ferrofluid seal size becomes very large because of the magnet edge effects. In our configuration, the thickness of the ferrofluid seal for which the shape changes is 0.6 mm. This element of information is very important since it gives clearly what quantity of ferrofluid sould be used to make a small seal. We have also represented the final three-dimmensional form of the ferrofluid seal with the same dimensions (r in = 25 mm. J = 1 T, h = 3 mm) in Fig 6. . Representation of growing ferrofluid volumes under a magnetic field from the warmest areas to the coolest ones: the warmest colors correspond to the region in space where the magnetic energy is the strongest

Shape of the crushed ferrofluid seal

We take now the non-magnetic cylinder into account. The device has to be dimensioned in such a way so that the cylinder crushes the seal, in order to reach the intended watertightness. Figure 7 shows where the ferrofluid goes when the inner non-magnetic cylinder crushes the ferrofluid seal. The phenomenon is axisymmetrical, as the inner cylinder is supposed here to be centered. As seen in Fig. 4, for small quantities of ferrofluid, the ferrofluid seal contour can be described in terms of an equation of ellipse. Therefore, when a small ferrofluid seal is crushed, the quantity of ferrofluid goes symmetrically under and above the initial ferrofluid seal by describing an equation of ellipse. The change of the ferrofluid seal shape generates an energy reduction. Table 3 presents the energy reduction when the cylinder, being first at a distance in length which equals 0.3 mm from r in , crushes radially the ferrofluid seal. The energy reduction is very important (68%) when the radial displacement of the inner cylinder is 0.2 mm.

CALCULATION OF THE STATIC CAPACITY OF THE FERROFLUID SEAL

The aim of this section is to determine the ferrofluid seal static capacity. Indeed, the pressure underneath the non-magnetic cylinder can be Consequently, a pressure gradient appears and the ferrofluid seal is thus deformed. Therefore, the calculation of the ferrofluid seal capacity is required. For this purpose, a second configuration is considered, which corresponds to the case when a cylindrical air gap appears in the seal along the cylinder because of an applied pressure on one side of the seal. This second configuration is represented in Fig. 8. Two steps are necessary to determine the capacity of the ferrofluid seal. The first step is the determination of the potential energy of the ferrofluid seal without any air gap in the seal. To do so, a numerical integration of ( 21) leads to a first value of the potential energy which is denoted E m (1). The used numerical method is the GaussKronrod method.

A three-dimensional representation of the ferrofluid seal perfored in shown in Fig 9 . The second step is the determination of the potential energy of the ferrofluid seal with the air gap. Again, a numerical integration leads to a second value of the potential energy denoted E m (2). The energy difference is denoted ∆E m and verifies [START_REF] Raikher | Orientational dynamics in magnetic fluids under strong coupling of external and internal relaxations[END_REF]:

∆E m = E m (1) -E m (2) (23)
The energy difference corresponds to the pressure work δW (P ) and satisfies [START_REF] Miyake | Sliding bearing lubricated with ferromagnetic fluid[END_REF]: ∆E m = δW (P ) = P Sd [START_REF] Miyake | Sliding bearing lubricated with ferromagnetic fluid[END_REF] where S is the surface area of the air gap and d is the thickness of the perforation. Consequently, the capacity P lim verifies [START_REF] Chang | A theoretical and experimental study of ferrofluid lubricated four-pocket journal bearing[END_REF]:

P lim = δW (P ) Sd (25)
A numerical application has been done with several ring inner radii (Fig. 10). We define the ferrofluid thickness as the thickness of the 4. Volume and Magnetic field corresponding to a given ferrofluid thickness ferrofluid between the outer ring permanent magnets and the inner cylinder. Figure 10 shows that the thinner the ferrofluid seal is, the higher the pressure gradient to which the ferrofluid seal can resist is. Table 4 gives both the volume corresponding to each ferrofluid thickness and the smallest magnetic field Hlim in the concerned volume.

RADIAL STIFFNESS OF THE FERROFLUID SEAL

Determination of the radial stiffness

This section presents a way of obtaining the radial stiffness of the ferrofluid seal. For this purpose, a third configuration is considered (Fig. 11), which corresponds to a decentered inner cylinder. Figure 12 illustrates the dependence of the radius with θ. The ferrofluid seal is thus more crushed when the cylinder comes closer to the magnets. The problem here is no longer axisymmetrical, but fully three-dimensional. Indeed, where the radial room for the ferrofluid is narrower, the seal becomes larger along the Oz axis, but the ferrofluid also moves towards the region of increasing room and so, the angle θ interferes. The way of obtaining the radial stiffness is done in two steps. The potential energy E m (1) is first calculated when the nonmagnetic cylinder is centered (Fig. 1) with [START_REF] Walker | Ferrohydrodynamics thrust bearings[END_REF]. Then, the potential energy E m (2) is calculated when the non-magnetic cylinder is not centered (Fig. 11). In the second configuration, the limits of the integrals which allow us to determine the potential energy depend on the angle θ. Consequently, the potential energy E m (2) is determined with [START_REF] Zhang | Static characteristics of magnetized journal bearing lubricated with ferrofluids[END_REF].

E m (2) = 2π 0 r in r 0 +e cos(θ) z 0 -z 0 e m (r, z)rdrdθdz (26)
where r 0 and e are determined by the equation of the circle which is not centered given by [START_REF] Osman | Static and dynamic characteristics of magnetized journal bearings lubricated with ferrofluid[END_REF].

r(θ) = r 0 + e cos(θ) (27)
where r 0 = 24.7 mm and e = 0.1 mm. Figure 12 shows how the inner cylinder is decentered. Consequently, the ferrofluid shape is not axisymmetrical since the inner cylinder is radially decentered.

As the ferrofluid shape is not axisymmetrical, the limits of the integrals in Eq. 26 depend on the angle θ. Thanks to the determination of E m (1) and E m (2), we can calculate the radial force F r with [START_REF] Shah | Anisotropic permeable porous facing abd slip velocity squeeze film in axially undefined journal bearing with ferrofluid lubricant[END_REF].

F r = E m (1) -E m (2) 2∆r (28)
where ∆r = 0.001 m is the radial distance in length for which the inner cylinder is moved away. The radial stiffness k r is determined by [START_REF] Cunha | A new boundary integral formulation to describe three-dimensional motions of interfaces between magnetic fluids[END_REF].

k r = F r ∆r ( 29
)
The numerical application gives k r = 5.6 N/mm. As a conclusion, for a ferrofluid seal whose thickness equals 0.3 mm, the radial stiffness is k r = 5.6 N/mm.

Decentering supported by the ferrofluid seal

The determination of the maximum possible decentering depends on the quantity of ferrofluid used in the seal. A numerical application is given here to illustrate that. The ferrofluid quantity we use corresponds to the case when its thickness equals 0.4 mm. With the inner cylinder, the air gap equals 0.3 mm. It is noted that the inner cylinder is radially moved away. As some ferrofluid is moved off from the place where the ferrofluid seal is crushed, the ferrofluid seal thickness increases at the opposite from the place where the ferrofluid seal is crushed. The moving off centre for which a perforation can appear in the seal corresponds to 0.14 mm. As a conclusion, in this configuration with r in = 25 mm, r out = 28 mm and h = 3 mm, the moving off centre supported by the ferrofluid seal equals 0.14 mm.

CONCLUSION

This paper has presented an analytical way of studying analytically the mechanical properties of a ferrofluid seal. This static study is based on the potential energy criterion. Moreover, the magnetic field created by the ring permanent magnets being very intense, the magnetic field created by the magnetic particles is omitted. The expression of the magnetic pressure has been determined analytically and has been plotted in order to see where the potential energy is the most important. The shape of the ferrofluid seal depends on the quantity of ferrofluid which is used and we have shown that for small quantities, the ferrofluid seal contour can be described in terms of an equation of ellipse. Moreover, if the thickness of the ferrofluid seal is greater than 6 mm, the ferrofluid seal spreads over the whole inner faces of the ring permanent magnets. It can be noted that the shape of the ferrofluid seal changes when the inner cylinder is taken into account. This change comes with an energy reduction which depends on the radial displacement of the cylinder, so on the seal crushing. The ferrofluid seal capacity has been determined and we have shown that it depends on the thickness of the ferrofluid seal. At last, a way of obtaining the radial stiffness of a ferrofluid seal is presented and we have found a radial stiffness which equals 5.6 N/mm. Furthermore, the moving off centre supported by the ferrofluid seal has been studied and we have found that it depends on the quantity of ferrofluid used in the seal. These results can be useful to design many engineering and industrial applications which use ferrofluid seals.

  MECHANICAL PROPERTIES OF A FERROFLUID SEAL: THREE-DIMENSIONAL ANALYTICAL STUDY BASED ON THE COULOMBIAN MODEL R. Ravaud, G. Lemarquand Laboratoire d'Acoustique de l'Universite du Maine, UMR CNRS 6613 Avenue Olivier Messiaen, 72085 Le Mans, France guy.lemarquand@univ-lemans.fr

Figure 2 .

 2 Figure 2. Magnetic field created by the two outer ring permanent magnets; r in = 25 mm, h = 3 mm, σ * = 1 T

e 4 √ e 4 1 -u 2 - η 4 √ -d-e 4 e 4

 412444 +du ], -d+e 4 d+e 4 d √ -d -e 4 (c + e 4 ) d(1+u) e 4 +du √ (b 4 ca 4 d)Π * -c+e 4 d+e 4 , i sinh -1 [ +du ], -d+e 4 d+e 4 d √ -d -e 4 (-c + e 4 ) d(1+u) e 4 +du √ 1 -u 2

3 Figure 3 .Figure 3

 333 Figure 3. Three-dimensional representation of the magnetic pressure of the ferrofluid seal

Figure 4 .

 4 Figure 4. Representation of growing ferrofluid volumes under a magnetic field from the warmest areas to the coolest ones; the warmest colors correspond to the region in space where the magnetic energy is the strongest

Figure 5

 5 Figure 5. Representation of growing ferrofluid volumes under a magnetic field from the warmest areas to the coolest ones: the warmest colors correspond to the region in space where the magnetic energy is the strongest

Figure 6 .Table 3 .

 63 Figure 6. Tree-dimensional representation of the ferrofluid seal volume; the ring inner radius equals 25 mm. J = 1 T, h = 3 mm

Figure 7 .

 7 Figure 7. Ferrofluid seal which is crushed

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. Geometry: two ring permanent magnets and an inner nonmagnetic cylinder; air gap in the seal

Figure 11 .Figure 12 .

 1112 Figure 11. Radially moved away cylinder, crushed ferrofluid