Emeline Perthame

Ching-Fan Sheu

David Causeur May

Signal identification in ERP data by decorrelated Higher Criticism Thresholding

Event-related potentials (ERPs) are intensive recordings of electrical activity along the scalp time-locked to motor, sensory, or cognitive events. A main objective in ERP studies is to select (rare) time points at which (weak) ERP amplitudes (features) are significantly associated with experimental variable of interest. The Higher Criticism Thresholding (HCT), as an optimal signal detection procedure in the "rare-and-weak" paradigm, appears to be ideally suited for identifying ERP features. However, ERPs exhibit complex temporal dependence patterns violating the assumption under which signal identification can be achieved efficiently for HCT. This article first highlights this impact of dependence in terms of instability of signal estimation by HCT. A factor modeling for the covariance in HCT is then introduced to decorrelate test statistics and to restore stability in estimation. The detection boundary under factor-analytic dependence is derived and the phase diagram is correspondingly extended. Using simulations and a real data analysis example, the proposed method is shown to estimate more efficiently the support of signals compared with standard HCT and other HCT approaches based on a shrinkage estimation of the covariance matrix.

Introduction

Event-related potentials (ERP) are voltage changes along the scalp time-locked to some physical or mental occurrence in the ongoing electrical brain activity recorded as electroencephalogram (EEG). Like functional magnetic resonance imaging (fMRI), ERPs are noninvasive instruments that directly reflect cortical neuronal activity. Unlike fMRI, ERPs studies provide better temporal resolution to chart the time course of mental processes and are less expensive to conduct. In basic research, ERPs offer a psychophysiological method for studying attentional processes, language, and memory functions, yielding information not available from behavioral studies alone. In clinical research, ERPs are one of several noninvasive biomarkers that have been proposed for evaluating neurological and psychiatric disorders such as Alzheimer's disease, amnestic mild cognitive impairment, attention deficit, hyperactivity disorder, among others.

For a typical trial in a study, ERP amplitudes are measured in milliseconds (ms) for up to one or more seconds with reference to the onset of an external event. ERP waveforms are notoriously noisy and highly variable, both within and between subjects, which explains why ERPs are usually averaged across trials of the same condition for the same subjects. To identify time points at which ERP amplitudes can reliably be linked to either stimulus (or response) events, researchers must shift, simultaneously, through thousands of features for significance testing. A balance must be struck between keeping a low false positive error rate while maintaining sufficient power for correct signal identification. How to achieve this objective for ERP data exhibiting arbitrarily strong temporal dependence is the focus of the present paper.

Searching for time points at which ERP amplitudes are significantly associated with a response variable can be seen as a signal identification issue in the "Rare and Weak" (RW) paradigm introduced by [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF]. Large-scale significance analysis of ERP waveforms is indeed based on a T -vector T = (T 1 , . . . , T T ) of test statistics for the collection of corresponding null hypotheses H 0,t of no association between the ERP measured at time t and the target variable. The RW paradigm provides a simple yet useful framework to investigate test procedures for detecting the presence of a nonzero signal. In the multiple testing setting, [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF] proposed Higher Criticism Thresholding (HCT) inspired by an idea of [START_REF] Tukey | T13 N: the higher criticism[END_REF] based on the significance of an overall body of tests. HCT is known to be effective to detect signals under independence. Indeed, closed-form detection bounds can be derived analytically and [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF] demonstrated that HCT attains the theoretically optimal decision limits. For the more challenging objective of selecting non-null features for classification or prediction, [START_REF] Donoho | Higher criticism thresholding: Optimal feature selection when useful features are rare and weak[END_REF] also demonstrated the superior performance of HCT with respect to multiple testing procedures controlling False Discovery Rate (FDR).

As reported in [START_REF] Causeur | A factor-adjusted multiple testing procedure for erp data analysis[END_REF], the pronounced dependence observed in ERP waveforms can induce a longrange regularity for the test statistics resulting in spuriously small p-values outside of the support of the signal and causing the non-null features to be misidentified. This instability in the ranking of p-values due to dependence has also been reported in high-dimensional genomic data analysis (see for example [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF], [START_REF] Ahdesmäki | Feature selection in omics prediction problems using cat scores and false non-discovery rate control[END_REF], [START_REF] Perthame | Stability of feature selection in classification issues for high-dimensional correlated data[END_REF]). Although HCT is known to be effective even when tests are weakly correlated (see [START_REF] Hall | Properties of higher criticism under strong dependence[END_REF]), its performance can still be improved when dependence is accounted for ( [START_REF] Ahdesmäki | Feature selection in omics prediction problems using cat scores and false non-discovery rate control[END_REF] and [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF]). For example, [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF] reported that the theoretical detection bounds derived in the RW framework are markedly affected by a strong dependence among the test statistics. Therefore, [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF] extended the RW framework by introducing the so-called innovated HCT (iHCT) and showed that iHCT restores the effectiveness of the HCT procedure under strong dependence.

Elsewhere, for feature selection in the Linear Discriminant Analysis (LDA) context, [START_REF] Ahdesmäki | Feature selection in omics prediction problems using cat scores and false non-discovery rate control[END_REF] applied HCT to Correlated Adjusted T-scores (CAT-scores), which are tests statistics decorrelated by a James-Stein shrinkage estimator of the covariance matrix (see [START_REF] Zuber | Gene ranking and biomarker discovery under correlation[END_REF]). [START_REF] Ahdesmäki | Feature selection in omics prediction problems using cat scores and false non-discovery rate control[END_REF] showed that performance of HCT is improved by this decorrelation using an inverse-square root of the covariance matrix. Alternatively, [START_REF] Perthame | Stability of feature selection in classification issues for high-dimensional correlated data[END_REF] extended the decorrelation method introduced in [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF] based on a regression factor model, to generalized linear modeling for a variable selection issue in supervised classification. [START_REF] Perthame | Stability of feature selection in classification issues for high-dimensional correlated data[END_REF] devised an adaptive algorithm to jointly estimate signal, covariance matrix and the class probabilities in order to efficiently decorrelate data for classification.

As demonstrated by [START_REF] Causeur | A factor-adjusted multiple testing procedure for erp data analysis[END_REF], the complex dependence pattern observed in the correlation structure of test statistics derived from ERP data can be well approximated using a factor decomposition, often with a moderate number of factors. The procedure proposed in this paper is expected to better select ERP features which are associated with treatment variable of interest. Moreover, it also takes advantage of simple and efficient algebraic tools to derive an inverse-square root of the covariance matrix. This article is organized as follows. Section 2 describes the oddball auditory ERP paradigm to motivate the subsequent modeling and data analysis issues. It also introduces signal detection methods in a multivariate linear settings for ERP experiments. A brief review of actual methods dealing with HCT under independence and dependence is presented in Section 3. The factor model framework is introduced in Section 4, leading to an extension of detection boundary under dependence and a Factor-innovated HCT method. The properties of the method are demonstrated through simulations and a real data study in Section 5. We conclude with a discussion in Section 6.

Signal detection in ERP experiments

We first consider hereafter the statistical testing for the mean comparison of I groups of ERP curves, with n i curves in group i. The n = n 1 + . . . + n I ERP curves are observed at timepoints in {t 1 , . . . , t T } on J subjects. The former signal detection issue is first motivated by an auditory oddball ERP study.

The auditory oddball experiment

In ERP studies, perhaps the most commonly used experimental task is the oddball paradigm ( [START_REF] Picton | The p300 wave of the human event-related potential[END_REF]). In this paradigm, typically two classes of stimuli are presented, one occurring frequently (standard) and the other occurring infrequently (target). The subject is required to distinguish between the two stimuli and to respond to the stimuli that are designated as targets.

An auditory ERP study was performed at Kaohshung Medical University in Taiwan to provide an illustrative data set for the present investigation. The task uses two pure tones of 500 Hz and 1,000 Hz. The former is presented 120 out of 150 trials, whereas the latter (target) is presented only for 30 trials. The order of tone presentation is random and the subject is asked to (silently) count the number of targets. At each of 4 electrode locations (FZ, C3, C4, & O1), ERP waveform was obtained from each of the two tone conditions for each of the J = 15 participants, each curve begins at -100 milliseconds (ms) and terminates at 399.5 milliseconds (ms) with two records per 1 ms. The stimulus onset is at 0 ms. For subsequent analysis, only the ERPs from the electrode location FZ will be used.

[28] and many other studies have demonstrated that an ERP waveform across the parieto-central area of the skull is usually observed around 300 ms (the so-called P300 component) and is larger after the target event. The question to be addressed is whether it is possible to select time points at which ERP features can reliably detect which one of the two tones was presented to the subject and whether these time points are indeed around 300 ms as expected. This verification is of fundamental importance if the P300 component is to be considered as an electrophysiological marker for further assessment of psychiatric and neurological disorders.

It is conjectured that brain activations would be different over different time points depending on whether participants listen tones of 500Hz or 1000Hz. The data consist of T = 799 time points measured for J = 15 subjects and I = 2 conditions. Figure 1 shows ERP curves for subjects 1 (blue) and 2 (orange) for condition Hz500 (plain) and Hz1000 (dashed) as an example. This figure illustrates the large variability among subjects.

Y ijt = x 0ij µ t + a it + ε ijt , (1)
where x 0ij is a r-profile of baseline covariates for subject j, which does not depend on condition i, and ε ijt are random errors, normally distributed with mean 0, standard deviations σ = (σ t 1 , . . . , σ t T ) and correlation matrix R. In the above auditory oddball experiment, x 0ij = (1, δ j ) where δ j is just the 0 -1 variable which takes the value 1 if the subject is j. As shown in Figure 2, the within-condition standard errors s = (s t 1 , . . . , s t T ) vary quite a lot along time, where s 2 t is the degree-of-freedom corrected mean squared error at time t.

Moreover, the image plot of the within-condition correlation matrix among ERPs in Figure 3 shows that it is both structured by a strong-autocorrelation component and a block pattern.

Model (1) is similarly expressed as follows in matrix notations:

Y t = X 0 µ t + Xa t + ε t ,
where Y t = (Y 11t , . . . , Y 1,n 1 ,t , Y 21t , . . . , Y 2,n 2 ,t , . . . , Y I1t , . . . , Y I,n I ,t ) , X 0 is the n × r matrix whose rows are x 0ij , a t = (a 2t , . . . , a It ) is the (I -1)-vector of condition effect parameters, X Figure 3: Within-condition correlations in the auditory ERP dataset is a n×(I -1) matrix whose kth column X k is 0, except the entries between indices n 1 +. . .+n k +1 and n 1 +. . .+n k+1 which are equal to 1, and ε t = (ε 11t , . . . , ε 1,n 1 ,t , ε 21t , . . . , ε 2,n 2 ,t , . . . , ε I1t , . . . , ε I,n I ,t ) .

Under the within-subject independence assumption introduced above, for all heteroscedasticity and correlation pattern (σ, R), the generalized least-squares (GLS) estimate of a t coincides with the maximum likelihood (ML) estimator under normality:

ât = S -1 xx S xyt , (2)
where S xx is the (I -1) × (I -1) empirical variance matrix derived from the design matrix of model ( 1):

S xx = 1 n X P 0 X,
where

P 0 = I n -X 0 (X 0 X 0 ) -1 X 0 .
Similarly, S xyt is the (I -1)empirical covariance vector:

S xyt = 1 n X P 0 Y t .
Finally, the residual variance σ 2 t is estimated by the residual degree of freedom-corrected mean square of residuals errors:

s 2 t = Y t (P 0 -P )Y t n -(r + I -1) , (3)
where P is the n × n orthogonal projection matrix:

P = XP 0 X P 0 X -1 P 0 X .
Signal detection Signal detection can be expressed as the statistical test of the whole-signal null hypothesis H 0 : for all t, a t = 0.

The type-I error rate α for the former test is also the probability of declaring erroneously that, for at least one time t, a t = 0. Controlling the type-I error rate for the global test on the whole curves can therefore be viewed as equivalent to controlling the Family-Wise Error Rate (FWER), namely the probability of at least one false rejection, in the simultaneous testing of the collection of null hypothesis H 0t : a t = 0.

For each of the former tests, the following F-tests F t are appropriate:

F t = n I -1 â t S xx ât s 2 t . Let p t = 1-G I-
P H 0 t n I -1 â t S xx ât s 2 t ≥ f * ≤ α,
where

f * = G -1 I-1,n-(r+I-1) (1 -p * ).
However, the FWER-controlling multiple testing procedures are known to be conservative and the Bonferroni correction is not adapted to dependence among tests.

Multivariate analysis of variance

An alternative approach is to derive a single test statistics for H 0 by gathering the curves into a single vector Y = (Y t 1 , Y t 2 , . . . , Y t T ) . The multivariate analysis of variance model ( 1) is now expressed in the following expanded form:

Y = [I T ⊗ X 0 ]µ + [I T ⊗ X]a + ε, = X0 µ + Xa + ε, ( 4
)
where ⊗ is the Kronecker matrix product, µ = (µ t 1 , . . . , µ t T ) and a = (a t 1 , . . . , a t T ) . The above expanded linear framework is not homoscedastic. Indeed,

Var(ε) = [D σ RD σ ] ⊗ I n = V ε ,
where D σ is the T × T diagonal matrix which diagonal terms are σ t 1 , . . . , σ t T .

Therefore,

Var(â gls ) = 1 n [D σ RD σ ] ⊗ S -1 xx . (5)
Recent papers (see [START_REF] Bugli | Functional anova with random functional effects: an application to event-related potentials modelling for electroencephalograms analysis[END_REF][START_REF] Smith | Regression-based estimation of erp waveforms: I. the rerp framework[END_REF][START_REF] Smith | Regression-based estimation of erp waveforms: Ii. nonlinear effects, overlap correction, and practical considerations[END_REF]) suggest analysis of variance F-tests based on the above presentation of the ERP curves. Furthermore, [START_REF] Bugli | Functional anova with random functional effects: an application to event-related potentials modelling for electroencephalograms analysis[END_REF][START_REF] Smith | Regression-based estimation of erp waveforms: I. the rerp framework[END_REF][START_REF] Smith | Regression-based estimation of erp waveforms: Ii. nonlinear effects, overlap correction, and practical considerations[END_REF] take advantage of this framework to introduce a smooth non-parametric model for µ and a, using B-splines or wavelets. This modifies the design matrix of the model and correspondingly reduces the number of regression coefficients. Indeed, if ϕ stands for the T × S matrix associated to the basis functions, e.g. B-splines, (such that ϕ is = φ s (t i ), s = 1, . . . , S; i = 1, . . . , T ), then X0 = ϕ ⊗ X 0 and X = ϕ ⊗ X in expression (4) and the parameters µ and a are respectively Srand S(I -1) vectors of regression coefficients on the basis functions.

In the former functional analysis of variance approach, the handling of dependence in V ε remains the same as in the multivariate linear model ( 4).

The Generalized Least Squares (GLS) method is here appropriate since it provides the minimum variance linear unbiased estimator of the regression parameters provided (σ, R) is known. The corresponding F-statistics for H 0 : a = 0 is given by the following expression:

F gls = nâ gls ([D 1/σ R -1 D 1/σ ] ⊗ S xx )â gls . ( 6
)
which null distribution is χ 2 (I-1)T . However, F gls cannot be calculated because its expression depends on the unknown variance parameters. To handle this issue, the analysis of variance F-tests in [START_REF] Bugli | Functional anova with random functional effects: an application to event-related potentials modelling for electroencephalograms analysis[END_REF][START_REF] Smith | Regression-based estimation of erp waveforms: I. the rerp framework[END_REF][START_REF] Smith | Regression-based estimation of erp waveforms: Ii. nonlinear effects, overlap correction, and practical considerations[END_REF] are either based on the usual homoscedastic assumption V ε = σ 2 I nT or an assumption of a lag-1 autoregressive process for R:

R ρ =      1 ρ ρ 2 . . . ρ T -1 ρ 1 ρ . . . ρ T -2 . . . . . . . . . ρ T -1 ρ T -2 ρ T -3 . . . 1     
, where ρ is the lag-1 autocorrelation.

Homoscedasticity assumption Let Ŷ denote the vector of fitted values of Y , using the above GLS estimates of the model parameters. The corresponding F-test for the null hypothesis H 0 : a = 0 under the homoscedasticity and independence assumption is given by the following expression:

F ols = n â gls (I T ⊗ S xx )â gls σ2 , (7)
where

σ2 = (Y -Ŷ ) (Y -Ŷ ) T (n -(r + I -1)) . ( 8
)
Heteroscedasticity assumption Note that, if the above homoscedasticity assumption is relaxed, then V ε = D σ 2 ⊗ I n and the corresponding F-statistics:

F s = nâ gls (D 1/s 2 ⊗ S xx )â gls , (9)
where s 2 = (s 2 t 1 , . . . , s 2 t T ) is given by expression (3), is explicitly derived from the sum of the individual F t statistics:

F s = (I -1) T i=1 F t i . (10)
Therefore, in the present situation, heteroscedasticity can be straightforward accounted for by considering that the null distribution of F s is the distribution of (I -1)× the sum of T independent F I-1,n-(r+I-1) variables. In the following, the former null distribution is denoted FI-1,n-(r+I-1) .

Autoregressive covariance assumption Under the lag-1 autoregressive covariance assumption, let us introduce the matrix L ρ,σ defined as follows:

L ρ,σ =       -ρ σt 1 1 σt 2 0 . . . 0 0 -ρ σt 2 1 σt 3 . . . 0 . . . . . . . . . 0 0 . . . -ρ σt T -1 1 σt T       . If ε * = (L ρ,σ ⊗ I n )ε stands for the (n -1)T -vector of innovations, then it is straightforward checked that Var(ε * ) = (1 -ρ 2 )I (n-1)T .
Therefore, starting from a non-zero initial estimate ρ0 of ρ to derive L ρ,σ , an update is obtained by the variance of the corresponding ε * :

ρ2 1 = 1 - n i=1 T j=1 ε * 2 it j (n -(r + I -1))(T -1) .
Plugging-in this updated estimate of ρ in L ρ,σ provides new decorrelated residuals ε * and in turn a new estimate ρ2 of ρ. This defines an iterative algorithm which converges to the so-called Feasible GLS estimator ρ of ρ.

A corresponding feasible version F s,ρ of the F gls -statistics is obtained by plugging-in the estimates of σ t and R ρ in expression (6):

F s,ρ = nâ gls (D 1/s R -1 ρ D 1/s ⊗ S xx )â gls . (11)
The exact null distribution of F s,ρ is untractable. The analysis of variance F-tests are usually based on the asymptotic approximation of the null distribution by a F (I-1)T,(n-(r+I-1))T degrees of freedom.

Finally, we suggest here to consider the F-test based on F s (see expression ( 9)). Noticing that

Cor(F t i , F t j ) ≈ n→+∞ ρ 2|i-j| ,
the null distribution of F s under the lag-1 autoregressive covariance assumption can be approximated by the distribution of (I -1)× the sum of lag-1 autocorrelated F I-1,n-(r+I-1) variables, with auto-correlation ρ 2 . In the following, this null distribution is denoted FI-1,n-(r+I-1) (ρ).

Impact of a model misspecification on signal detection

In order to compare the different signal detection strategies described in the present section, we simulate datasets which dimensions and distribution mimic those obtained from the real auditory oddball dataset. The original dataset contains 30 ERP curves observed in two conditions, 15 curves in each condition on the time interval [0;400ms], at the frequency of 1 observation every half millisecond. We focus here on the significance of the difference curve. In the present simulation settings, it is assumed under H 1 that the expected ERP curve are constantly zero in condition 1, whereas it describes the waveform plotted in Figure 4 in condition 2. In the next section, we will show that, in the Rare-and-Weak paradigm of [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF], the present situation falls into the so-called estimable region in terms of sparsity and weakness of the signal in the phase diagram under independence.

1000 datasets are generated, under H 0 and under H 1 , with n = 30 rows and T = 800 columns, with standard deviation profiles s and a lag-1 autoregressive correlation matrix, with autocorrelation ρ = 0.99. Note that the estimated auto-correlation in the auditory oddball experiment is 0.997. On each simulated dataset, the following F-tests are implemented:

1. The F-test F ols (see expression [START_REF] Causeur | ERP: Significance analysis of Event-Related Potentials data[END_REF]) based on the homoscedasticity and independence assumption, with the null distribution F T,(n-2)T ;

2. The same F-test as above calculated in the model including a B-spline smoothing of the regression coefficients. The corresponding F-test is implemented using the functions bam and anova.gam in the R package mgcv with the default option for the choice of the smoothing parameter (minimizing the Generalized Cross Validation criterion). Note that the null distribution of the F-statistics is still a Fisher distribution, whose degrees-of-freedom accounts for the smoothness of the fit by means of the trace of the smoothing matrix;

3. The F-test F s,ρ (see expression [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF]) with the null distribution F T,(n-2)T ;

4. The same F-test as above calculated in the model including a B-spline smoothing of the regression coefficients. The corresponding F-test is also implemented in the function bam, using arguments especially designed to introduce a lag-1 autocorrelation prior;

5. The F test F s (see expression 9) with the null distribution FI-1,n-(r+I-1) (ρ). The calculation of the p-values are based on Monte-Carlo estimation of the former null distribution.

Figure 5 displays the empirical null probability distribution function of the p-values for five Fstatistics. The plots confirms that the main impact of the model misspecification here is a too liberal F ols -test, not controlling the type-I error rate. The concern is even more obvious when the fit is based on regression spline smoothing. This is consistent with the observation reported in many papers, including [START_REF] Bugli | Functional anova with random functional effects: an application to event-related potentials modelling for electroencephalograms analysis[END_REF], that the confidence bands based on the homoscedasticity and independence assumption are strikingly too narrow. Another conclusion is that the approximate null distribution of F s,ρ is obviously wrong, leading to a too liberal test, both without and with spline smoothing. The F-test based on the sum of individual F-tests, with an ad-hoc correction of the null distribution, turns out to control the type-I error rate at the desired level. Correspondingly, Figure 6, which shows the empirical non-null probability distribution functions of the same five F-tests, demonstrates the lesser ability of F s to detect the signal.

Finally, the above simulation study shows that, provided we can efficiently derive the null distribution of the F-test under the dependence assumption, a signal detection strategy based on the individual F t statistics can help achieving a correct control of the type-I error rate. The Higher Criticism Thresholding method presented hereafter is also exclusively based on the vector of F t statistics. Note that, in the above simulation study, the ERP time dependence is supposed to be described by a lag-1 autoregression dependence structure, which is far from true (see Figure 3 for the image plot of the correlation matrix calculated with the auditory oddball dataset).

In the following, to be consistent with the framework defined by [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF], we will focus on the situation where I = 2. Consequently, we will prefer dealing with the T-tests T = (T t 1 , . . . , T t T ) deduced by taking the signed square root of the F-tests F = (F t 1 , . . . , F t T ) . The Higher Criticism (HC) method is designed to detect signals, namely to provide a single p-value for the significance of a signal, based on the vector of significance tests for the individual coordinates of this signal. It was initially proposed by [START_REF] Tukey | T13 N: the higher criticism[END_REF] in his course notes at Princeton University and brought up to date by [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF]. HC is designed to be optimal in a rare and weak model for the signal, defined in [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF] as a sparse normal mixture model for z-scores. In this paper, we use the equivalent notation of [14]

T = µ + E , E ∼ N (0, R), ( 12
)
where R is set to the identity matrix and µ is a sparse T -vector with a proportion of non-null features 0 ≤ ε T ≤ 1 and non zero entries set to the same signal amplitude A T ≥ 0. Note that the normality assumption introduced above holds for most ERP studies in which the tests for the association between the ERPs and the response variable can be handled by t-tests for the significance of a single parameter ( [START_REF] Causeur | A factor-adjusted multiple testing procedure for erp data analysis[END_REF], [START_REF] Guthrie | Significance testing of difference potentials[END_REF]). The RW model supposes a subtle situation in which the signal is rare, as ε T = T -β with β ∈ ( 1 2 , 1) and weak A T = 2r log(T ), 0 < r < 1. When parameters (β, r) are known, explicit detection boundaries can be calculated for the RW setup under independence, which separate the space (β, r) into undetectable, detectable ( [START_REF] Ingster | Minimax detection of a signal for p n balls[END_REF]) and estimable ( [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF]) regions and are summarized in a so-called phase diagram (see Figure 7): if r > ρ * D (β), the signal is detectable which means that the sum of type I and type II errors of the Neymann-Pearson Likelihood Ratio test (LRT) tends to 0 when the number of tests tends to infinity. If r > ρ * E (β), the signal is not only detectable but the support of non null coordinates is identifiable. If (β, r) are out of these bounds, the signal is undetectable so that the sum of type I and type II errors for any test tends to 1 as the number of tests tends to infinity. The bounds ρ * D of detectability and ρ * E of estimability have the following expressions:

ρ * D (β) = β -1 2 if 1 2 < β ≤ 3 4 (1 - √ 1 -β) 2 if 3 4 < β < 1, , ρ * E (β) = β.
In the example introduced in the simulation study of the previous section, the difference curve between two conditions shows a waveform. Therefore, the corresponding curve of r coefficients is not constant on intervals of nonzero signals as in the above RW setup, as shown in Figure 8. However, note that a proportion of 3% (β = 0.52) of the signal has a r coefficient larger than 0.7.

The corresponding RW situation (β = 0.52, r = 0.7) falls into the estimable region in the phase diagram displayed in Figure 7. Signal detection by HC starts from the observation that the presence of a signal generates a discrepancy between the empirical distribution of the p-values and the theoretical null uniform distribution U[0; 1]. Indeed, it is expected that, for non-zero signal coordinates, the corresponding t-th increasingly ordered p-value p (t) is such that p (t) t/T . This discrepancy is measured by a Higher Criticism function HC(p (t) , t). The various HC methods are based on different definitions of the HC objective function: [START_REF] Donoho | Higher criticism thresholding: Optimal feature selection when useful features are rare and weak[END_REF] and [START_REF] Klaus | Signal identification for rare and weak features: higher criticism or false discovery rates[END_REF] standardize the difference t/T -p (t) by the standard-deviation of ordered p-values under the null, V ar(p (t) ) = t T (1 -t T ) and define the HC objective function as:

HC(p (t) , t) = √ T t/T -p (t)
t/T (1 -t/T ) .

Alternatively, [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF] and [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF] suggest to standardize the difference t/T -p (t) by p (t) (1 -p (t) ) and define the HC objective function as:

HC(p (t) , t) = √ T t/T -p (t) p (t) (1 -p (t) )
.

Using arguments from the theory of empirical processes, [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF] shows that

HC * √ 2 log log T → 1, in Probability,
where HC * is the maximum of HC(p (t) , t). It is deduced that the type-I error rate of any signal detection rule of the form HC * ≥ (1 + a) √ 2 log log T , for a > 0, tends to 0 as T tends to infinity. Moreover, [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF] shows that the former tests are optimal in the sense that, for any Rare-Weak situation (r, β) within the detectable region, the type-II error rate of the tests based on HC * also tends to 0 when T tends to infinity. This property is known as the optimal adaptivity of Higher Criticism. In order to obtain a testing strategy which controls the type-I error rate, [START_REF] Cai | Optimal detection of heterogeneous and heteroscedastic mixtures[END_REF] suggest to derive the value of a from a Monte-Carlo estimation of the null distribution of HC * . We have implemented the former method in the simulation study introduced in the previous section. Figures 9 and10 display the empirical null (resp. non-null) probability distribution function of the p-values for the HC statistics and the F σ-statistics, which accounts for dependence. The plots show the two methods have similar performance.

Note that, in practice, the HC function is often only maximized over a subset I of features with small p-values:

HC * = max t∈I HC(p (t) , t).
Thus, another major variation among the different HCT methods is on the definition of the subset I: for [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF] and [START_REF] Donoho | Higher criticism thresholding: Optimal feature selection when useful features are rare and weak[END_REF], I = {t, 1/T ≤ t/T ≤ α 0 } where α 0 ∈ [0, 1] is a preset proportion of small pvalues. Alternatively, [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF] proposes to avoid the extremely small p-values by taking I = {t, 1/T ≤ p (t) ≤ α 0 }. At last [START_REF] Klaus | Signal identification for rare and weak features: higher criticism or false discovery rates[END_REF] proposes a global maximization I = {1, . . . , T }. In the above simulation of the HC statistics, the value α 0 = 0.1 has been chosen, to be consistent with the recommendation in [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF].

In the next sections, the standard HCT refers to the following criterion:

HC * = max 1/T ≤t/T ≤α 0 √ T t/T -p (t) t/T (1 -t/T ) (13)
which appears to be the most widely used in the literature.

In a feature selection objective for supervised classification, [START_REF] Donoho | Higher criticism thresholding: Optimal feature selection when useful features are rare and weak[END_REF] highlight the link between maximizing the HC function and minimizing the error rate of the classifier implemented on the features selected as follows: if the maximum value HC * of the HC function is reached at index t * , then the subset of selected features is given by {t, p t ≤ p (t * ) }. [START_REF] Donoho | Higher criticism thresholding: Optimal feature selection when useful features are rare and weak[END_REF] demonstrates that the former Higher Criticism Threholding (HCT) method outperforms methods based on the False Discovery Rate (FDR) thresholding. However, [START_REF] Klaus | Signal identification for rare and weak features: higher criticism or false discovery rates[END_REF] shows the equivalence between HCT and the definition of a cutoff on the local FDR.

Empirical non-null pdf of p-values p-value

Higher Criticism Sum F Theoretical uniform

HCT under dependence

Even if the optimality of HCT is known to be robust to weak dependence, the simulation study above suggests that it would be improved by explicitly taking into account dependence. The existing variants of HC accounting for dependence are mainly based on a decorrelation step of the vector T of test statistics. Initially proposed by [START_REF] Zuber | Gene ranking and biomarker discovery under correlation[END_REF] in a multiple testing framework but also used as a variable selection procedure in a supervised classification issue in [START_REF] Ahdesmäki | Feature selection in omics prediction problems using cat scores and false non-discovery rate control[END_REF], the Correlation-Adjusted T-scores are defined as:

τ adj = R -1/2 τ,
where τ are standard t-scores arising from a multiple testing procedure. The correlation matrix R is estimated by a James-Stein estimator R γ = γI m + (1 -γ)R as proposed by [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF]. Interestingly, [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF] provide an analytical expression for the estimation of γ. The inverse square root of R γ is computationally feasible in high dimension by observing that the matrix Z = 1 γ R γ can be decomposed as Z = I m + U M U where M is a definite positive symmetric matrix and U an orthonormal basis matrix. Hence,

Z α = I m -U (I r -(I r + M ) α )U ( 14
)
where r is the rank of the empirical correlation matrix. Even if it is efficient in practice, especially for many applications in genomic data analysis, the estimator proposed by [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF] seems to fail to catch the dependence structure of ERP data. Indeed, the right panel of Figure 11 displays an image of the residual correlation matrix presented in Figure 3, which shows that the blocks of auto-correlation are underestimated.

Similarly, [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF] propose the innovated HCT (iHCT) based on a preliminary whitening of T . To this end, they introduce a matrix U m such as U m RU m = I m and apply standard HCT to decorrelated test statistics in the manner of [START_REF] Zuber | Gene ranking and biomarker discovery under correlation[END_REF]. They recommend to use the inverse of the Cholesky factorization of the matrix R, which is computationally unstable in high dimension. As the support of the signal is slightly shifted by this linear transformation, the authors propose to restore the support by smoothing techniques applied on U m , in the context of an auto-regressive dependence pattern.

Factor-Innovated Higher Criticism Thresholding

The proposed method is conceptually similar to innovated HCT and CAT-scores except that the decorrelation is based on a factor model. In many areas in which high-throughput devices are used to observe complex systems, such as neuroscience (see [START_REF] Causeur | A factor-adjusted multiple testing procedure for erp data analysis[END_REF]) or genomics (see [START_REF] Blum | A factor model to analyze heterogeneity in gene expression[END_REF], [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF], [START_REF] Leek | A general framework for multiple testing dependence[END_REF], [START_REF] Sun | Multiple hypothesis testing adjusted for latent variables, with an application to the AGEMAP gene expression data[END_REF]), this model is now widely used to implement data reduction techniques. In these papers, the authors suggest to cope with dependence by adjusting data from the effects of a moderate number of latent variables, which catch the dependence structure.

Whitening using latent factors

Factor model assumes the conditional independence of the variables given latent factors. Indeed, conditionally on latent factors, Expression (1) becomes:

Y ijt = x 0ij µ t + a it + b t z ij + e ijt , (15)
where z ij = (z ij1 , . . . , z ijq ) is normally distributed with mean 0 q and variance I q and e ij = (e ij,t 1 , . . . , e ij,t T ) are independent random error vectors such that e ij ∼ N T (0 T , Ψ Σ ) where Ψ Σ is a diagonal matrix of specific variances described hereafter. q is the number of latent factors, q T . The algorithms developed in a large number of papers (see [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF], [START_REF] Causeur | A factor-adjusted multiple testing procedure for erp data analysis[END_REF] [26], [START_REF] Leek | Capturing heterogeneity in gene expression studies by surrogate variable analysis[END_REF], [START_REF] Perthame | Stability of feature selection in classification issues for high-dimensional correlated data[END_REF]) are designed for the estimation of the latent factors in order to adjust data from their effects. In this paper, we intentionally focus on another implication induced by the rank-reduced factor structure of the covariance matrix Σ. Indeed, factor model assumes equivalently that the covariance matrix Σ can be decomposed as follows:

Σ = Ψ Σ + B Σ B Σ ,
where B Σ is a T × q matrix of loadings, which describes the common dependence shared by the variables and Ψ Σ = diag(Ψ Σ,1 , . . . , Ψ Σ,T ) is a diagonal matrix, which defines the specific variance. A similar factor structure holds for residual correlation matrix:

R = Ψ + BB , where Ψ Σ = D σ 2 Ψ, Ψ = diag(Ψ 1 , . . . , Ψ T ) and B Σ = D σ B.
The parameters (Ψ, B) are estimated by an EM-algorithm proposed by [START_REF] Sheu | Accounting for time dependence in largescale multiple testing of event-related potential data[END_REF] and implemented in the R package ERP available online on the CRAN ( [START_REF] Causeur | ERP: Significance analysis of Event-Related Potentials data[END_REF]). The authors provide an adaptive estimation of the dependence structure under the factor model [START_REF] Ingster | Some problems of hypothesis testing leading to infinitely divisible distribution[END_REF] based on a simultaneous estimation of the signal and covariance matrix. Moreover, the method is designed to take advantage of a prior knowledge on time intervals where the signal does not occur in order to denoise the least-squares estimation of the signal. The number of factors q is estimated by the variance inflation criterion proposed by [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF].

The left panel of Figure 11 illustrates that factor modeling is flexible enough to recover the dependence pattern observed in the auditory oddball ERP dataset. The Figure shows an image plot of the estimation of correlation matrix presented in Figure 3 by a factor model with q = 6 factors.

Detection boundary

Deriving exact boundaries in an arbitrary dependence pattern is impossible since the possibly nonhomogeneous dependence along time can generate a different detectability and estimability of the signal for some subset of features. By the way, [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF] do not propose a closed-form phase diagram under dependence but explicit lower and upper bounds for the detection boundaries under some dependence patterns. In this section, we also propose to extend the phase diagram presented in Section 3 to the dependent case in a factor model framework.

Provided the test statistics are linear transformations of the response variables, which is generally the case for t-tests, the covariance structure of the test statistics is inherited from the residual correlation of those response variables. If Z stands for the same linear transformation applied to the n × q matrix Z, which rows are z ij , then model [START_REF] Guthrie | Significance testing of difference potentials[END_REF] for the test statistics T becomes:

T = µ + BZ + E , E ∼ N (0, Ψ).
Equivalently,

Ψ -1/2 (T -BZ ) = Ψ -1/2 µ + E * , E * ∼ N (0, I T ). ( 16
)
The above equation offers a RW setup similar to model ( 12), in which the decorrelated test statistics Ψ -1/2 (T -BZ ) are independent with equal variance 1, given the factors Z. Upper and lower bounds for the detection boundary of a signal with constant amplitude A T = 2r log(T ) for a proportion ε T of signal coordinates are deduced. Indeed, depending on the distribution of the specific variances Ψ t along the features, the amplitude of the transformed signal Ψ -1/2 µ can take any values between γ 0 = A T / √ Ψ max = 2r log(T ) and γ 0 = A T / √ Ψ min = 2r log(T ), where Ψ min = min t (Ψ t ), Ψ max = max t (Ψ t ), r = r/Ψ min and r = r/Ψ max . The conditions on (r, β) for detectability established under independence can be invoked here, with two different values r and r for the strength parameter r, providing respectively the lower and upper bounds for detectability. First, if r satisfies the following condition:

r > Ψ max ρ * D (β), (17)
then all the signal is detectable. Moreover, if r does not satisfy the following condition,

r > Ψ min ρ * D (β), (18)
then the signal can not be detected. Between these two bounds, detection of the signal is uncertain, depending on the correspondence between the profile of specific variances and the support of the signal.

The same reasoning for estimability gives similar forms for estimation boundaries. At last, we propose to define a factor RW setup in which the sparsity of the signal is given by

ε T = T -β with β ∈ (1/2, 1),
and the strength of the signal by

A T = 2rlog(T ) with 0 < r < Ψ max .
Under independence, the signal strength is set to 2r log(T ) with 0 < r < 1 to make signal detection challenging. Indeed, 2 log(T ) is the expectation of the maximum of T independent normal variables with mean 0 and variance 1. Note that signals with amplitudes 2r log(T ), with Ψmax < r ≤ 1, are considered as weak under independence but not under dependence. Consistently, this observation is also reported by [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF], who explains that correlated designs are actually favorable for signal detection with procedures accounting for dependence. Moreover, this is illustrated by Figure 12, which shows the detectability and estimability phase diagrams obtained with the estimated specific variances in the auditory oddball experiment. Indeed, the plot confirms that both the detectability and the estimability regions are wider under dependence.

The bounds in Expressions ( 17) and ( 18) are consistent with the ones established by [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF]. Indeed, Theorems 3.1 and 4.2 of the former paper give the same expression as above for the detection bounds, where γ 0 and γ 0 are defined as follows:

γ 0 = lim inf T →+∞ max √ T ≤k≤T - √ T R -1 kk , γ 0 = lim sup T →+∞ max √ T ≤k≤T - √ T R -1 kk ,
where (R 11 , . . . , R kk , . . . , R T T ) stands for the diagonal of matrix R. In the present situation, the following sandwich inequality holds:

γ 0 = 1 max(Ψ) ≤ diag(R -1 ) ≤ 1 min(Ψ) = γ 0 ,
which gives expressions ( 17) and (18).

Factor innovated HCT

As in [START_REF] Zuber | Gene ranking and biomarker discovery under correlation[END_REF] and [START_REF] Hall | Innovated higher criticism for detecting sparse signals in correlated noise[END_REF], we propose a factor-innovated HCT (F-iHCT) which consists in applying standard HCT to decorrelated test-statistics. In the factor RW setup defined above, we define the following square root L of the correlation matrix R, such as R = LL :

L = (I m -U [(I q + [I q + D 2 ] 1/2 ) -1 + I q ] -1 U )Ψ -1/2 ,
where U and D comes from the singular value decomposition of the standardized loadings Ψ -1/2 B = U DV . Note that this formula only requires inversion and square root computation of diagonal matrices. The decorrelated test statistics are therefore deduced:

T * = L T
and the corresponding p-values are denoted p * t . HCT is then applied on the collection of p-values p * = (p * 1 , . . . , p * T ).

Simulation results and real data analysis

The properties of the proposed method are now compared to standard HCT and to some other methods based on decorrelation by innovations or by adjustment for effect of latent variables. The comparison is made through a simulation study and an application on the auditory oddball experiment.

Simulation study

Simulation settings The properties of F-iHCT are now investigated through simulations. 1,000 datasets with dimensions 30 × 799 are generated according to a multivariate normal distribution.

Both the correlation structure and the within-condition variances are estimated from the auditory oddball ERP data introduced in section 2 (see Figure 3 and Figure 2). This simulation plan mimics the observed data on the oddball experiment by dimensions and covariance structure except that the true signal is known. Each dataset is split into two balanced groups. The normal distribution has expectation zero for the first 15 subjects (group 1) and the expectation for the 15 last subjects (group 2) is plotted on According to the RW setup, the present combination of r and β characterizes a not very sparse signal, with a weak to strong strength.

Methods Four methods are compared in this simulation study and are described hereafter. As in [START_REF] Donoho | Higher criticism thresholding: Optimal feature selection when useful features are rare and weak[END_REF] and [START_REF] Donoho | Feature selection by higher criticism thresholding achieves the optimal phase diagram[END_REF] the variable selection step by different versions of HCT, with two possible values for α 0 (α 0 = 0.1 or 0.5), is followed by a supervised classification by Diagonal Discriminant Analysis on the subset of selected variables. The following methods are compared:

1. Variable selection by standard HCT on raw p-values, classification by Naives Bayes (see [START_REF] Bickel | Some theory for Fisher's Linear Discriminant function, naive Bayes, and some alternatives when there are many more variables than observations[END_REF]) denoted by Standard HCT;

2. Variable selection by HCT on CAT-scores ( [START_REF] Zuber | Gene ranking and biomarker discovery under correlation[END_REF][START_REF] Ahdesmäki | Feature selection in omics prediction problems using cat scores and false non-discovery rate control[END_REF]), classification by diagonal Shrinkage Discriminant Analysis (SDA, see [START_REF] Ahdesmäki | Feature selection in omics prediction problems using cat scores and false non-discovery rate control[END_REF]) denoted by CAT-scores;

3. Variable selection by Factor-innovated HCT, classification by conditional Bayes classifier (proposed by [START_REF] Perthame | Stability of feature selection in classification issues for high-dimensional correlated data[END_REF]) denoted by F-iHCT; ), classification by conditional Bayes classifier (see [START_REF] Perthame | Stability of feature selection in classification issues for high-dimensional correlated data[END_REF]) denoted by AFA.

For all the methods described above, the proportion of signal recovery, called precision, the false discovery rate (FDR), the number of selected features and the prediction error rate are computed. For all datasets, variable selection and estimation of classification rule are performed on training data (including the optimization of meta-parameters) and prediction error is computed on testing data.

Results When α 0 is well specified regarding to the proportion of null hypothesis (α 0 = 0.1), Figure 14 shows that selection by CAT-scores appears to be the most efficient to catch weak signals, with both the smallest FDR and the largest precision for small amplitudes of signal. Even if CATscores does not achieve the best performance for large signal strengths, the FDR, precision and number of selected variables are remarkably stable. Standard HCT seems robust to dependence as the method performs well in term of FDR but its precision is small regarding methods based on decorrelation. Moreover, the number of selected variables is also small, which suggests that HCT is conservative under dependence. Lastly, classification by Naive Bayes fails as the error rates are the largest for weak to moderate strengths of signal. Variable selection and classification procedures based on the factor model assumption (AFA and F-iHCT) provide the best results both in terms of false positive, recovery of the signal and prediction error. FDR turns out to be small for moderate to high signal strengths and a correct power of signal identification is achieved.

As shown on Figure 15, all methods are affected by a misspecification of α 0 parameter (α 0 = 0.5). CAT-scores method selects too many variables so it achieves a good precision at the cost of large FDR. F-iHCT and AFA perform well in classification despite the number of false positives. Classification and precision rates of standard HCT are also affected by this wrong choice of α 0 .

ERP data analysis

The 4 methods compared in the simulation study are now applied on the auditory oddball ERP data presented in Section 2. It is reminded that the purpose associated to the oddball experiment is to predict a new label from ERP curves. For each method, the number of selected features and the prediction error are computed. As the number of observations is small, the classification error is computed by leave-one-out cross-validation (CV).

Table 1 shows the number of selected features by the 4 compared selection methods for different values of α 0 . Whatever the value of α 0 , the 4 methods select the same number of features. The difference between all methods lies on the indices of selected features as they do not selected the same time points, as shown as an example on Figure 17.

Figure 16 presents the cross-validated error rates for several values of α 0 . For values of α 0 larger than 0.15, standard HCT is stable and performs rather well. For more sparse models, standard HCT reaches larger error rates and is improved by decorrelation methods based on a factor model assumption (AFA and F-iHCT). The performance of the CAT-scores method varies slightly depending on α 0 . The curves of F-iHCT and AFA are erratic for values of α 0 smaller than 0.125 but they stabilize when α 0 increases. For values of α 0 larger than 0.275, F-iHCT and AFA appear to be the most effective methods as they perfectly classify data. Nevertheless, one can notice that for equal error rates, the two methods do not select the same features as shown on the bottom of Figure 17.

Figure 17 shows the curve of the mean difference among the two groups and the time points selected by the 4 compared methods for α 0 = 0.125 (top) and α 0 = 0.275 (bottom). These values of α 0 are chosen because they provide two levels of sparsity but comparably small CV error. As expected, time points after 300ms are selected by all methods which is consistent with the literature but time points around 100 ms also appear to be significant. To conclude, this study highlights the importance of the choice of α 0 in HCT, since it has a significant impact on the sparsity and classification rates of the models for all methods. The method denoted by AFA which performs decorrelation by adjustment of covariates for the effect of latent factors seems to be the most suitable method in this example, even if F-iHCT also achieves interesting classification rates.

Discussion and conclusion

This article addresses signal detection and identification in Event-Related-Potentials (ERP) data. It is motivated by an ERP study in the oddball paradigm, where two classes of stimuli are presented to subjects, one occurring frequently (standard) and the other occurring infrequently (target). The aim is to identify time intervals of ERP curves which could be markers for the rare-frequent difference and to derive a classification rule. When the selection statistics are independent, Higher Criticism is known to be optimal as a signal detection method in the Rare-and-Weak setup initially introduced by [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF]. Higher Criticism Thresholding is moreover efficient to estimate the support of the signal. The Rare-and-Weak setup is conceptually adapted to the ERP signal detection issue. However, selection statistics in ERP signal identification issues are characterized by a strong and complex dependence structure.

In ERP data analysis, signal detection is usually handled by F-tests for the overall nullity of the signal on the whole time interval of observation. Under heteroscedasticity and independence, we show that the corresponding Generalized Least Squares (GLS) F-test is expressed as a sum of the individual F-tests for each feature, which gives an explicit expression for the null distribution. Similarly, under a dependence pattern structured by a lag-1 auto-correlation, we show that the former GLS F-test can also be implemented, provided the null distribution accounts for the auto-correlation between the individual F-tests. A simulation study in which dependence among simulated features is structured by a strong lag-1 auto-correlation, demonstrates that the former test controls the type-I error rate, which is far from true for other F-tests, such as those obtained in Functional Analysis of Variance approaches.

The present paper proposes a variant of the HCT procedure which takes advantage of a factor model assumption to decorrelate the test statistics. Indeed, this framework provides algebraical tools to derive an inverse square root of the correlation matrix involved in the computation of innovations. Moreover, similarly to [START_REF] Ingster | Some problems of hypothesis testing leading to infinitely divisible distribution[END_REF] and [START_REF] Donoho | Higher criticism for detecting sparse heterogeneous mixtures[END_REF], a phase diagram under a general dependence assumption is deduced and the Factor Innovated HCT, a decorrelated HCT based on innovations, is proposed.

The method is assessed by a more intensive simulation study, in which the dependence structure of simulated datasets mimics the observed correlation pattern in the auditory oddball ERP dataset. Estimated signal q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q standard HCT q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q CAT-scores q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q F-iHCT q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q AFA 0 100 200 300 400 Estimated signal q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q standard HCT q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q CAT-scores q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q F-iHCT q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q AFA

Figure 1 :

 1 Figure1: ERP curves for subjects 1 (blue lines) and 2 (orange lines) of the auditory oddball experiment in conditions Hz500 (plain) and Hz1000 (dashed)

Figure 2 :

 2 Figure 2: Within-condition standard errors in the auditory ERP dataset

Figure 4 :

 4 Figure 4: Expected difference ERP curve in the simulation study

Figure 5 :

 5 Figure 5: Empirical null probability distribution function of the F-test p-values under lag-1 autoregressive correlation

Figure 6 : 1

 61 Figure 6: Empirical non-null probability distribution function of the F-test p-values under lag-1 autoregressive correlation

Figure 7 :Figure 8 :

 78 Figure 7: Phase diagram of undetectable, detectable and estimable regions under independence

Figure 9 :

 9 Figure 9: Empirical null probability distribution function of the p-values of the HC and F σ statistics under lag-1 autoregressive correlation

Figure 10 :

 10 Figure 10: Empirical non-null probability distribution function of the p-values of the HC and F σ statistics and under lag-1 autoregressive correlation

Figure 11 :

 11 Figure 11: Estimation of within-condition correlation matrix by AFA algorithm ([23]) under a factor model assumption with 6 factors (left panel) and by the shrinkage estimator of [22] (right panel)

Figure 12 :

 12 Figure 12: Detection (upper plot) and identification (lower plot) boundaries for the auditory oddball ERP data. The grey line is the boundary under independence and black lines are upper and lower bounds of detectability (top) and identification (bottom) of a signal under dependence.

Figure 13 .

 13 The difference curve is therefore a waveform with various amplitudes and the indices of non null features are in [150ms, 200ms]. 1, 000 training datasets are generated for each signal strength. Eight corresponding testing data of size 1000 × 799 with two balanced groups are also generated according to the same simulation plan for a prediction purpose. The RW model parameters for this simulation plan are ε T = 12% and the maximum amplitude of signals is expressed as A T = 2rlog(T ) with r taking 8 equally distributed values in [0.004; 0.688].

4 .Figure 13 :

 413 Figure 13: Simulation study -Signal amplitude curves along time

Figure 14 :

 14 Figure 14: Results of the simulation study depending on signal strength and α 0 = 0.1: False Discovery Rate (top left), Precision (top right), Number of selected features (bottom left), Prediction error (bottom right).

FigureFigure 16 :

 16 Figure Results of the simulation study depending on signal strength and α 0 = 0.5: False Discovery Rate (top left), Precision (top right), Number of selected features (bottom left), Prediction error (bottom right).

Figure 17 :

 17 Figure 17: Real data study -Signal estimation (grey line) and significant time points (blue points) selected by standard HCT, HCT performed on CAT-scores, factor innovated HCT and HCT performed on p-values provided by AFA for α 0 = 0.125 (top) and α 0 = 0.275 (bottom) on auditory ERP experiment

  on the p-values for the simultaneous tests of H 0t , t = {t 1 , . . . , t T }. In the present situation of independent tests, one of the most famous method is the Bonferroni correction: choosing p * = α/T ensures that the Family-Wise Error Rate (FWER) is lower than α:

	1,n-(r+I-1) (F t ), where G I-1,n-(r+I-1) (.) is the probability distribution function
	of the Fisher distribution with (I -1) and (n -(r + I -1)) degrees of freedom, denote the corre-
	sponding p-value. Many multiple testing methods can be used to determine a rejection threshold p

*

Table 1 :

 1 Real data study -Number of selected time points for ERP auditory experiment

	α 0	0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
	Number of variables	40	80	120 160 200 240 280 320 360 400

This study illustrates that methods based on decorrelation of HCT under a factor model assumption perform well both in selection and classification. Indeed, F-iHCT and AFA, which is based on data adjustment for latent factors, achieve similar results. When applied on the oddball ERP dataset, all compared methods (standard HCT, CAT-scores, F-iHCT and AFA) select time points around 300 ms as expected but lower cross-validated error rates are achieved by the two methods based on a factor model assumption. Both application on simulations and on real data reveal the sensitivity of HCT to the choice of α 0 parameter, which determines the sparsity of the model.

The above promising results in signal identification should now be exploited in the construction of optimal signal detection strategies in ERP study. Indeed, provided the null distribution of the F-iHC statistics can be estimated, it could be interesting to compare the subsequent testing method to the class of F-tests.


Kun Jiang

Alessandro Corrêa Victorino

Ali Charara

Real-Time Estimation of Vehicle's Lateral Dynamics at Inclined Road Employing Extended Kalman Filter

In order to enhance road safety, this article presents a robust methodology for estimation of vehicle's lateral dynamics states, including velocities, sideslip angles, accelerations and tire forces. Awareness of vehicle dynamics states is essential for most intelligent vehicle control systems. However some of these states are not available by direct sensor measurement for technical and economical reasons. Even those available sensors are always coupled with noises or errors. In order to become reliable, many intelligent systems use very expensive sensors, which are not available for ordinary passenger cars. Therefore, it is quite interesting to have an accurate estimation of all dynamics states based on low-cost sensors. The main contribution of this article is providing a reliable estimator of vehicle dynamics states in presence of sensor noises and irregular roads. A best-wheel selection approach is applied in order to have a better understanding of vehicle's driving situation and enable us to select the best models for estimating dynamics states. The proposed estimator is based on the algorithm of extended Kalman Filter, which could minimize the variance of estimation errors. The estimator is implemented in our experimental vehicles and tested at banked track. Experimental results validate and prove the feasibility of this approach.

I. INTRODUCTION

Accurate knowledge of vehicle dynamics states (such as velocities, yaw rate, sideslip angle and wheel ground forces) is essential for advanced driver assistance systems (ADAS), such as adaptive cruise control (ACC), electric stability control (ESC). The control of longitudinal and lateral dynamics stability can be realized by longitudinal speed control and yaw rate control. The longitudinal speed and yaw rate can be measured directly, thus the control signal could be obtained easily [START_REF] Zhu | Real-time yaw rate prediction based on a non-linear model and feedback compensation for vehicle dynamics control[END_REF][START_REF] Klomp | Longitudinal velocity and road slope estimation in hybrid electric vehicles employing early detection of excessive wheel slip[END_REF]. However, a system directly controlled by measurement will be probably unstable. The measurement may contain lots of undesired noises and errors. Moreover, the measurement usually cannot cover all the dynamics states of vehicle. In this situation, due to incorrect and inadequate comprehension of vehicle dynamics states, the control system is not able to make the best decisions to ensure vehicle's safety.

Consequently, the objective of this paper is to develop a robust estimator to provide more accurate and more detailed information about vehicle dynamics states than the direct measurement. The dynamics states to be observed include longitudinal speeds, lateral speeds, yaw rate, tire forces and the sideslip angle of vehicle.

Vehicle longitudinal speed sometimes could be regarded as wheel speed, which is measured by wheel speed sensors. Nevertheless, wheel speeds could become unrelated to the vehicle velocity when driving in slippery conditions. In order to compensate the errors caused by wheel slip, GPS signal is used to correct the data of wheel speed sensor as a fusion of data [START_REF] Dahmani | Vehicle dynamic estimation with road bank angle consideration for rollover detection: theoretical and experimental studies[END_REF][START_REF] Gadola | Development and validation of a Kalman filter-based model for vehicle slip angle estimation[END_REF][START_REF] Li | Reliable vehicle sideslip angle fusion estimation using low-cost sensors[END_REF][START_REF] Ouahi | Simultaneous observation of the wheels' torques and the vehicle dynamic state[END_REF][START_REF] Jihan | Vehicle sideslip and roll parameter estimation using GPS[END_REF]. Besides the longitudinal speed, sideslip angle is also an important parameter to describe vehicle dynamics. The higher the sideslip angle is, the more easily drift out happens for a vehicle, especially on a low friction road. As a result, the sideslip angle has been used in many control system to improve stability. Two common techniques to estimate vehicle sideslip angle are to directly integrate inertial sensors and to use a physical vehicle model [START_REF] Gadola | Development and validation of a Kalman filter-based model for vehicle slip angle estimation[END_REF][START_REF] Li | Reliable vehicle sideslip angle fusion estimation using low-cost sensors[END_REF].

a) Integration: According to kinematic motion of vehicle, the derivative of sideslip angle is a function of lateral acceleration and vehicle velocities, as shown in Equation [START_REF] Jiang | Adaptive estimation of vehicle dynamics through RLS and Kalman filter approaches[END_REF]. The advantage of the direct integration method is based on the independence of tire properties, road friction and other vehicle parameters. Nevertheless, the sensor error will be accumulated in the integration process especially for the lateral acceleration signal, which is easily influenced by some internal and external factors.

βcdg = a y v x -ψ (1)
where a y is the lateral acceleration, V x is the longitudinal velocity, ψ is the yaw rate and β cdg is the sideslip angle of vehicle at center of gravity (COG).

b) Physical model: The physical model based approach is to estimate the sideslip angle according to the physical vehicle dynamics models. The role of the vehicle dynamics model is to extract the useful information (such as the sideslip angle) from the measurements. The accuracy of this approach is greatly affected by the quality of the dynamics models. In the application, the commonly-used dynamics model is usually a greatly simplified description of the complex vehicle behaviors. The variation of road condition and other vehicle parameters is not taken into account. As a consequent, the model errors is hardly avoidable. A typical vehicle dynamics model is the bicycle model, illustrated in Figure 1.

In order to circumvent the shortcomings of direct integration method and the model-based method, it is necessary to employ sensor fusion. In this way, the dynamics state of vehicle is estimated by the combination of physical model and integration method. Furthermore, during the sensor fusion process, different vehicle dynamics models and sensors can be combined to provide a robust estimation. As published in [START_REF] Klomp | Longitudinal velocity and road slope estimation in hybrid electric vehicles employing early detection of excessive wheel slip[END_REF], the longitudinal speed is estimated by combining wheel speed sensors and wheel torque sensors. The sensor fusion process can be generally concluded as a process to attribute different weightings or confidence to the measurement of each sensor. The Kalman filter is the most widely-used algorithm to realize the sensor fusion for estimation of vehicle dynamics states. Extended Kalman filter is implemented in [START_REF] Chen | Sideslip angle estimation using extended Kalman filter[END_REF] to estimate tire/road forces. Dual Kalman filter is proposed in [START_REF] Wenzel | Dual extended Kalman filter for vehicle state and parameter estimation[END_REF] for the estimation of side slip angle and steering stiffness. Particle filter methodology is employed for the sideslip angle in [START_REF] Wang | Real-Time Experimental Validation of Nonlinear Observer for Vehicle Dynamics Parameters Estimation : A Laboratory Vehicle Description[END_REF].

Until today, developing an accurate estimators for vehicle dynamics states is still a very challenging research topic. That is due to two aspects: 1) the need of an accurate and also simple dynamics model; 2) the need of a robust and time-efficient data processing technique to eliminate unpredictable sensor errors. The main contribution of this paper is to improve the performance of our estimator in both of the two aspects. For the aspect of vehicle modeling, we employed a non-linear tire model and the double track vehicle model, illustrated in Figure 1. Moreover, we have taken into account the vehicle's pitchroll movements and road inclination (both road bank angle and slope angle). For the other aspect, we proposed a best-wheel selection algorithm and combined with the Extended Kalman filter to facilitate the sensor fusion process. In general, we proposed a robust algorithm for real-time estimation of vehicle dynamics states in presence of sensor errors, road angles and tire slips. To evaluate the performance of our estimators, the proposed algorithm is implemented in our experimental vehicle DYNA, a Peugeot 308 [START_REF] Wang | Real-Time Experimental Validation of Nonlinear Observer for Vehicle Dynamics Parameters Estimation : A Laboratory Vehicle Description[END_REF], as shown in Figure 4.

This paper is organized as follows. Section 2 presents the construction of our proposed estimator. Section 3 introduces the experimental vehicles. Then, experiment results are described in Section 4. Finally, concluding remarks and future perspectives are given in Section 5.

II. ROBUST ESTIMATION OF VEHICLE'S LATERAL DYNAMICS STATES

In the literature, the estimator of dynamics states is developed based on the assumption that the inertial sensors are always available. In our context, the sensor failure may occurs. Therefore, the value of acceleration, yaw velocity and linear speed should be estimated with models rather than obtained directly by the measurement. For the simplicity and efficiency, the vehicle dynamics model is divided into 5 subsystems, represented by the 5 blocks in Figure 2. For each of the subsystem, we have developed a Kalman filter to optimize the estimation. The 5 blocks are the pitch-roll angle estimator, the accelerations estimator, the road angle estimator, the velocity estimator and the sideslip angle estimator. The advantage of using 5 Kalman filters is to make the observer more robust, even in presence of sensor errors. The estimation of pitch-roll angle, road angle and velocities are presented in our previous work [START_REF] Jiang | Adaptive estimation of vehicle dynamics through RLS and Kalman filter approaches[END_REF]. In this article, we focus on the robust estimation of lateral dynamics, especially the side slip angles and lateral tire forces. To clearly present the modification we proposed, we concluded our work as two parts: the best wheel selection algorithm and the EKF based sensor fusion algorithm.

A. Best-wheel selection

To begin with, we remind that the sideslip angle β cog is defined by the Equation [START_REF] Zhu | Real-time yaw rate prediction based on a non-linear model and feedback compensation for vehicle dynamics control[END_REF]. Then according to this definition and the double track model, the tire slip angle at each wheel (α f l , α f r , α rl , α rr ) can be calculated by Equation [START_REF] Klomp | Longitudinal velocity and road slope estimation in hybrid electric vehicles employing early detection of excessive wheel slip[END_REF].

β cog = arctan( v y v x ) ( 2
)
α f l = -δ l + β cog + L1 ψ vx α f r = -δ r + β cog + L1 ψ vx α rl = β cog -L2 ψ vx α rr = β cog -L2 ψ vx (3)
where v x , v y are the lateral and longitudinal speed at COG, δ l , δ r are the steering angle of front left and right wheel respectively, the index f l refers to the front left wheel.

Direct measurement of β cog or v y is not available in the ordinary cars. However, the tire slip angle (α f l , α f r , α rl , α rr ) could be estimated with tire models, as introduced in [5] [START_REF] Jihan | Vehicle sideslip and roll parameter estimation using GPS[END_REF]. Thus the problem of estimation of sideslip angle at COG ( β cog ) is transferred into the estimation of tire slip angle at each wheel. Theoretically, all the four tire slip angles can be computed with the tire models, but when the tire slip angle is large, the tire model become very non-linear and less accurate. Therefore, it is better to choose a best-wheel among the four wheels to compute the sideslip angle at COG β cog . The best wheel in our context means the wheel with the smallest lateral side slip. The mathematical criterion of selecting the bestwheel is represented by Equation ( 4)

β best =    δ -L1 ψ vx + α f L2 ψ vx + α r if |δ| -L ψ vx < -ε if |δ| -L ψ vx > ε (4)
where β best is the sideslip angle at COG estimated by the best wheel approach, α f is the slip angles of front tires , α r is the slip angle of rear tires, δ is the average steering angle, ε is a constant parameter.

When the front wheels are laterally slipping, the vehicle will be under-steering, δ -L ψ vx > 0. As a result, the sideslip angle at COG should be computed with rear wheel slip angles. Otherwise, when the rear wheels are laterally slipping, the vehicle will be over-steering, δ -L ψ vx < 0 , sideslip angle should be calculated by front wheels. The best wheel approach could eliminate the estimation errors caused by tire slip, due to that the slipping wheel is not used in the estimation. Only the wheel under small slip is chosen for further estimation.

B. Extended Kalman filter based sensor fusion

Even though the tire slip angle of the best-wheel is small, it is still challenging to accurately calculate its value, as many factors could influence the tire performance, such as temperature or humidity. In the literature, the linear tire model is employed. In order to take into account the non-linearity and the impact of vertical load, we chose the non-linear tire model, Dugoff model, as shown in Equation ( 5) and illustrated by Figure 3.

F y = C β β 1+β f (λ) f (λ) = λ(2 -λ) ; λ 1 λ ; λ > 1 (5)
where λ indicates the tire slip situation, for more details, the reader may refer to [START_REF] Wang | Real-Time Experimental Validation of Nonlinear Observer for Vehicle Dynamics Parameters Estimation : A Laboratory Vehicle Description[END_REF]. The Figure 3 also explains why we should chose the small slip wheel. When the excessive slip occurs, the tire force can hardly reflect the value of slip angle.

The static tire lateral forces are obtained based on the yaw dynamics and bicycle model, as expressed by Equation [START_REF] Li | Reliable vehicle sideslip angle fusion estimation using low-cost sensors[END_REF].

  ψI z m v a y m v a x   =   L 1 sin δ L 1 cos δ -L 2 sin δ cos δ 1 cos δ -sin δ 0     F xf F yf F yr   +cov(noise) (6)
However, the tire lateral forces are not generated instantaneously. Thus the transient tire model is employed to further improve the accuracy of estimation, depicted by Equation [START_REF] Ouahi | Simultaneous observation of the wheels' torques and the vehicle dynamic state[END_REF].

  Ḟ xf Ḟ yf Ḟ yr   =   0 0 0 0 p 1 0 0 0 p 2     F xf F yf -F yf F yr -F yr   + cov(noise) (7)
where F y f and F y r are static lateral tire forces calculated by Dugoff models, p 1 , p 2 are estimator parameters and configured by experiments.

Substituting the Equation (5, 6, 7) into (4) generates the sideslip angle at CDG, which is noted as β best . The advantage of best-wheel approach is the ability of estimating sideslip angle without integration. However, the estimation is based on simplified tire models, the model errors should not be ignored. Therefore, we employed the EKF algorithm to combine all these estimation methods and models as a sensor fusion. The Kalman filter for estimating sideslip angle can be written as:

βcdg vy = 0 0 0 0 β cdg v y + 1 0 ( ay-g sin ϕr vx -ψ)
+cov(noise)

β best 0 = 1 0 v x -1 β cdg v y + cov(noise) (8)
where a y , ϕ r , v x , ψ, ϕ r are estimated in former blocks, and therefore used as constant values in this block.

III. EXPERIMENTAL VEHICLE

A. Vehicle description

In this section we present the sensors and onboard acquisition system architecture of our experimental vehicle. The experimental vehicle DYNA, a Peugeot 308sw as shown in Figure 4, is instrumented by the laboratory HEUDIASYC UMR 7253 CNRS at Compiègne, France. This experimental platform is dedicated to validate the algorithm of estimating vehicle dynamics and the embedded real-time systems. All experimental data in this article is acquired by this system. The sensors we used could be generally classified as two categories, sensors used for validation and sensors used for input of the estimator. Kistler RoaDyn S625 wheel force transducers: It's able to measure all the tire-road contact forces and wheel torques in three dimensions. Four wheel force sensors are fixed at each wheel. These sensors are very expensive for ordinary cars. It is noted that they are used only for reference.

A scenario record camera is used to register the vehicle trajectory.

Sensors used for input measures: CROSSBOWVG700AB: It combines MEMSIC's high performance fiber optic gyros with silicon micro-machined (MEMS) accelerometer technology. It could provide a highly accurate measurement of Vertical Gyro (VG) and Inertial parameters.

CORRSYS-DATRON HT500: it is a non contact distance sensor. It provides measurement of the deflection between chassis and ground. They are installed respectively at four corners of vehicle body.

Available data on CAN bus: wheel rotation velocity, engine speed, yaw rate, brake pressure, lateral acceleration from the ESP, steering wheel angle. And a low-cost GPS.

These devices are located in the trunk of the car with the electric circuit system as shown in Figure 5. The monitoring equipment is located on the back left seat. The operator can manipulate the configuration of the system with a monitor and keyboard. Sensor data should be sent to a computer which has installed the estimation algorithm. However, some of these sensors cannot be directly connected to the computer. Therefore we developed an acquisition system based on the UEI PowerDNA Ethernet DAQ Cube to fulfill the task. It is capable of acquiring 48 analog channels using a 24-bit converter [START_REF] Wang | Real-Time Experimental Validation of Nonlinear Observer for Vehicle Dynamics Parameters Estimation : A Laboratory Vehicle Description[END_REF]. It consists of a computer UEI powered by a 200 MHz PowerPC processor running a real-time operating system Xenomai. The Xenomai is a real-time development framework cooperating with the Linux kernel. The software provides data acquisition by using the analog cards AI-217, shown in Figure 6. Digital filters with different cut-off frequency are installed in the cube to cut off high frequency noise. The acquisition frequency of the cube is 4 kHz, while the data is sent back in 200 Hz.

B. Software modules

This car is equipped with an industrial PC in the trunk. A particular prototyping environment, PACPUS 1 , is developed in C/C++ for the real time estimation system. Following the principle of component-oriented design, PACPUS provides users more versatility in their developments. Different components of PACPUS are developed in our work to fulfill different tasks. For instance, we develop the component "CubeClient" for the task of communication between the cube and the PC. The component is devoted to the task of sending the request, receiving the UDP packet and then decoding the UDP packet. Then the data is sent to the component of vehicle dynamics estimation. Similarly, the components for management of the other sensors are also constructed in our work. The estimation algorithm is developed in the .DLL form as a real-time application. The framework PACPUS employs the Qt API for graphical interfaces and can be integrated with other development environment for multi-sensor fusion. The schema of PACPUS is shown in Figure 7. Figure 9 shows the result of estimated road angles. The slalom behavior introduced disturbance into the road angle estimation. By using only inertial sensors, the estimation with an error of 5 degree is acceptable. Further improvement can be achieved by fusion with a precise 3D map. Figure 10 shows the estimation of yaw rate and heading angle in presence of GPS errors. The estimated lateral dynamic states are illustrated in Figure 11. The average estimation error is about 3%.

B. Experimental results of slalom test at banked track

Another test is conducted at banked high speed track. In the testing circuit, there are three parallel tracks: high, middle and low track. According to the information provided by UTAC CERAM, we know the bank angles at high, middle and low track are about 40 degree, 30 degree and 15 degree respectively. In order to hold stable at banked track, the vehicle speed reached 135 km/h, as shown in Figure 12. A slalom behavior is implemented by continuously changing from high track to low track. The estimation results of road angles are shown by Figure 13. When t=80s-90s, the vehicle keeps stable at low track, the estimation result is satisfactory. When the vehicle begins to do the slalom behavior, road angle estimation is disturbed but still correctly indicate the variation of road geometry. The lateral tire forces and sideslip angle are illustrated in Figure 14. The average estimation error is about 9%.

V. CONCLUSIONS AND PROSPECTS

This paper has presented a new algorithm to estimate vehicle dynamics states, in presence of road angles and sensor errors.

To reduce the dynamics model errors, we employed the non linear tire model and modified the double track model by taking into account the road inclination. We also developed an algorithm to detect excessive wheel slip and select the best wheel for estimation of side slip angle. In this way, the estimator stay robust even one wheel is undergoing excessive slip. Furthermore, we employed the extended Kalman filter to combine the estimation results of different approaches, which can be regarded as a process of sensor fusion. Experimental results are presented to evaluate the performance of the new observer.

In the future study, we will focus on improving estimation accuracy with a precise digital map.

Figure 1 .

 1 Figure 1. Bicycle Model and Double track model

Figure 2 .

 2 Figure 2. Algorithm overview

Figure 3 .

 3 Figure 3. Dugoff lateral force model characteristics

Figure 4 .

 4 Figure 4. Experimental vehicle DYNA, equipped with different sensors

Figure 5 .

 5 Figure 5. Embedded electronics and computer in the trunk

Figure 6 .

 6 Figure 6. The structure of the acquisition system

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Software architecture of acquisition system and estimation system

Figure 10 .Figure 11 .

 1011 Figure 10. Estimation of heading angle and yaw rate in presence of GPS errors

Figure 12 .Figure 13 .

 1213 Figure 12. Description of Test 2: slalom test at banked track

Figure 14 .

 14 Figure 14. Estimation of lateral dynamics at banked track: lateral forces of front wheels, lateral forces of rear wheels and side slip angle

The authors are with Sorbonne universités, Université de Technologie de Compiègne, Laboratoire Heudiasyc, CNRS UMR 7253, funded in part by the Region Picardie, VERVE project, MS2T LABEX and ROBOTEX. kun.jiang@hds.utc.fr, acorreau@hds.utc.fr, acharara@hds.utc.fr

The framework PACPUS is an open source with free license CeCILL-C. It is available at https://devel.hds.utc.fr/software/pacpus/wiki