Azadeh Moghtaderi

Patrick Flandrin

Pierre Borgnat

Trend Filtering via Empirical Mode Decompositions

Keywords: Empirical mode decomposition, trend filtering, adaptive data analysis, monthly mean carbon dioxide cycle, seasonality

The present work is concerned with the problem of extracting low-frequency trend from a given time series. To solve this problem, the authors develop a nonparametric technique called empirical mode decomposition (EMD) trend filtering. A key assumption is that the trend is representable as the sum of intrinsic mode functions produced by the EMD. Based on an empirical analysis of the EMD, the authors propose an automatic procedure for selecting the requisite intrinsic mode functions. To illustrate the effectiveness of the technique, the authors apply it to simulated time series containing different types of trend, as well as real-world data collected from an environmental study (atmospheric carbon dioxide levels at Mauna Loa Observatory) and from a large-scale bicycle rental service (rental numbers of Grand Lyon Vélo'v).

Introduction

Many real-world time series exhibit a "composite" behavior, in the sense that such a time series can be decomposed into a superposition of two "components." Typically one of these components can be classified as "trend," while the other component is classified as "fluctuation." (Note that the word "residual" is sometimes used instead of "fluctuation." In this paper, however, the term residual will be used in the context of empirical mode decomposition; see Section 2.) The problem of effecting such a decomposition, and classifying the resulting components as trend or fluctuation, is called the trend filtering problem (or trend estimation problem). Solving this problem is desirable, since an analysis of the trend component of a time series can often yield valuable information, e.g., for prediction. An obvious initial barrier to solving the trend filtering problem is that the terms "decomposition," "trend," and "fluctuation" are context-dependent. Indeed, given a time series generated by a particular physical system, it may be clear (based on physical intuition) how to solve the trend filtering problem. In the absence of physical intuition, it may still be possible to solve the trend filtering problem, provided one makes an ad hoc definition of trend; see [START_REF] Alexandrov | A review of some modern approaches to the problem of trend extraction[END_REF]. Such definitions may require extra assumptions concerning the nature of the time series.

A common ad hoc definition of trend is that of a "long-term change in the mean" [START_REF] Chatfield | The analysis of time series: An introduction[END_REF][START_REF] Alexandrov | A review of some modern approaches to the problem of trend extraction[END_REF]. This definition can lead to approaches which attempt to turn the trend filtering problem into one of regression. For example, it may be reasonable to assume that the time series has a trend component described by a low-degree polynomial. The coefficients of this polynomial can then be estimated by a standard polynomial regression; we again refer to [START_REF] Alexandrov | A review of some modern approaches to the problem of trend extraction[END_REF] for a more comprehensive discussion. Other approaches exist which do not impose such a strict model on the trend. For instance, nonparametric trend filtering assumes that the fluctuation possesses generic stationarity properties, and that the trend can be found by an ad hoc smoothing operation applied to the entire time series, e.g., using the Henderson filter [START_REF] Henderson | Note on graduation by adjusted average[END_REF] or the Hodrick-Prescott filter [START_REF] Hodrick | Postwar U.S. business cycles: An empirical investigation[END_REF]. Yet another possibility is to interpret the trend estimation problem in the frequency-domain sensefor instance, one can assume the trend is represented by a particular set of low-frequency (possibly polynomials or unit root) oscillations. This turns the trend filtering problem into a bona fide filtering problem. Viewed in this way, it may be profitable to use Wiener-Kolmogorov filtering [START_REF] Pollock | Wiener-Kolmogorov filtering frequency-selective filtering and polynomial regression[END_REF] to solve the trend filtering problem. Finally, it is worthwhile to mention that generalized "trend cycles," defined as a "short-term trend [that] generally includes cyclical fluctuations," have also been studied [START_REF] Alexandrov | A review of some modern approaches to the problem of trend extraction[END_REF]. Deciding if a trend cycle should be considered as trend (e.g., in any of the above senses) depends on the application and of course the observation scale.

In this paper, we introduce a novel approach to solving the trend filtering problem. We call this approach empirical mode decomposition trend filtering. It is philosophically similar to the "low-frequency approach" described in the preceding paragraph. Indeed, empirical mode decomposition trend filtering is based on the following definition: Trend is that component of a time series which is "slowly varying" in the sense that it is represented by the "slowest" intrinsic mode functions produced by the empirical mode decomposition (EMD). Recall [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] that the EMD is an algorithm which decomposes a time series into a finite additive superposition of "intrinsic mode functions," or IMFs. The IMFs are computed in an iterative fashion-each iteration produces an IMF which is "rapidly varying" relative to the residual time series. Thus our decomposition into components is effected by the EMD. The remaining question is "Which of the IMFs produced by the EMD should be deemed the slowest?" It is precisely this question which is addressed by EMD trend filtering. In particular, we attempt to answer this question by examining certain properties of the IMFs' energies and zero crossing numbers; these properties were first reported in (Flandrin et al., 2004b;[START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise, and Hurst exponent estimation[END_REF]. We give evidence which supports the fact that certain changes in these properties characterize the tipping point between trend and fluctuation.

It must be mentioned that the use of the EMD to solve the trend filtering problem has already been proposed in the literature. However, such work has either relied on an a priori model for the fluctuation (Flandrin et al., 2004a), or has considered the trend as being the final residual time series produced by the EMD [START_REF] Wu | On the trend, detrending, and variability of nonlinear and nonstationary time series[END_REF]. In a sense, using the EMD to solve the trend filtering problem shares common features with singular-spectrum analysis applied to the same problem [START_REF] Vautard | Singular-spectrum analysis in nonlinear dynamics with applications to paleoclimatic time series[END_REF][START_REF] Ghil | Interdecadal oscillations and the warming trend in global temperature time series[END_REF][START_REF] Vautard | Singular-spectrum analysis: A toolkit for short, noisy chaotic signals[END_REF]. This is because the SSA also effects a decomposition into oscillatory components. Like the EMD-based method proposed by [START_REF] Wu | On the trend, detrending, and variability of nonlinear and nonstationary time series[END_REF], a possible approach to solving the trend estimation problem using SSA is to identify the trend as the lowest-frequency oscillatory component. Other possibilities are to look for oscillatory components with prescribed smoothness or monotonicity properties; see [START_REF] Alexandrov | A review of some modern approaches to the problem of trend extraction[END_REF].

The rest of the paper is organized as follows. In Section 2, we briefly review some background material concerning the EMD. In Section 3 we state what trend means in the context of this paper. In Section 4, we describe the details of EMD trend filtering. The performance of EMD trend filtering is demonstrated in Sections 5, 6 and 7 through analyses of simulated and realworld time series. The concluding remarks are made in Section 8. Finally, the extension of the EMD trend filtering method to multiplicative models is provided in the Appendix.

The Empirical Mode Decomposition

The empirical mode decomposition (EMD) is an algorithm which decomposes a time series into a finite additive superposition of oscillatory components, each of which is called an intrinsic mode function (IMF); see [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. The EMD does not rely on any technical assumptions concerning the nature of the time series; note that this includes modelling assumptions. The basic idea is that IMFs are computed subject to two requirements: First, the number of local extrema and number of zero crossings of each IMF vary by at most one. Second, the mean of the upper and lower envelopes of each IMF should be identically equal to zero, where the envelopes are computed by means of a fixed interpolation scheme. (In the numerical results presented in this paper, we have confined ourselves to the use of cubic spline interpolation.) The IMFs are computed by means of an iterative scheme. This scheme however depends on a stopping criterion which guarantees that the requirements above are satisfied within a given tolerance while at the same time each extracted IMF is meaningful in both its amplitude and frequency modulations; we again refer to [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] for details.

To make this intuitive description more precise, let X = {X t } t≥0 be a (real, discrete-time, stochastic) process, and let X = (X 0 , X 1 , . . . , X N -1 ) be a realization of X. (These assumptions illustrate a notational convention that is used throughout the rest of the paper, namely that time series of length N are written in bold typeface and are regarded as elements of the Euclidean space R N .) As an initialization step, set i = 1 and ρ 0 = X . The EMD computes the IMFs of X using the following algorithm.

(1) Identify the local maxima and local minima of ρ i-1 .

(2) Together with the chosen interpolation scheme, use the maxima and minima from step (1) to compute the upper and lower envelopes of ρ i-1 . (3) Determine the local trend, denoted Q i , as the mean of the upper and lower envelopes from step (2). (4) Compute the local fluctuation, denoted h = X -Q i .

(5) If h is not an IMF, in the sense that it does not satisfy the two requirements described in the beginning of this section, then increment i and go to step (1) with ρ i-1 = h. [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] call this the sifting process; it is this process which depends on the stopping criterion.) (6) If h is an IMF, in the sense that it satisfies the two requirements described at the beginning of this section, then the ith intrinsic mode function of X is M i = h, and the ith residual is ρ i = X -M i . Increment i, then go to step (1).

The algorithm halts when the ith residual has no further oscillations, in the sense that it has no local maxima or local minima. We denote by I the largest index for which M i is defined. Then

X = I i=1 M i + ρ I .
(1)

In this decomposition, M 1 through M I can be thought of as containing a "spectrum" of local oscillations in X , with the shortest-period (highest frequency) oscillations represented in M 1 and the longest-period (lowest frequency) oscillations represented in M I . The computational complexity of the algorithm depends on X , the chosen interpolation scheme, and the stopping criterion. However, the algorithm usually halts in a reasonably small number of steps. For example, it is known (Flandrin et al., 2004a) that if X is a broadband process (a broadband process includes a relatively wide range (or band) of frequencies), then the decomposition produced by the EMD has an almost dyadic filter-bank structure, typically with I ≈ log 2 N . Moreover, it is known that the sifting process typically halts after some tens of iterations [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF].

Trend in EMD

As discussed in Section 1, the term "trend" is meaningless and has to be made more precise in order to be useful. In this section we state what we mean by trend in this paper and in the context of EMD trend filtering. To begin with, let us introduce some notation.

Let Y = (Y 0 , Y 1 , . . . , Y N -1 ) be a realization of a process Y = {Y t } t≥0 , and let C = (C 0 , C 1 , . . . , C N -1 ) ∈ R n be a trend component. Assume also that Y is a broadband process with a continuous spectrum. From Y and C we may form two new time series: The first is Y + C, the additive mix of Y and C; the second is CY, the multiplicative mix of Y and C. (Here the multiplication is being performed componentwise.) In either case, we say that Y is fluctuation of the mix. Now let X be the additive or multiplicative mix of Y and C. The question we wish to answer is: "Solely given X as data, under what conditions should it be possible to accurately estimate C from X ?" To do so, we must constrain the trend and fluctuation of the mix in some fashion. We take the following pragmatic approach that is based on properties of EMD.

Recall that in EMD, the successive IMFs are oscillations going from high frequency to low frequency, and that this property is valid locally in time (there is not necessarily a global separation of spectrum of successive IMFs) [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. A loose "definition" of a trend in this paper is that C is locally slowly varying as compared to Y. Hence, a pragmatic way of satisfying this is that the trend should be obtained as the sum of the last few IMFs and the residual extracted from X .

Let us now turn the attention to some properties of fluctuation of the mix which can also define (in contrast) the trend. First, and in agreement with Flandrin et al. (2004b) and [START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF], the mean frequency of the successive IMFs of broadband processes decrease, similarly to constant-Q filter-banks, with a factor near 2. This will be the first criterion studied in Section 4.1 by estimating the mean frequency from the number of zero crossings of IMFs. Second, the finding of [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise, and Hurst exponent estimation[END_REF] is that the "energy" of the IMFs of many broadband processes decreases as the index of the IMFs increases. This has been first reported and demonstrated in [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise, and Hurst exponent estimation[END_REF] for fractional Gaussian noise (fGn) processes [START_REF] Embrechts | Selfsimilar Processes[END_REF] which are convenient models for generic broadband processes. We will provide more discussion in Section 4.2 for the validity of these characteristics. An explicit assumption in our work is that the fluctuation Y contaminating the trend C have such energy profile. This does not exclude situations with a substantial energy increase downwards low frequencies, as is the case for fGns with Hurst exponent H > 1/2. Indeed the decreasing energy condition does not apply directly to the broadband processes, but to their IMFs. In practice, given the previously mentioned dyadic structure for the IMF spectra, processes Y with power spectra diverging as f -α at the zero frequency are admissible provided that α < 1.

In the presence of a trend, the prescription used in this paper is that the IMF index which shows a rupture in the two properties described above separates the trend from the fluctuation. It follows from this prescription that a trend in the present work is neither restricted to be monotonic nor to be some polynomial functions. The trend in this work can however contain oscillations while in [START_REF] Wu | On the trend, detrending, and variability of nonlinear and nonstationary time series[END_REF] only the residual of EMD was deemed a trend, hence constraining it to have no oscillations at all.

In the following two sections, we will describe in details the properties discussed above and their abilities in separating trend and fluctuation.

EMD Trend Filtering

Let X be the additive mix of Y and C, where these entities are given as in the previous section. As described there, our goal is to accurately estimate C from X . This section is devoted to describing EMD trend filtering which can be used to obtain such an estimate.

The following notation and terminology will be employed throughout this section. Let M i be the IMFs of X , where 1 ≤ i ≤ I, and let i * be such that

C i * = I i=i * M i + ρ I (2)
is the best approximation to C in the Euclidean metric. We call i * the best index and C i * the best approximation of C. Estimating C is equivalent to estimating the best index. If i * is an estimate of i * , then we denote by C i * the corresponding estimate of C. EMD trend filtering, described over the course of the next three subsections, actually consists of three approaches to estimating i * . These are called, respectively, the ratio, energy, and energy-ratio approaches.

The extension of the trend filtering method to multiplicative mixes including all simulations can be found in the Appendix.

Ratio approach

In this subsection we describe the first approach to estimate i * , which is based on an empirical property of the zero crossing numbers of IMFs.

Let us first establish some additional notation. For a given time series, the zero crossing number of its ith IMF is denoted by Z i , and let us define

R i = Z i-1 /Z i for i ≥ 2. (This is well-defined since Z i ≥ 1; see Section 2.)
Of course, R i depends fundamentally on the given time series; since the particular time series is always clear from context, we suppress this dependence. We call R i the ith ratio of the zero crossing numbers (ith RZCN). It has been observed by Flandrin et al. (2004b) and [START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF] that if the time series under study is a realization of a generic broadband process, the approximation R i ≈ 2 holds.

Let us first support this observation. We considered 20 broadband processes of the following types: 17 fGn processes with H = 0.1, 0.15, 0.2, . . . , 0.9, two stationary AR(2) processes, and a nonstationary AR(2) process with time-dependent coefficients. For each process in the collection, we simulated B = 10000 realizations of length N = 2000, then computed the IMFs of each realization along with their zero crossing numbers. Denoting the ith RZCN of its bth realization by R i,b , where 2 ≤ i ≤ I b , and setting 1 displays this empirical distribution, and supports the contention that R i ≈ 2. In fact, this distribution is approximately Gaussian with mean 2. Furthermore, it is evident from Fig. 1 that apart from the expected peak at 2, we also observe several smaller but visible peaks at higher values. These peaks appear to be due to the presence of high-order IMFs; indeed, these slowly oscillating modes have small zero crossing numbers. Because RZCNs are calculated as the ratio of two integers, if the numerator is a small number, then the distribution of the elements of R will have peaks at integer or rational values such as 2, 5/2, 3, 4/3, etc. Hence, RZCNs with integer or rational values for small denominators have slightly higher expected probabilities than neighbouring values.

R b = (R 2,b R 3,b • • • R I b ,b ), we then computed the empirical distribution of the elements of R = ( R 1 R 2 • • • R B ). Fig.
Generically, the approximation R i ≈ 2 fails for i near the best index i * . This observation is supported by the following data. For each broadband process in the collection and using its realizations, we constructed 10000 additive mixes, using C 3 (displayed in Fig. 5) as a trend. We then computed the IMFs of each mix along with their RZCNs and set R b and R as described earlier. The left-hand plot in Fig. 2 displays the empirical distribution of the elements of R for additive mixes, and supports the contention that R i ≈ 2 fails. In fact this empirical distribution is non-Gaussian as its side peaks grows taller in comparison with the distribution shown in Fig. 1. The problem however is that it is not yet clear whether R i ≈ 2 fails around i * .

To clarify this, we proceeded with further simulations. For each broadband process in the collection, we used the IMFs obtained for each mix and used the knowledge of C 3 to evaluate the best index i * (see Section 5 for details.)

For each mix, we then computed the best approximation of the fluctuation by eliminating those IMFs whose indices are greater than or equal to i * (we call this detrending the mix.) We set R b and R for the remaining IMFs and computed the empirical distribution of the elements of R. This distribution, shown in the right-hand plot in Fig. 2, is Gaussian with mean 2 as was the case in Fig. 1. We therefore conclude that R i ≈ 2 fails around i * .

Based on what we described above, we propose to estimate i * by choosing i * to be the smallest index i for which R i is "significantly different from 2". We refer to this as the ratio approach. The results of our simulations for broadband processes suggest that in order to conclude whether or not R i is "significantly different from 2", a common threshold test can be used. For 0 ≤ p ≤ 100, we therefore compute p% and (100 -p)% significance level of the empirical distribution shown in Fig. 1 as the left threshold and the right threshold respectively. At the end, any RZCN outside of the appropriate right and left thresholds is considered significantly different from 2.

A weakness of the ratio approach is that, since selection of the left and right thresholds is based on empirical results, it is always possible that for a given p, the smallest i for which R i appears significantly different from 2 is a false detection. In Section 5, we will discuss how to select an optimum p.

Energy approach

In this subsection we describe the second approach to estimate i * , which is based on an empirical property of the so-called "energy" of the IMFs.

To describe this property, we need to establish some additional notation. Let {Z t } t≥0 be an arbitrary process. For a given time series which is a realization of {Z t }, we define the energy of its ith IMF, denoted G i , by

G i N -1 t=0 |M i t | 2 , 1 ≤ i ≤ I.
Assume now that we have B different time series obtained from {Z t }. Given the bth time series, 1

≤ b ≤ B, if G i,b
denotes the energy of its ith IMF, the averaged energy of its ith IMF is defined by G i broadband process in the collection, we then computed G i and observed that G i increases at the best index i * . The right-hand plot in Fig. 4 displays log 2 G i computed for detrended additive mixes only up to i = 5. This is because for some examples, i * > 5 but for the majority i * = 5. Based on the above discussion, identifying the smallest index i ≥ 2 such that G i > G i-1 evaluates i * . This approach is called the energy approach.

1 2 3 4 5 6 7 8 -4 -3 -2 -1 0 1 2 3 4 5 6 i log 2 G i 1 1.5 2 2.5 3 3.5 4 4.5 5 -4 -3 -2 -1 0 1 2 3 i log 2 G i
As for the ratio approach, one could think of looking for significant increases which would be based on some statistical information about the dispersion of energy of each IMF. This viewpoint has been considered first for white Gaussian noise in [START_REF] Huang | A confidence limit for the empirical mode decomposition and Hilbert spectral analysis[END_REF] and further generalized in Flandrin and Gonçalves ( 2004) and Flandrin et al. (2004a), even in a detrending perspective. The limitation however is that the associated confidence intervals depend strongly on some prior knowledge about the spectra of broadband processes. This is the main reason that we do not follow such direction, as we are interested in a procedure which is not model-dependent.

A limitation with the energy approach is that one is often given a single time series to use for trend estimation. Computation of energy based on only one time series may cause an increase in G i when i = i * .

Energy-ratio approach

To overcome limitations of the previous approaches, we introduce the last and most important approach to estimate i * . As described, the energy and ratio approaches are confronted with possible false detections of the smallest index which does not associate with the trend. Since the criteria proposed by the energy and ratio approaches to evaluate i * are independent, the number of false detections can be reduced by combining these two approaches.

To be more precise, for each 2 ≤ i ≤ I, we compute each index i such that G i > G i-1 . For a fixed p, we also evaluate every index i where R i is significantly different from 2. We then evaluate i * to be the smallest common index in both approaches. This approach is called the energy-ratio approach.

Performance Evaluation of the EMD Trend Filtering; Evaluation of an optimum p

We follow two main goals in this section. The first goal is to evaluate the overall performance of the EMD trend filtering. The second goal is to empirically evaluate an optimum p which can improve the performance of the energy-ratio approach in comparison with the energy and ratio approaches. In order to do so, we use 10 simulated examples including 6 additive and 4 multiplicative mixes such that

X k =            C k + Y k , 1 ≤ k ≤ 3 C k Y k , 4 ≤ k ≤ 6 (C k-3 -1) + Y k-3 , k = 7 C k-3 + Y k-3 , 8 ≤ k ≤ 9 (C k-9 + 1)Y k-9 , k = 10.
(3)

In order to construct the above mixes, we use the following.

Let Y k = {Y k t } t≥0 , 1 ≤ k ≤ 6, be 6 generic broadband processes such that for 1 ≤ k ≤ 2, we have where {ζ t } and {ξ t } are two independent white noise processes with variance 10 4 , and for 3 ≤ k ≤ 6, we have 4 fGn processes with Hurst exponents 0.7, 0.5, 0.15, and 0.75 respectively. Let

Y 1 t = 0.8Y 1 t-1 -0.4Y 1 t-2 + ζ t , and
Y 2 t = 0.2Y 2 t-1 + 0.5Y 2 t-2 + ξ t , 500
Y k = (Y k 0 , Y k 1 , . . . , Y k N -1 ) be a realization of Y k . Now, let us assume that C k = (C k 0 , C k 1 , . . . , C k N -1
), 1 ≤ k ≤ 6, are 6 trends where for 1 ≤ k ≤ 4, we have 4 randomly constructed trends using peacewise linear and cubic spline techniques and for 5 ≤ k ≤ 6, we have

C 5 t = 2 -e -(t-1000) 2 2×400 2
, and C 6 t = 1.5 + cos(2πf s t), f s = 0.002. 3). We denote the bth realization of the kth example by b k . In order to achieve the goals described earlier in this section, we started with the following computations. We applied EMD to X b k (or log |X b k | for multiplicative mixes) in order to extract its IMFs. Denote I, M i , and ρ I by I b k , M i b k , and ρ I b k respectively.

For each i † ∈ {1, 2, . . . , I b k }, we computed

C b k i † = I b k i=i † M i b k + ρ I b k , (4)
and the Euclidean distance (ED) between We then estimated C k using three different trend filtering methods. The methods we used are the Hodrick-Prescott (HP) filter [START_REF] Hodrick | Postwar U.S. business cycles: An empirical investigation[END_REF], the singular-spectrum analysis (SSA) [START_REF] Vautard | Singular-spectrum analysis: A toolkit for short, noisy chaotic signals[END_REF] and the EMD trend filtering using ratio, energy and energy-ratio approaches. We denoted the trend estimates obtained above by C b k m where the letter m indicates the type of trend filtering. For simplicity, we selected m to be r, g, and gr to refer to the ratio, energy, and energy-ratio approaches respectively. Since the ratio and energy-ratio approaches are dependent on p, we denoted the estimates for these methods by C 1,2, 3, and 4 report all the averaged EDs and ENs computed using these parameters.

C k and C b k i † , denoted E b k i † . The best index i * is that i † which results in minimum E b k i † , denoted E b k i * . Clearly, C b k i * is the best approximation of C k . Let Y b k i * = X b k -C b k i * . Here Y b k i * is
To evaluate the performance of the EMD trend filtering which was the first goal in this section, we make two attempts. The first attempt is to compare the best approximation of C k obtained from the EMD trend filtering with estimates obtained from the HP filter and the SSA. In order to do so, for each k, we compared the reported EDs from the second column of Table 1 with those from the third to sixth columns. Since these EDs are comparable, we conclude that the EMD trend filtering performs similarly to the HP filter and the SSA. Note that since both HP filter and the SSA are dependent on free parameters, the quality of their performance can vary in comparison with the EMD trend filtering. This is clear from the reported EDs in Table 1.

k E k i * E k hp E k hp E k ssa E k ssa E Y k i * E Y k 1 0.
The second but also necessary attempt we make is to compare the fluctuation of each mix with the best approximation of the fluctuation. This is done by comparing the averaged ENs reported in the last two columns of Table 1. The fact that these two columns are comparable is an indication that the EMD trend filtering is a well-performed method in estimating the trend.

Recall the second goal in this section which is to empirically evaluate an optimum p which makes the energy-ratio approach to perform better than the energy and ratio approaches. We should note that by using the term optimum here, we really mean within the given examples.

In order to obtain such p, we used the averaged EDs reported in Tables 2, 3 and 4. Comparing the values reported in Tables 2 and3 shows that for majority of p and k, the averaged EDs associated with the energy-ratio approach are smaller than those for the ratio approach. This means that the energyratio approach performs better than the ratio approach regardless of p. As a result, selection of p should only depend on how the energy-ratio approach compares with the energy approach. We therefore compare the averaged EDs reported in Table 3 350 26.50 26.97 27.76 3.111 7.195 8.149 27.58 32 3.012 2.246 1.436 27.90 28.47 29.30 3.317 7.422 8.594 29.15 34 3.201 2.295 1.511 29.11 29.75 30.62 3.518 7.643 8.965 30.50 35 3.300 2.315 1.546 29.71 30.46 31.37 3.624 7.736 9.158 31.27 40 3.694 2.384 1.695 32.73 33.70 34.41 4.145 8.220 9.892 34.54 45 3.979 2.426 1.794 34.23 35.19 35.81 4.573 8.564 10.34 36.21 it reaches its minimum (denoted p k 2 and marked with a star in Table 3) and increases again but it does not exceed E k g (at least for maximum p = 45.) This observation indicates that first of all, an optimum p is not unique as it depends strongly on the type of example. Second, there is a wide range of p values which make the energy-ratio approach to perform better that the energy approach. We therefore select an optimum p, denoted p * to be such that p * > max k p k 1 and also p * < max k p k 2 . We therefore can select any 16 < p * < 26. In this paper, we use p * = 18.

Simulated Examples

In this section, we demonstrate the performance of the energy-ratio approach in estimating i * via two simulated examples. The examples use here are the additive mix X 1 and the multiplicative mix X 5 introduced in Section 5 for further analysis.

The notation used in this section is exactly the same as in Section 5 except that since we only work with one time series of each mix, we replace b k in the notation with k. For the ratio and energy-ratio approaches, we use p * = 18.

Simulated example 1

Recall Y 1 and C 1 from Section 5. Let

Y 1 = {Y 1 0 , Y 1 1 , . . . , Y 1 N -1 } be a realization of Y 1 . Set the additive mix X 1 = Y 1 + C 1 for N = 2000.
We apply EMD to X 1 and extract its IMFs and obtain I = 10.

Using the IMFs obtained for X 1 , we first compute C 1 i † for 1 ≤ i † ≤ 10 as in Eq. (4). We then compute the EDs between C 1 and C 1 i † , denoted E 1 i † , and the EDs between X 1 and C 1 i † , denoted E Y 1 i † . These are reported in the first two rows of Table 5. Based on these reported values, we can see that since i † = 8 results in minimum E 1 i † , we conclude that i * = 8. An additional support for this selection is that E Y 1 8 is the closest value to the EN of Y 1 which is 5.493. We now want to compare the performance of the ratio, energy and energy-ratio approaches in estimating i * .

Looking at the energy of the IMFs of X 1 , we observe that the IMF indices for which G i > G i-1 are i = {6, 8, 9, 10}. Based on the energy approach, we evaluate i * = 6 which is the smallest observed index in this case. Looking at the RZCN of each IMF on the other hand, we observe that the IMF indices for which R i is significantly different from 2 are i = {4, 5, 7, 8, 9, 10}. Based on the ratio approach, we evaluate i * = 4. Finally, the energy-ratio approach and2: The first two rows are associated with X 1 , where 1 ≤ i † ≤ 10. The last two rows are associated with X 5 , where 1 ≤ i † ≤ 11. evaluates i * = 8 as the smallest common IMF index between the energy and ratio approaches.

It is clear from above that the energy-ratio approach has performed excellently in estimating i * by eliminating the false detections in the ratio and energy approaches. Fig. 6 displays the energy and ratio approaches together with the estimated trend using i * .

Simulated example 2

Recall Y 5 and C 5 from Section 5. Let Y 5 = {Y 5 0 , Y 5 1 , . . . , Y 5 N -1 } be a realization of Y 5 . Set the multiplicative mix X 5 = C 5 Y 5 for N = 2000. We apply EMD to log |X 5 | and extract its IMFs and obtain I = 11. Similarly to the previous example, we use the IMFs obtained for log |X 5 | to first compute log C 5 i † for 1 ≤ i † ≤ 11. We then compute the EDs between log C 5 and log C 5 i † , denoted E 5 i † , and the EDs between log |X 5 | and log C 5 i † , denoted E Y 5 i † . These are reported in the last two rows of Table 5. Based on these reported values, we can see that since i † = 10 results in minimum E 5 i † , we conclude that i * = 10. An additional support for this selection is that E Y 5 10 is the closest to the EN of log |Y 5 | which is 47.38. We now compare the performance of the ratio, energy and energy-ratio approaches in obtaining i * .

Looking at the energy of the IMFs of log |X 5 |, we observe that i = {7, 10, 11}. Based on the energy approach, we evaluate i * = 7. Looking at the RZCN of each IMF on the other hand, we observe that i = {2, 3, 8, 9, 10, 11}. Based on the ratio approach, we evaluate i * = 2. Finally, the energy-ratio approach evaluates i * = 10 as the smallest common IMF index between the energy and ratio approaches.

It is clear from above that the energy-ratio approach has performed excellently in estimating i * by eliminating the false detections in the ratio and energy approaches. Fig. 7 displays the energy and ratio approaches together with the estimated log-trend using i * .

Real-World Examples

In this section we demonstrate the performance of the energy-ratio approach via two real-world examples. The first example is the monthly mean carbon dioxide (CO 2 ) data from Mauna Loa and the second example is the Grand Lyon-Vélo'v bicycle rental data from the city of Lyon in France.

Monthly mean CO 2 at Mauna Loa

In this section, we analyze the monthly mean CO 2 data collected from March 1958 to March 2010 and measured at Mauna Loa observatory in Hawaii (Available via FTP:ftp://ftp.cmdl.noaa.gov/ccg/co2/trends /co2 mlo.txt. The authors have received permission from Dr. Pieter Tans in order to use this data.) The left-hand plot in Fig. 8 displays the monthly mean CO 2 data at Mauna Loa. After removing the averaged seasonal cycle expected in this data, a trend is obtained. This trend is given at the URL together with the data, and it will serve as a reference for a comparison with the result from EMD trend filtering. For more information on the known seasonal cycle and trend calculation see the URL provided above. The right-hand plot in Fig. 8 displays the one year cycles of the monthly mean CO 2 data after removing the expected trend together with their average. We call this average the expected annual cycle.

We now use EMD trend filtering for monthly mean CO 2 data in order to estimate its underlying trend. Applying EMD to this data, we obtain I = 3 and following the energy-ratio approach, we evaluate i * = 3. The Left: Expected trend together with the estimated trend using EMD trend filtering. Since these two trends look very similar, the smaller plot is made to display only a small portion of these trends. Right: Yearly cycles of the detrended monthly mean CO 2 data using the estimated trend together with their average (dark black line). The dashed line displays the difference between the expected and estimated annual cycles.

left-hand plot in Fig. 9 displays the estimated trend plotted together with the expected trend obtained from removing the seasonal cycle. Since these two trends look very similar, the smaller plot is made to display only a small portion of these trends. It is clear that the estimated trend from the EMD trend filtering is only a smoother version of the expected trend.

After subtracting the estimated trend from the data, we divide the detrended data into one year cycles and then average over all cycles to obtain the estimated annual cycle. The right-hand plot in Fig. 9 displays all the one year cycles of the monthly mean CO 2 data after removing the estimated trend together with the estimated annual cycle. The dashed line in Fig. 9 displays the difference between the expected and estimated annual cycles. This difference confirms the strong similarities between the two annual cycles.

Grand Lyon-Vélo'v

In this section we analyze the data from Vélo'v, the community shared bicycle program that started in Lyon in May 2005 (for more information, see http://www.velov.grandlyon.com.) The program Vélo'v is a major initiative in public transportation, in which bicycles are proposed to rental by anyone at fully automated stations in many places all over the city, to be returned at any other station. Such a community shared system offers both a new and versatile option of public transportation, and a way to look into the movements of people across the city. In order to understand the dynamics of this system, a question is to estimate and model the evolution in time of the number of rentals made throughout the city [START_REF] Borgnat | Studying Lyon's Vélo'v: A statistical cyclic model[END_REF]. The left-hand plot in Fig. 10 displays the raw data which is the number of hourly rentals for two years of activity of the Vélo'v system from December 2005 to December 2007. (The authors would like to thank JCDecaux for providing access to this data.)

The number of rentals contains cyclical patterns over the day (e.g., more activities during the day, mainly at specific rush hours, than during the night) and the week (e.g., more activities during week-days than week-ends). It also contains superimposed fluctuations due to external contingencies (e.g., rain or holidays) and a general multiplicative trend over the months [START_REF] Borgnat | Studying Lyon's Vélo'v: A statistical cyclic model[END_REF]. We apply EMD trend filtering to this data in order to estimate underlying multiplicative trend. We obtain I = 12 and using the energyratio approach we evaluate i * = 10 which we use to estimate the trend. The estimated trend is displayed in the left-hand plot in Fig. 10 where superimposed over the data.

This trend is meaningful for the data, and can be related to, and explained by, two effects: (i) the system was expanded in 2005 and 2006 at the same time it was already in exploitation, hence, there is a long-term increase of the hourly rentals over the two years of data, (ii) because of seasonal effects, the use of Vélo'v is smaller during winter, and also during the main summer holidays; this causes several drops of the trend, during winter and also summer holidays. Using the estimated trend, detrended Vélo'v data are obtained by dividing the number of hourly locations by the estimated trend. This is displayed in the right-hand plot in Fig. 10. The detrended data is, visually, more stationary than the raw data. This allows a good estimation of the cyclic pattern over the week of the number of hourly rentals. Fig. 11 displays the weekly cycles of the Vélo'v data after removing the estimated trend, and the average over all the weeks.

The estimate of the average usage of the Vélo'v bicycles as a function of time in the week, is meaningful in that it reveals the main features of the Vélo'v activity: during week-days, there are three sharp peaks of rentals in the morning, noon and the end of the afternoon; during the week-ends, there are small peaks at noon, and smooth and large peaks during the afternoon.

Finally, let us note that here the multiplicative trend estimation procedure was applied to a case where the underlying process that the trend multiplies to is not actually a broadband process: it is more specifically a periodic process (with clear periods of one week and one day) with added fluctuation. Nevertheless, the procedure is able to find the relevant multiplicative trend describing the evolutions at the scale of the seasons, and that is used to detrend the data. This is believed to be due mostly to the fact that fluctuations have typical periodic scales (one day or one week) which are much smaller than the typical scale of evolution (several months) of the trend, making of this scale separation a prerequisite that might be more important than the existence of a broadband spectrum in a stricter sense.

Conclusion

An automated method has been proposed to filter the trend in a time series, whose principle is to extract the lowest frequency intrinsic mode functions (IMFs) via empirical mode decomposition (EMD). The core of the method is to decide which IMFs belong to the trend, on the rationale that a trend causes both a departure of the ratio of zero crossing numbers from 2, and an increase of the energy contained in the low-frequency IMFs, as compared to the expected behavior of broadband processes. Combining both criteria, the procedure was shown to work well on several examples with additive or multiplicative trends. We emphasize that the approach is fully data-driven (as is EMD) and, besides the parameters of the decomposition itself, the EMD trend filtering method has only one free parameter which is the level of significance p.

Many numerical examples were reported to illustrate the robustness of the EMD trend filtering and its potential interests have been further illustrated on two real-world examples: the CO 2 data which displays an additive trend, and the Vélo'v data which shows a multiplicative trend. In both cases, filtering of the trends allows us to propose an estimation of the cycle inside the data (annual cycle for the CO 2 data, weekly and daily cycles for the Vélo'v data) that compares favorably to existing methods both for the extracted trends and estimated cycles. A strength of the method is that it works, even if the fluctuations above the trend do not follow exactly a priori behaviors for the fluctuation that where used to design empirically the test (displaying for instance oscillatory behaviors more than the assumed broadband behavior.) This is related to its character as a fully data-driven and model-free approach.

A perspective of this work would be to go beyond trend-filtering and use the same type of approach to group together IMFs obtained by EMD in several signals describing a trend, then the major cycles, and finally the rapid fluctuation. This would be an interesting asset for the model-free decomposition of processes.

Acknowledgement

The authors would like to thank the anonymous reviewers and Associate Editor for their helpful comments and suggestions. Most of the work reported here was completed during the postdoctoral stay of Azadeh Mogh- log C 3 to evaluate the best index i * . We then used i * to detrend each logtransformed mix and recomputed R b , G i,b , R, and G i for the remaining IMFs. The bottom left-hand plot in Fig. 13 displays the empirical distribution of the elements of R and the bottom right-hand plot displays log 2 G i both for detrended log-transformed multiplicative mixes. The comparison between the left-hand and the right-hand plots in Fig. 13 indicates that R i ≈ 2 fails and G i increases both at the best index i * .

These simulations altogether validate the fact that, after log-transformation, multiplicative mixes have properties similar to the additive mixes. Hence, the energy-ratio approach is expected to operate in a similar way to estimate the log-transformed trend, with the exception that the left and right thresholds used for the ratio approach are different. The appropriate left and right thresholds are now p% and (100 -p)% significance level of the empirical distribution shown in the left-hand plot in Fig. 12.

Figure 1 :

 1 Figure 1: Empirical distribution of the elements of R for broadband data: Computed for 10000 realizations of 20 broadband processes in the collection. Each line of different type and color associates with a broadband process in the collection.

Figure 2 :

 2 Figure 2: Empirical distribution of the elements of R for additive mixes: Left: Computed for additive mixes obtained by adding C 3 from Fig. 5 and realizations of broadband processes in the collection. Right: Computed for detrended additive mixes.

Figure 3 :

 3 Figure 3: log 2 G i : Computed for 10000 realizations of 20 broadband processes in the collection.

Figure 4 :

 4 Figure 4: log 2 G i : Left: Computed for additive mixes. Right: Computed for detrended additive mixes and displayed only up to i = 5.

Figure 5 :

 5 Figure 5: Trends used in simulated examples: C k for 1 ≤ k ≤ 6.

Fig. 5

 5 Fig. 5 displays C k for 1 ≤ k ≤ 6 when N = 2000.For each k, we created B = 10000 realizations of length N = 2000 of Y k and constructed the mixes for each realization following Eq. (3). We denote the bth realization of the kth example by b k . In order to achieve the goals described earlier in this section, we started with the following computations. We applied EMD to X b k (or log |X b k | for multiplicative mixes) in order to extract its IMFs. Denote I, M i , and ρ I by I b k , M i b k , and ρ I b k respectively.

  the best approximation of the fluctuation Y b k . We computed the Euclidean norm (EN) of Y b k i * and Y b k and denoted them by E Y b k i * and E Y b k respectively.

  b k ,p m . After all, we computed the ED between C k and each trend estimate and denoted them by E b k m (or E b k ,p m .) For each k, we then averaged all the ENs and EDs computed above over B realizations and denoted them by E k m (or E k,p m ). In this paper, in order to obtain E k hp , we used two free parameters of 10 5 and 5 × 10 5 , for E k ssa , we used the window lengths of 100 and 200, and for E k,p r and E k,p gr , we used 27 fixed p where 0 ≤ p ≤ 45. Tables

Figure 6 :

 6 Figure6: EMD trend filtering for simulated example 1: Top left: The energy approach. The small circles are log 2 G i for 1 ≤ i ≤ 10 and the small triangles mark those indices i ≥ 2 where G i > G i-1 . Bottom left: The ratio approach. The small circles are log 2 R i for 2 ≤ i ≤ 10 and the small triangles mark those indices i where R i is significantly different from 2. The dashed lines are the averaged left and right thresholds of the distribution shown in Fig.1when p * = 18. Top right: X 1 vs. C 1 gr . Bottom right: C 1 (dashed line) vs. C 1 gr (solid line).

Figure 7 :

 7 Figure 7: EMD trend filtering for simulated example 2: Top left: The energy approach. Bottom left: The ratio approach. The dashed lines are the averaged left and right thresholds of the distribution shown in the left-hand plot in Fig. 12 when p * = 18. Top right: log |X 5 | vs. log C 5 . Bottom right: log C 5 (dashed line) vs. log C 5 (solid line).

Figure 8 :

 8 Figure 8: Monthly mean CO 2 data and the expected annual cycle: Left: Monthly mean CO 2 data from March 1958 to March 2010. Right: Yearly cycles of the detrended data using the expected trend together with their average (dark black line).

Figure 9 :

 9 Figure9: Estimated trend and annual cycle for the monthly mean CO 2 data: Left: Expected trend together with the estimated trend using EMD trend filtering. Since these two trends look very similar, the smaller plot is made to display only a small portion of these trends. Right: Yearly cycles of the detrended monthly mean CO 2 data using the estimated trend together with their average (dark black line). The dashed line displays the difference between the expected and estimated annual cycles.

Figure 10 :

 10 Figure 10: Vélo'v raw and detrended data: Left: The raw Vélo'v data together with the estimated trend using EMD trend filtering. Right: Detrended Vélo'v data.

Figure 11 :

 11 Figure 11: Weekly detrended Vélo'v data: The weekly cycles of the detrended Vélo'v data and their average (dark black line).

Figure 13 :

 13 Figure 13: Empirical distribution of the elements of R and log 2 G i for logtransformed multiplicative mixes: Left: Computed for log-transformed multiplicative mixes obtained by multiplying 1 + C 3 from Fig. 5 and realizations of broadband processes in the collection. Right: Computed for detrended log-transformed multiplicative mixes.

Table 1 :

 1 Performance evaluation of the EMD trend filtering: For 1 ≤ k ≤ 10, the second to sixth columns report the average over B of the EDs between C k and respectively C b k i * , C b k hp with free parameter 10 5 , C b k hp with free parameter 5×10 5 , C b k ssa with window length 100, and finally C b k ssa with window length 200. The last two columns are the average over B of the ENs of Y b k i * and Y b k respectively.

	822 0.840 0.697 0.818 0.627 5.475 5.493
	2 0.887 0.930 0.785 0.920 0.729 2.854 2.922
	3 0.642 0.655 0.581 0.652 0.542 2.334 2.398
	4 4.808 6.014 4.898 3.924 3.070 49.73 49.63
	5 4.803 6.396 5.212 4.144 2.863 49.75 49.65
	6 8.513 7.315 6.123 6.265 13.75 49.45 49.61
	7 0.752 0.898 0.733 0.871 0.624 7.392 7.400
	8 0.631 0.393 0.286 1.651 0.211 17.32 17.31
	9 4.594 4.314 3.850 4.273 3.715 12.47 13.18
	10 6.369 7.028 5.732 3.924 3.070 49.84 49.66

  with those reported in Table4. For each k, we select the smallest p in Table3such that E

		p	E 1,p r	E 2,p r	E	3,p r	E	4,p r	E 5,p r	E 6,p r	E	7,p r	E	8,p r	E 9,p r	10,p E r
		1	5.346 5.819 3.441 5.421 6.350 21.62 1.524 0.655 21.43 9.033
		3	2.987 3.400 1.807 5.457 5.877 17.44 1.092 0.663 10.86 8.262
		5	1.709 1.618 0.917 6.037 6.115 13.20 1.020 0.863 6.927 7.972
		7	1.694 1.601 0.916 7.206 7.185 11.72 1.083 1.237 6.847 8.607
		8	1.474 1.341 0.838 8.150 8.031 12.31 1.109 1.514 6.558 9.297
		9	1.379 1.137 0.770 9.142 8.805 12.29 1.137 1.792 5.977 9.913
		10 1.366 1.101 0.759 10.18 9.822 12.76 1.191 2.101 5.752 10.73
		11 1.336 1.065 0.750 10.94 10.63 12.99 1.253 2.419 5.550 11.50
		12 1.371 1.068 0.760 11.88 11.54 13.73 1.335 2.766 5.506 12.37
		13 1.412 1.067 0.770 12.88 12.59 14.45 1.404 3.122 5.456 13.28
		14 1.453 1.078 0.782 13.84 13.62 15.23 1.491 3.482 5.463 14.16
		15 1.513 1.091 0.797 14.72 14.58 15.94 1.577 3.815 5.465 14.95
		16 1.574 1.112 0.814 15.44 15.35 16.54 1.671 9.125 5.508 15.61
		17 1.639 1.156 0.837 16.37 16.32 17.35 1.772 4.438 5.559 16.63
		18 1.706 1.211 0.858 17.24 17.14 18.08 1.884 4.739 5.629 17.43
		19 1.772 1.277 0.880 18.10 18.02 18.94 1.986 5.019 5.720 18.26
		20 1.856 1.379 0.912 18.74 18.75 19.65 2.097 5.314 5.832 18.98
		22 2.011 1.580 0.983 20.31 20.43 21.25 2.319 5.861 6.135 20.72
		24 2.174 1.763 1.059 21.91 22.08 22.97 2.503 6.264 6.535 22.49
		26 2.375 1.943 1.156 23.48 23.67 31.54 2.695 6.599 7.083 24.26
		28 2.597 2.076 1.258 24.98 25.33 26.15 2.907 6.931 7.639 25.85
		30 2.812 2.178 1.
									k,p gr < E	k g and denote it by p k 1 . We display
	E	k,p k 1 gr	in bold in Table 3 and we have p k 1 ∈ {11, 5, 10, 1, 3, 11, 16, 1, 13, 9}.
	For each k, we observe that for p > p k 1 , E	k,p gr varies (mostly decreases) until

Table 2 :

 2 Averaged EDs for ratio approach: For each k and for 27 selected fixed 1 ≤ p ≤ 45, this table reports the average over B of the EDs between C k and C

	b k ,p r	.

Table 3 :

 3 Averaged For each k, the bold averaged EDs associate with the smallest p where E

	k,p gr < E k g and the

EDs for energy-ratio approach: For each k and for 27 selected fixed 1 ≤ p ≤ 45, this table reports the average over B of the EDs between C k and C b k ,p gr .

Table 4 :

 4 Averaged EDs for energy approach: For each k, this table reports the average over B of the EDs between C k and C

	b k

g .

Table 5 :

 5 Search for i * in examples 1

	i †	1	2	3	4	5	6	7	8	9	10	11
	E 1 i † E Y 1 i † E 5 i † E Y 5 i †	5.49 2.1e-15 3.32 4.64 5.17 5.32 5.38 5.48 5.50 6.61 11.2 4.38 3.21 2.19 1.57 1.2 0.69 0.63 3.80 9.75 47.4 36.6 27.9 19.7 14.2 10.4 8.25 5.05 3.68 2.49 10.4 --7.5e-14 34.8 40.3 43.9 45.5 46.4 46.9 47.2 47.2 47.2 48.7

This document is a collaborative effort. * Corresponding author Email addresses: azadeh@mast.queensu.ca (Azadeh Moghtaderi), patrick.flandrin@ens-lyon.fr (Patrick Flandrin), pierre.borgnat@ens-lyon.fr (Pierre Borgnat)

B B b=1 G i,b .It is shown in[START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise, and Hurst exponent estimation[END_REF] that if the time series under study are realizations of a generic broadband process, then G i is a decreasing sequence in i. This results were obtained by studying fGn processes. This observation is also supported by the following data. Recall the broadband processes and their realizations from Section 4.1. We computed the IMFs of each realization along with G i,b and G i . Fig.3displays log

G i for 20 broadband processes in the collection. The result of this simulation supports the idea that G i is a decreasing sequence in i when computed for broadband data.Our key observation is that, generically, G i increases for i near the best index i * . This observation is supported by the following data. Recall the additive mixes obtained in Section 4.1. We computed the IMFs of each mix along with G i,b and for each broadband process in the collection, we computed G i . The left-hand plot in Fig.4displays log 2 G i computed for additive mixes. For each broadband process in the collection, we observe that G i increases at some i but we cannot yet determine whether or not it has occurred around i * . To clarify this, we detrended each mix as described in Section 4.1 and then recomputed G i,b using the remaining IMFs. For each

of the IMFs in the additive case are still valid for the multiplicative ones. To validate such properties, we proceeded by using the same simulations which were proposed in Sections 4.1 and 4.2 for broadband data and additive mixes.

We first recall 20 broadband processes and their realizations from Section 4.1. We computed the IMFs of the log-transform of the absolute value of each realization along with their zero crossing numbers and energy. We set R b and G i,b and then R and G i as described in Sections 4.1 and 4.2. The left-hand plot in Fig. 12 displays the empirical distribution of the elements of R and it supports the contention that R i ≈ 2. In fact, this distribution is approximately Gaussian with mean 2. The right-hand plot of Fig. 12 on the other hand displays log 2 G i and supports the contention that energy is a decreasing sequence in i for log-transformed broadband data.

Similar to the additive mixes, the approximation R i ≈ 2 expects to fail and also G i expects to increase for i near the best index i * . These observations are supported by the following data. For each broadband process in the collection and using its realizations, we constructed 10000 multiplicative mixes, using 1 + C 3 (displayed in Fig. 5) as a trend. We then computed the IMFs of the log-transformed absolute value of each mix and set R b , G i,b , R, and G i as described earlier. The top left-hand plot in Fig. 13 displays the empirical distribution of the elements of R and the top right-hand plot displays log 2 G i both for log-transformed multiplicative mixes. Similarly to Sections 4.1 and 4.2, for each broadband process in the collection, we used the IMFs obtained for each log-transformed mix and used the knowledge of


Hoel Le Capitaine
email: hoel.lecapitaine@univ-lr.fr

Carl Frélicot
email: carl.frelicot@univ-lr.fr

An optimum class-rejective decision rule and its evaluation

Keywords: bayesian classification, decision rules, loss structure, reject option, risk minimization

Decision-making systems intend to copy human reasoning which often consists in eliminating highly non probable situations (e.g. diseases, suspects) rather than selecting the most reliable ones. In this paper, we present the concept of class-rejective rules for pattern recognition. Contrary to usual reject option schemes where classes are selected when they may correspond to the true class of the input pattern, it allows to discard classes that can not be the true one. Optimality of the rule is proven and an upper-bound for the error probability is given. We also propose a criterion to evaluate such class-rejective rules. Classification results on artificial and real datasets are provided.

I. INTRODUCTION AND MOTIVATION

Let Ω = {ω 1 , ..., ω c } be a set of classes with prior probabilities p(ω i ) and assume that, given a pattern x in a feature space X, the conditional densities f (x|ω i ) and the mixture density f (x) = c i=1 p(ω i ) f (x|ω i ) can be calculated (or estimated). The probability that x belongs to class ω i is given by

p(ω i |x) = p(ω i ) f (x|ω i ) f (x) . (1)
In the sequel, we assume that we use a Bayesian classifier, and focus our study on posterior probabilities. If risk(x) denotes the risk of taking a wrong decision for x, the so-called Bayes decision rule, which assigns x to the class of maximum posterior probability [START_REF] Chow | An optimum character recognition system using decision functions[END_REF], minimizes the error probability defined by e = X risk(x) f (x) dx.

(2)

However, when some classes are not known or when classes overlap in X, this rule may not be always efficient. Chow modified it so that one can reject a pattern if its maximum posterior probability is less than a user-defined threshold t [START_REF] Chow | An optimum character recognition system using decision functions[END_REF]. Chow's rule minimizes (2) as a function of t which specifies the reject probability and allows to assign x to n ∈ {1, c} classes. Ha proposed a class-selective rule which minimizes [START_REF] Ha | The optimum class-selective rejection rule[END_REF] for an average number of selected classes n [START_REF] Ha | The optimum class-selective rejection rule[END_REF]. It allows to assign x to n ∈ [1, c] classes. In recent years, research on reject option has focused on binary classification, see e.g. [START_REF] Yuan | Classification methods with reject option based on convex risk minimization[END_REF], [START_REF] Grandvalet | Support vector machines with a reject option[END_REF], neglecting the multi-class case, with the notable exception of [START_REF] Tax | Growing a multi-class classifier with a reject option[END_REF], where ad-hoc transformations are operated on both distance and density base classifiers.

In this paper, the multi-class case is considered without requiring a specific transformation on output models. We are interested in applications, e.g. diagnosis of diseases or police investigation, where it can be more efficient to first discard less reliable cases, e.g. diseases or suspects, instead of selecting most reliable ones. To this aim, we introduce the concept of class-rejective rules: classes are eliminated and x is assigned to the remaining ones. In order to obtain good performances with respect to the costs associated with specific actions, one may reject as much classes as possible, while preserving a low error rate. Moreover, such a rule is able to reject all the c classes if none of them corresponds to x, resulting to an usual distance reject scheme [START_REF] Dubuisson | A statistical decision rule with incomplete knowledge about classes[END_REF].

II. DECISION RULE AND OPTIMALITY

Let α j = α j (x) be the action of selecting, for x, a subset S j of classes among the 2 c possible ones, and L(α j |ω i ) the loss incurred by action α j when the true class of x is ω i . The overall risk R to be minimized with respect to actions α j is the average expected loss

R = X R(α j |x) f (x) dx (3)
where R(α j |x) is the conditional risk defined by

R(α j |x) = c i=1 L α j (x)|ω i p(ω i |x). (4)
We define a loss structure modeling the specifical needs of the class-rejection problem. Let R j be the subset of r j rejected classes and set L(α

j |ω i ) = L e (α j |ω i ) + L r (α j )
where L e (α j |ω i ) and L r (α j ) denote respectively the loss of eliminating the true class ω i and the loss of having eliminated a few number of classes, when taking action α j . We propose to define both losses as:

L e (α j |ω i ) = C e if ω i ∈ R j 0 otherwise (5) L r (α j ) = C r (c -r j ) (6)
where C e and C r are user-defined positive costs depending on the application but not on the decision rule. Thus, (4) can be rewritten as

R(α j |x) = C e ωi∈Rj p(ω i |x) + C r (c -r j ). (7)
Since p(x) > 0 for all x, minimizing (3) is equivalent to minimizing (4) or ( 7) over all possible numbers of rejected classes r ∈ [0, c] and we have to solve

min r∈[0,c] min rj =r R(α j |x) . (8)
For a fixed r, it is obvious that the solution of ( 8) is such that R j contains the r classes of lowest posterior probabilities.

Let us use the decreasing ordered sequence of posterior probabilities such as p(ω

(i) |x) ≥ p(ω (i+1) |x), ∀i = 1, c.
Then, the conditional risk ( 7) can be defined as a function of the r worst classes by

R(r|x) = C e c i=c-r+1 p(ω (i) |x) + C r c -r . (9)
Setting t = C r /C e and simplifying by Ce, (8) becomes

min r∈[0,c] R(r|x) = min r∈[0,c] c i=c-r+1 p(ω (i) |x) + t (c -r) . (10)
The first term in the right-hand side of ( 10) is (not strictly) convex because it sums increasing values and the second term t(cr) is a (not strictly) convex function, so R(r|x) is convex. Therefore, if there is a unique solution r , it satisfies:

R(r -1|x) > R(r |x) (11) R(r + 1|x) > R(r |x). (12)
Using [START_REF] Asuncion | UCI machine learning repository[END_REF], it results in p(ω (c-r +1) |x) < t and p(ω (c-r ) |x) > t. Since p(ω (i) |x)'s are ordered, it ensues the optimal number of rejected classes:

r (x, t) = max r∈[0,c] r : p(ω (c-r+1) |x) < t (13)
or, equivalently,

r (x, t) = min r∈[0,c] r : p(ω (c-r) |x) > t (14)
with p(ω (c+1) |x) = 0 and p(ω (0) |x) = 1 by convention.

What about several solutions for (10)? Since R(r|x) is convex, multiple solutions are necessary consecutive and taking the highest (respectively the lowest) only change '>' to '≥' in [START_REF] Highleyman | Linear decision functions, with application to pattern recognition[END_REF] (respectively ( 12)) and '<' to '≤' in [START_REF] Vanderlooy | A critical analysis of variants of the auc[END_REF] (respectively '>' to '≥' in ( 14)).

With this formulation, the risk associated to the optimal number of rejected classes r can be written as

risk(x, r ) = c i=c-r +1 p(ω (i) |x) (15)
Special values of the solution can be emphasized:

• r = 0: x is associated to all classes as it would be using Chow's rule [START_REF] Chow | An optimum character recognition system using decision functions[END_REF],

• r = c -1: x is rejected from all classes except one, so the rule reduces to the Bayes decision rule and the risk ( 15) is the Bayes one risk(x) = 1 -p(ω (1) |x), • r = c: x is rejected from all (known) classes, i.e. it is distance rejected in the sense of [START_REF] Dubuisson | A statistical decision rule with incomplete knowledge about classes[END_REF]. As a final proposition, let us give an upper-bound of the induced probability error e(t). From ( 13) and the ordering of p(ω (i) |x)'s, we have

t ≥ p(ω (c-r +1) |x) ≥ ... ≥ p(ω (c) |x). (16)
Therefore,

c i=c-r +1 t ≥ p(ω (c-r +1) |x) + ... + p(ω (c) |x) ⇔ r × t ≥ risk(x, r ). (17)
Integrating both sides, we finally get

t X r (x, t) f (x) dx ≥ X risk(x, r ) f (x) dx ⇔ t × r (t) ≥ e(t) (18)
where r is the average number of rejected classes.

III. EVALUATION OF CLASS-REJECTING RULES

In order to evaluate a decision rule whose parameters are Θ, one generally use: the correct classification C(Θ), the misclassification (or error) E(Θ) and the reject R(Θ) probabilities or rates. Chow introduced the error-reject (ER) trade-off (both should be minimized, see also [START_REF] Hansen | The error-reject tradeoff[END_REF]) and proposed to analyze the so-called (E, R)-curve (E(Θ) vs. R(Θ)) for all possible values of Θ = t to find the optimal or an operational value for t [START_REF] Chow | An optimum character recognition system using decision functions[END_REF]. When comparing rules, the least the AU C (Area Under Curve) is, the better the rule, see [START_REF] Fawcett | An introduction to roc analysis[END_REF] for its description in ROC analysis. For classselective rules, e.g. Ha's rule [START_REF] Ha | The optimum class-selective rejection rule[END_REF], an object is considered as misclassified if the class it is issued from does not belong to the subset of selected ones. Then, the (E, R) trade-off is replaced by the (E, n) trade-off, where n is the average number of selected classes, and the AU C(E, n) can be used to compare rules [START_REF] Capitaine | Aggregation operators for similarity measures. application to ambiguity in pattern recognition[END_REF].

For the new kind of class-rejective rules we introduce, we set that an object is misclassified if the class it is issued from belongs to the subset of rejected classes, and we introduce the (E, r) trade-off which consists in jointly minimizing E(Θ) and maximizing the average number of rejected classes r(Θ). Given such rules, the best one is the rule which has the least AU C(E, r):

min available rules AU C(E, r) = min available rules D E r (t) dt (19)
where D is the definition domain for t. For a single Bayesian classifier, we compare the (E, r) curves, which are the counterparts of Chow's error-reject trade off curves, obtained using the available rules.

IV. EXPERIMENTAL RESULTS

With respect to AU C(E, r), we compare the new Optimal Class-Rejective (OCR) rule to what we call the Bottom-r Ranking (BrR) rule. By analogy to the top-n ranking rule used by Ha [START_REF] Ha | The optimum class-selective rejection rule[END_REF], it simply consists in setting a constant number of rejected classes r for the whole data set, i.e. posterior probabilities are not considered. The error rate is then computed by considering that the r lowest probabilities correspond to the r rejected classes.

Both rules are tested on fourteen datasets whose characteristics (number n of patterns, number p of features, number c of classes, degree of overlap) are reported in Table I. Eleven datasets are from the UCI ML-Repository [START_REF] Asuncion | UCI machine learning repository[END_REF], three are the following synthetic ones:

• DH consisting of two overlapping Gaussian classes with different covariance matrices according to the Highleyman distribution, each composed of 800 observations in R 2 , see [START_REF] Highleyman | Linear decision functions, with application to pattern recognition[END_REF]. • D1 containing 2000 points drawn from two Gaussian seven-dimensional distributions of 1000 points each with means v 1 = t (1, 0, ..., 0) and v 2 = t (-1, 0, ..., 0), and equal covariance matrices Σ 1 = Σ 2 = I.

• D2 which contains 4000 points drawn from four Gaussian bi-dimensional distributions of 1000 points each. Their respective means are v 1 = t (1, 1), v 2 = t (1, -1), v 3 = t (-1, 1) and v 4 = t (-1, -1), and their covariance matrices are all equal

Σ 1 = Σ 2 = Σ 3 = Σ 4 = I.
Classes composing each real dataset are also assumed to be Gaussian. Note that this does not limitate the OCR rule which only need posterior probabilities to be computed using (1) whatever the class-conditional density model. The threshold t is varied within the range [0, 1] for the OCR rule and the number of rejected classes is discretely increased from 0 to c -1 (to skip the distance reject case) for the BrR rule. For both rules, error rates are computed through a resubstitution procedure, and the performance is given in Table II where best scores are indicated in bold.

These results show that the proposed OCR rule always outperforms the BrR rule.

Looking at the amount of overlap between classes, one can go deeper into the analysis of the relative scores considering their ratio AU C(BrR)/AU C(OCR) and their difference The reason is that for better separated classes, there are more patterns for which the maximum posterior probability is much higher than the others, so the associated risk (and implicitly the error) is lower for the OCR rule whereas this is not taken into account for the BrR rule. While the ratios are lower for datasets presenting a strong overlap, the differences are much more important (except for Yeast) than for datasets presenting better separated classes. This means that the less the classes overlap, the less the benefit is, as it is for any decision rule including a reject option, e.g. [START_REF] Chow | An optimum character recognition system using decision functions[END_REF], [START_REF] Dubuisson | A statistical decision rule with incomplete knowledge about classes[END_REF], [START_REF] Ha | The optimum class-selective rejection rule[END_REF].

Figure 1 shows the (E, r) curves obtained on the three artificial datasets. One can see that for every average number of rejected classes, the error rate obtained with the proposed OCR rule is lower than the one obtained with the BrR rule. Obviously, when the number of rejected classes is equal to c -1 (respectively 1,1 and 3 for D1, D2 and DH), the error rate is the same for both rules as one could expect. As pointed out in the previous section, this case corresponds to the Bayes decision rule, so that both rules are the same one.

V. CONCLUSION

In this paper, we propose the new concept of classrejective decision rule for pattern recognition. Instead of selecting the most probable classes, it allows to discard (reject) the less probable ones. Its optimality with respect to the error probability for an average number of rejected classes is proven as well as an upper-bound to the error probability. A measure to evaluate such rules is introduced, similarly to the one used to evaluate class-selective decision rules. Classification performance obtained on artificial and real datasets show that it outperforms the rule which consists in rejecting a constant number of classes. The proposed rule can be used with any classifier, provided posterior or equivalent values (e.g. membership degrees from a fuzzy classifier) are available. Future research will concern the study of a new mixed class-selective-rejective decision rule which should jointly optimize the number of selected and rejected classes over the user-defined threshold definition domain. We also plan to use this optimum decision rule for outliers detection, i.e. patterns that do not match any of the known classes so that they must be distance rejected. This problem, as well as its variant for support vector machines, will be addressed using hinge loss minimization [START_REF] Bartlett | Classification with a reject option using a hinge loss[END_REF]. It will be studied by taking into account new results on AUC variants [START_REF] Vanderlooy | A critical analysis of variants of the auc[END_REF].

Figure 1 .

 1 Figure 1. (E, r) curves of both OCR and BrR rules curve obtained on the three artificial datasets (from left to right: DH, D1 and D2).

Table I DESCRIPTION

 I OF THE TESTED DATASETS.

	Dataset	n	p	c	Overlap
	DH	1600	2	2	slight
	D1	2000	7	2	moderate
	D2	4000	2	4	strong
	Ionosphere	351	34	2	strong
	Forest	495411 10	2	strong
	Vowel	528	10 11 moderate
	Digits	10992	16 10	slight
	Thyroid	215	5	3	slight
	Iris	150	4	3	slight
	Pima	768	9	2	strong
	Statlog	6435	36	6	moderate
	Glass	214	9	6	strong
	Yeast	1484	8	10	strong
	Page-blocks	5473	10	5	moderate
		Table II
	AU C(E, r) FOR THE OPTIMAL CLASS-REJECTIVE AND THE BOTTOM-r
		RANKING RULES.
	Dataset	OCR Rule B-rR Rule
	DH	0.0053		0.0262
	D1	0.0406		0.0727
	D2	0.1794		0.2373
	Ionosphere	0.1644		0.3205
	Forest	0.1309		0.1577
	Vowel	0.0374		0.0473
	Digits	0.0098		0.0319
	Thyroid	0.0064		0.0302
	Iris	0.0092		0.0200
	Pima	0.1146		0.1628
	Statlog	0.0959		0.1490
	Glass	0.1133		0.1472
	Yeast	0.7572		0.7716
	Page-blocks	0.0906		0.0987
	AU C(BrR) -AU C(OCR). The ratios are much more
	important for datasets presenting a slight overlap of classes
	{DH,Digits,Thyroid,Iris} than for datasets presenting a
	strong overlap {D2,Ionosphere,Forest,Pima,Glass,Yeast}.

ACKNOWLEDGMENTS

Support for part of this work was provided by the French ministry of research through a PhD grant.


Hoel Le Capitaine
email: hoel.lecapitaine@univ-lr.fr

Carl Frélicot
email: carl.frelicot@univ-lr.fr

On selecting an optimal number of clusters for color image segmentation

Keywords: cluster validity, clustering methods, color image segmentation, overlap and separation measures

This paper addresses the problem of region-based color image segmentation using a fuzzy clustering algorithm, e.g. a spatial version of fuzzy c-means, in order to partition the image into clusters corresponding to homogeneous regions. We propose to determine the optimal number of clusters, and so the number of regions, by using a new cluster validity index computed on fuzzy partitions. Experimental results and comparison with other existing methods show the validity and the efficiency of the proposed method.

I. INTRODUCTION

Region-based image segmentation consists in partitioning an image into non-intersecting regions such that pixels within a region are homogeneous and pixels from adjacent regions are not. Fuzzy clustering algorithms have become popular, e.g. fuzzy c-means (FCM) [START_REF] Bezdek | Pattern Recognition with fuzzy objective function algorithm[END_REF]. Like any unsupervised algorithm, it requires the number of clusters to be set by the user. This makes FCM unsuited for segmenting image databases unless this number is automatically determined. Cluster validity indexes (CVIs) come to this end and once an optimal number of clusters for a dataset is obtained, the corresponding fuzzy partition is selected as the optimal one [START_REF] Bezdek | Some new indexes of cluster validity[END_REF], [START_REF] Kim | On cluster validity index for estimation of the optimal number of fuzzy clusters[END_REF], [START_REF] Pal | On cluster validity for the fuzzy cmeans model[END_REF]. Usual FCM and its spatial version [START_REF] Chuang | Fuzzy c-means clustering with spatial information for image segmentation[END_REF] we use are presented in section II. In section III, we recall basic definitions on fuzzy aggregation operators [START_REF] Klement | Logical, Algebraic, Analytic, and Probabilistic Aspects of Triangular Norms[END_REF] and present a CVI to be used in an automatic segmentation process. Experimental results are given in section IV, comparison to other unsupervised segmentation algorithms is provided. We finally conclude and draw some perspectives in section V.

II. FUZZY CLUSTERING AND IMAGE SEGMENTATION

Clustering is an instance of unsupervised classification which aims at finding a structure of groups in a set of n p-dimensional points X = {x 1 , ..., x n }. In this framework, the label vectors u k = u(x k ) do not exist and clustering algorithms can be used to obtain them from X. For instance, the fuzzy c-means (FCM) algorithm [START_REF] Bezdek | Pattern Recognition with fuzzy objective function algorithm[END_REF] partitions X into c > 1 clusters by minimizing the following objective function

J m (U, V ) = n k=1 c i=1 u m ik ||x k -v i || 2 (1)
where u ik is the membership degree of x k to the i th cluster represented by its centroid

v i ∈ R p . Centroids are gathered into a (c×p) matrix V = [v 1 , ..., v c ]. Degrees u ik are subject to c i=1 u ik = 1 for all x k in X and to 0 < n k=1 u ik < n (∀i = 1, • • • , c).
In addition, they are elements of the fuzzy c-partition matrix U (c × n). The so-called fuzzifier m > 1 is a weighting exponent which makes the resulting partition more or less fuzzy [START_REF] Pal | On cluster validity for the fuzzy cmeans model[END_REF]. The higher m is, the softer the clusters' boundaries are. Minimization of J m (U, V ) is obtained by iteratively and alternatively updating (U, V ). Application to image segmentation consists in taking X as the entire set of pixels x k in a color space, say R 3 , and obtain r regions given by (U, V ), where r ≥ c since several non adjacent regions in the image can have the same centroid in the color space. Unfortunately, FCM fails when noisy pixels are present so a filtering process is required. An alternative is to add spatial constraints in order to smooth heterogeneous pixel values effect within a region, e.g. the simple sFCM algorithm [START_REF] Chuang | Fuzzy c-means clustering with spatial information for image segmentation[END_REF] where the degrees of neighboring pixels are used to modify U during the updating step as follows :

u ik = u p ik h q ik c j=1 u p jk h q jk ( 2
)
where h is a spatial function and (p, q) are parameters which control the relative importance of u jk and h jk so that the smoothing effect due to the neighborhood increases with q.

The general form of h is given by h ik = j∈N (x k ) u ij , where N (x k ) is the set of neighboring pixels, e.g. a square window centered on x k . Note that sFCM reduces to FCM when (p, q) = (1, 0).

III. THE CLUSTER VALIDITY INDEX

For the CV application, we are interested in aggregation functions that map a collection u in [0, 1] to a value in [0, 1], formally a vector u = t (u 1 , • • • , u c ) → A(u). Among the frequently used aggregation operators, one finds the class triangular norms (t-norms) and conorms (t-conorms). They have been introduced to characterize the general multivalued logic AND and OR operations and are widely used in fuzzy logic and fuzzy set theory to implement conjunctive and disjunctive operators respectively. A t-norm is a commutative, associative and monotone function having for neutral element 1. Alternatively, its dual operator, the t-conorm, is a commutative, associative and monotone function ⊥ having for neutral element 0. Examples of norms couples are given in Table I for two operands, see [START_REF] Klement | Logical, Algebraic, Analytic, and Probabilistic Aspects of Triangular Norms[END_REF] for a large survey, including parameterized families. The dual couple is generally associated with a strong fuzzy negation defined as N (v) = 1 -v and mentioned as the triple ( , ⊥, N ).

Table I USUAL TRIANGULAR NORM COUPLES Standard u 1 S u 2 = min(u 1 , u 2 ) u 1 ⊥ S u 2 = max(u 1 , u 2 ) Algebraic u 1 A u 2 = u 1 u 2 u 1 ⊥ A u 2 = u 1 + u 2 -u 1 u 2 Łukasiewicz u 1 L u 2 = max(u 1 + u 2 -1, 0) u 1 ⊥ L u 2 = min(u 1 + u 2 , 1)
Given a fuzzy clustering (U, V ) of X, each u ik defines the similarity of the object x k to the prototype v i , or the degree to which x k satisfies the i th group description. Using the membership vector

u k = t (u 1k , • • • , u ck ) (k th column of U )
, the standard t-conorm (max operator) is commonly used to select the most appropriate group the object should be assigned to. However, such an exclusive partitioning is not efficient because lower values interact with the greatest value, meaning that x k satisfies more than one group description. Therefore, it is adequate to aggregate the values in u k in order to assess to which extend x k belongs or not to several clusters and how many. We propose to use the fuzzy OR operator (fOR-l for short) defined in [START_REF] Mascarilla | A k-order fuzzy or operator for pattern classification with k-order ambiguity rejection[END_REF] for supervised classification with reject options purpose. This operator evaluates degrees of similarity at a given order (l) by combination of t-norms. Let P be the power set of C = {1, 2, ..., c} and P l = {A ∈ P : |A| = l} where |A| denotes the cardinality of subset A, then the fOR-l associates

to u k a single value l ⊥(u k ) ∈ [0, 1] defined by: l ⊥(u k ) = l ⊥ i=1,••• ,c u ik = A∈P l-1 ⊥ j∈C\A u jk (3)
It must be viewed as some kind of generalization of the "l th highest" value, l ∈ C. Using standard t-norms, l ⊥(u k ) is exactly the l th highest element of u k (see [START_REF] Mascarilla | A k-order fuzzy or operator for pattern classification with k-order ambiguity rejection[END_REF]).

Validating the provided clustering (U, V ) of X consists in assessing whether the resulting partition reflects the data structure or not. Due to the unsupervision, the user does not have any prior knowledge on c whereas it is required by clustering algorithms such as FCM. Most of works on CV focus on the number of clusters problem. A reliable CVI for FCM must consider both compactness and separation of the resulting fuzzy c-partition U . If only the former is considered, the best partition consists of as many clusters as points. On the other hand, if only the latter is considered, the best partition consists of a single cluster. So it is well established that a CVI has to combine both measures to be efficient. Separation measures generally take distances between centroids into account but it is seldom sufficient to interpret the geometrical structure of the data, and therefore the separation between clusters, see [START_REF] Kim | On cluster validity index for estimation of the optimal number of fuzzy clusters[END_REF] for examples on the limitations of these approaches. We propose to use, for each point x k , two measures that overcome these drawbacks: a fuzzy overlap measure which evaluates the degree of overlap of a specified number (l) of fuzzy clusters and a fuzzy separation measure which indicates the degree of overlap of the most probable cluster, i.e. the one corresponding to the highest membership degree, with respect to the c -1 other ones. A low value of this latter measure will denote a large separation of the most probable cluster from the others. Despites the ability to deal with overlapping clusters is now considered to be a major quality of a CVI [START_REF] Bouguessa | An objective approach to cluster validation[END_REF], most of works are based on an intuitive representation of overlap. An overlap measure between l fuzzy clusters for each point x k of X described by its membership vector u k can be obtained by [START_REF] Kim | On cluster validity index for estimation of the optimal number of fuzzy clusters[END_REF]. By successively computing l ⊥(u k ) for different values of l, we get a combination of l-order overlap degrees for x k . In order to determine the overall degree of overlap for a given point, it is adequate to determine which orders induce high overlap. The most satisfied order(s) can be quantified by the fuzzy disjunction of the l-order overlap measures (l = 2, c), so we define the overall overlap measure of a point x k by:

O ⊥ (u k , c) = 1 ⊥ l=2,c   l ⊥ i=1,c u ik   (4)
In [START_REF] Bezdek | Some new indexes of cluster validity[END_REF], Bezdek and Pal show that inter cluster separation plays a more important role in cluster validity than diameters. We propose to introduce such a measure by quantifying the fuzzy separation of each point x k with 1 ⊥(u k ), which is the overlap measure within one fuzzy cluster, i.e. its separation from the other fuzzy clusters, since u k sums up to one. This aggregation corresponds to the fuzzy disjunction of membership degrees for a given point x k , which selects the most probable cluster. We define the fuzzy separation of x k with respect to the c clusters as:

S ⊥ (u k , c) = 1 ⊥      1 ⊥ i=1,c u ik , • • • , 1 ⊥ i=1,c u ik c-1 times      (5)
Since a small value of the overlapping degree O ⊥ (u k , c) and a large value of the separation one S ⊥ (u k , c) indicates that the highest membership degree of x k corresponds to a well separated and not overlapping cluster, we finally define the l-order Overlap and Separation Index (OSI) taking values in [0, 1] as the average value of the ratios of both measures:

OSI ⊥ (c) = 1 n n k=1 O ⊥ (u k , c) S ⊥ (u k , c) (6)
Running sFCM for c in a range [c min , c max ] and selecting the value c best which minimizes ( 6) allows to segment an image in r best corresponding homogeneous regions.

IV. COLOR IMAGE SEGMENTATION RESULTS

Various studies have shown that many color spaces proposed for computer graphic applications are not well adapted to image processing. As pointed out in [START_REF] Angulo | Color segmentation by ordered mergings[END_REF], a convenient representation should yield distances and provides independence between chromatic and achromatic components. For this reason and comparison purpose, we use the CIElab color space, see [START_REF] Mignotte | Segmentation by fusion of histogram-based kmeans clusters in different color spaces[END_REF] for details and comparisons. We test the proposed segmentation approach on the Berkeley segmentation database [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] which consists in 300 color images. For each image, a set of benchmark results provided by human observers is given. Performance is evaluated by computing the consistency of the segmentations against the human ones using some classical criteria: -the Probabilistic Rand Index (PRI) which counts the number of pixels whose labels are the same for both segmentations, -the Variation of Information (VoI) which averages the conditional entropy of one segmentation given the other one, -the Global Consistency Error (GCE) which evaluates to which extent one segmentation can be viewed as a refinement of another. Many unsupervised image segmentation algorithms have been proposed so far. Since comparison cannot be done with all of them, we choose to compare the proposed method to three couples of methods: 1) Two algorithms that are widely used in the segmentation literature, which complexity is comparable to the CV approach: the Mean-shift [START_REF] Comaniciu | Mean-shift: A robust approach toward feature space analysis[END_REF] and N-cuts [START_REF] Shi | Normalized cuts and image segmentation[END_REF]. Meanshift is a robust effective mode-seeking algorithm based on nonparametric kernel density estimation. The optimal number of modes, depending on some resolution parameters, is automatically performed so segmentation consists in associating the pixels to the closest local mode. N-cuts (Normalized-cuts) is a graph-theoretic method taking each pixel as a node and connecting each pair of pixels by an edge so that the weight on the edge reflects the likelihood that the two pixels belong to the same region. Segmentation is obtained by cutting the edges less than an optimal value, solving an eigenvalue problem. For both algorithms, internal parameters are set to their default values.

2) Two more sophisticated and somehow new methods: the Compression-based Texture Merging [START_REF] Ma | Segmentation of multivariate mixed data via lossy coding and compression[END_REF] and Fusion of Clustering Results [START_REF] Mignotte | Segmentation by fusion of histogram-based kmeans clusters in different color spaces[END_REF]. CTM is based on local features assumed to be drawn from a mixture of (possibly degenerate) Gaussians. Component parameters are estimated using the lossy minimum description length and the segmentation is the one that minimizes the number of bits needed to code the segmented data, subject to a given distortion. FCR fuses segmentation maps obtained through the k-means clustering of histograms computed on various color spaces.

3) sFCM associated with two popular CVis: the Kwon (K) and Xie-Beni (XB) indexes, see [START_REF] Kim | On cluster validity index for estimation of the optimal number of fuzzy clusters[END_REF] for definitions.

For the CV approach, sFCM is used with a square window 7 × 7 for the neighborhood and (p, q) = (0, 2) so that the smoothing effect is large enough to remove small pixel variations and the computational time is reasonable.

Results reported in Table II show that the proposed method (sFCM + OSI ⊥ S ) obtains the overall best mean rank, the best score for VoI, and it outperforms the other CV approaches (group 3) whatever the performance measure. It also gives better results than the non clustering methods (group 1) for two measures. Compared to the more highly sophisticated methods (group 2), it gives quite satisfactory in-between scores. Note that for sFCM associated with the presented CVI, only the results obtained using the standard norms which aggregate cautiously are given. Using algebraic and Łukasiewicz norms sometimes produce over-and undersegmented results respectively, because they involve input values interactions and are Archimedean ( (u 1 , u 1 ) < u 1 ) whereas standard ones are not.

V. CONCLUSION

A non-parametric segmentation method based on a cluster validity index, provided fuzzy partitions of pixels in a given color space is proposed. It is associated to a spatial version of the fuzzy c-means algorithms in order to take into account neighboring pixels. The index allows to select a number of clusters that are as separated and non overlapping as possible and to use it to segment an image in homogeneous regions. Results show the efficiency of the proposed approach, compared to other unsupervised segmentation algorithms. We intend to adapt it to the color space selection problem to compare it to other sophisticated methods. In a future work, we will intensively study how mathematical properties of triangular norms affect the selected partition, and therefore the segmentation.

Figure 1 .

 1 Figure 1. Images segmented by sFCM+OSI ⊥ S .