Dominique Pastor
email: dominique.pastor@enst-bretagne.fr

Asmaa Amehraye
email: asmaa.amehraye@enst-bretagne.fr

Algorithms and Applications for Estimating the Standard Deviation of AWGN when Observations are not Signal-Free

Keywords: Binary hypothesis testing, decision, estimation, likelihood theory, multivariate normal distribution, speech denoising

Consider observations where random signals are randomly present or absent in independent and additive white Gaussian noise (AWGN). By using a recently established limit theorem, we introduce a new estimator for the estimation of the noise standard deviation when the signals are less present than absent and have unknown probability distributions.

The bias, the consistency and the minimum attainable mean square estimation error of the estimator we propose are still unknown. However, the experimental results that are presented are very promising. First, when the Minimum-Probability-of-Error decision scheme for the non-coherent detection of modulated sinusoidal carriers in independent AWGN is tuned with the outcome of our estimator instead of the true value of the noise standard deviation, the Binary Error Rate tends to the optimal error probability when the number of observations is large enough. Second, given some speech signal corrupted by independent AWGN, our estimator can be used to estimate the noise standard deviation so as to adjust the standard Wiener filtering of the noisy speech. The objective performance measurements obtained by so proceeding are very close to those achieved when the Wiener filtering is tuned with the true value of the noise standard deviation.

I. INTRODUCTION

In many signal processing applications, observations can be modelled by d-dimensional random vectors that result from the random presence of signals in independent and additive white Gaussian noise (AWGN).

Often, very little is known about the signals or about most of their describing parameters. The probability distributions of the signals are sometimes partially or definitely unknown. This issue is typically met in Electronic (Warfare) Support Measure (ESM) systems faced with non-cooperant communications. This is also the case with sonar systems that receive signals resulting from noise generated by motors and hull vibrations transmitted through a fluctuating environment. Echoes received by a radar system are not always easy to model because they derive from a convolution between a known transmitted pulse and an unknown environment.

In addition to our lack of prior knowledge about the signals, the noise standard deviation may be unknown as well. In order to process the observations, it can then be necessary to estimate it.

For instance, Constant False Alarm Rate (CFAR) systems, standardly used in radar processing, perform an estimate of the noise standard deviation in order to detect radar targets with a false alarm rate close to some prespecified value. Basically, CFAR systems select observations considered as signal-free and called "reference cells". The estimation of the noise standard deviation is then achieved on the basis of these reference cells. The design of a CFAR system is no easy task (see [START_REF] Minkler | The Principles of Automatic Radar Detection in Clutter[END_REF], [START_REF] Levanon | Radar Principles[END_REF]).

A second example is the denoising of speech signals corrupted by independent AWGN. Indeed, most standard filtering techniques must be adjusted with an estimate of the noise standard deviation when the latter is unknown. A popular approach involves using a Voice Activity Detector (VAD). The VAD decides, with some error rate, that some periods of time are signal-free and the noise standard deviation is then estimated on the basis of these periods of time. Another approach is to achieve an estimate of the noise standard deviation by computing the smallest eigenvalues of the noisy speech autocorrelation matrix. Depending on the matrix size, the computation of the eigenvalues may prove unstable.

In contrast with such approaches mostly based on the physics of the problem, we hereafter propose an estimator relying on statistics only. This estimator is based on a recent theoretical result established in [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF]. A simplified version of this result will however be sufficient to derive our estimator.

II. PAPER ORGANIZATION

After introducing some notations and pieces of terminology, section III presents the theoretical result, namely proposition 3.1, from which the estimator proposed in this paper is derived. Proposition 3.1 and the more general result [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF]Theorem 1] concern any sequence of mutually independent d-dimensional real random observations where signals, with unknown probability distributions but with norms larger than or equal to some known value, are randomly present or absent with unknown probabilities of presence less than or equal to one half in independent AWGN.

Our estimator is then presented in section IV. It applies to sets of non signal-free observations where signals have unknown probability distributions and unknown probabilities of presence less than or equal to one half. No prior knowledge on a lower bound for the norms of the signals is any longer required.

The bias, the consistency and the mean square error of this estimator are still unknown. Nevertheless, several experimental results are presented in section V. These experimental results concern two-dimensional observations where the signals are uniformly distributed on the circle centred at the origin with known radius α. The two components of every signal can be regarded as the in-phase and quadrature components of a modulated sinusoidal carrier with amplitude α and phase uniformly distributed in [0, 2π].

An application to speech denoising is also presented. Given some speech signal corrupted by independent AWGN, section VI describes how the estimator introduced in this paper makes it possible to estimate the noise standard deviation and adjust the standard Wiener filtering of the noisy speech signal.

Section VII concludes this paper by suggesting several theoretical developments, extensions to speech processing as well as applications in radar processing and electronic warfare.

III. THEORETICAL BACKGROUND

In what follows, only one probability space (Ω, M, P) is considered and every random vector or variable is assumed to be defined for every ω ∈ Ω by setting this random vector or variable to 0 on any negligible subset where it could be undefined. As usual, we write (a-s) for "almost surely".

Let S stand for the set of all the sequences of ddimensional real random vectors. Given a positive real value σ 0 , an element X = (X k ) k∈N of S will be called a d-dimensional white Gaussian noise (WGN) with standard deviation σ 0 if the random vectors X k , k = 1, 2, . . ., are mutually independent and identically Gaussian distributed with null mean vector and covariance matrix σ 2 0 I d where I d stands for the d × d identity matrix. For the sake of brevity, we will henceforth write independent instead of mutually independent.

We define the minimum amplitude of an element S = (S k ) k∈N of S as the supremum a(S) of the set of those α ∈ [0, ∞] such that, for every natural number k, S k is larger than or equal to α (a-s):

a(S) = sup {α ∈ [0, ∞] : ∀k ∈ N, S k ≥ α (a-s)} . ( 1
)
If f is some map of S into R, we say that the limit of f (S) is ∈ R when a(S) tends to ∞ and write that lim a(S)→∞ f (S) = if, for any positive real value η, there exists some α 0 ∈ (0, ∞) such that, for every α ≥ α 0 and every S ∈ S with a(S) ≥ α, we have |f (S) -| ≤ η.

Given ν ∈ [0, ∞), let L ν (Ω, R d ) stand for the set of those d-dimensional real random vectors Y such that E[ Y ν ] < ∞.
In the sequel, we will encounter the set

∞ (N, L ν (Ω, R d )) of those elements S = (S k ) k∈N of S such that S k ∈ L ν (Ω, R d ) for every k ∈ N and sup k∈N E[ S k ν ] is finite.
In what follows, 0 F 1 is the generalized hypergeometric function (see [4, p. 275]). Given ρ ∈ [0, ∞), ξ(ρ) is the unique positive solution for x in the equation

0 F 1 (d/2; ρ 2 x 2 /4) = exp(ρ 2 /2). (2)
Given q ∈ [0, ∞), Υ q stands for the map defined for every x ∈ [0, ∞) by

Υ q (x) = x 0 t q+d-1 exp(-t 2 /2)dt. (3)
The following result is a particular case of a more general limit theorem established in [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF]. The statement considered in the present paper is sufficient with respect to our purpose. Given any random vector Y and any real number τ , I( Y ≤ τ ) stands for the indicator function of the event { Y ≤ τ }.

Proposition 3.1: Let Y = (Y k ) k∈N be some element of S such that, for every k ∈ N, Y k = ε k S k + X k where S = (S k ) k∈N , X = (X k ) k∈N and ε = (ε k ) k∈N are an element of S, some d-dimensional WGN with standard deviation σ 0 and a sequence of random variables valued in {0, 1} respectively.

Assume that (A1) for every k ∈ N, S k , X k and ε k are independent; (A2) the random vectors Y k , k ∈ N, are independent; (A3) the probabilities of presence P({ε k = 1}), k ∈ N, are less than or equal to one half and the random variables ε k , k ∈ N, are independent;

(A4) there exists some ν ∈ (0, ∞] such that S ∈ ∞ (N, L ν (Ω, R d )).
Given two non-negative real numbers r and s such that 0 ≤ s < r ≤ ν/2, some natural number m and any pair (σ, T ) of positive real numbers, define the random variable ∆ m (σ, T ) by

∆ m (σ, T )= m k=1 Y k r I( Y k ≤ σT ) m k=1 Y k s I( Y k ≤ σT ) -σ r-s Υ r (T )
Υ s (T ) .

(4) Then, σ 0 is the unique positive real number σ such that, for every

β 0 ∈ (0, 1], lim a(S)→∞ lim m ∆ m (σ, βξ(a(S)/σ)) ∞ = 0 (5) uniformly in β ∈ [β 0 , 1].
In the foregoing statement, U models a sequence of observations where, for every given k ∈ N, S k stands for some possible random signal and ε k models the possible occurrence of S k in the background of AWGN modelled by X.

In proposition 3.1, neither the probabilities of presence nor the signals are required to be identically distributed. Moreover, it is worth mentioning that the convergence criterion exhibited by this result derives from a corollary of Kolmogorov's classical strong limit theorem and not from usual generalizations of the central limit theorem such as the Lindeberg and the Lyapounov theorems.

Note that, for τ ∈ [0, ∞), the ratio

m k=1 Y k r I( Y k ≤ τ )/ m k=1 Y k s I( Y k ≤ τ )
is defined everywhere for the following reason. Let

x
m k=1 |x k | r / m k=1 |x k | s = 0 if (x 1 , . . . , x m ) = (0, . . . , 0).
Proposition 3.1 concerns positive solutions of Eq. ( 5) only because σ = 0 trivially satisfies this equation regardless of the specific convergence involved. Straightforwardly, Eq. ( 5) is also satisfied for all σ ∈ [0, ∞) when r = s ≥ 0. This explains why it is assumed that r > s ≥ 0.

IV. THE ESSENTIAL SUPREMUM, THE MODIFIED AND

THE COMPLEX ESTIMATES In this section, the notations used so far are kept with exactly the same meaning. On the basis of the previous theoretical result, we present several algorithms for the estimation of the noise standard deviation. We begin with an algorithm already introduced in [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF]. This algorithm requires prior knowledge of the minimum amplitude of the signals. We then propose a method that does not need such prior knowledge. This method is mainly heuristic. We particularize it to the case of two-dimensional real random vectors, or equivalently, complex random variables. This particular case will serve to address the applications treated in sections V and VI.

A. The Essential Supremum Estimate

According to Eq. ( 5), given β 0 ∈ (0, 1] and some positive real value η, there exists α 0 ∈ (0, ∞) that satisfies the following property: for every real number α larger than or equal to α 0 , every element S of ∞ (N, L ν (Ω, R d )) whose minimum amplitude is larger than or equal to α, and every β ∈ [β 0 , 1], the absolute value of the difference between the random variable

m k=1 Y k r I( Y k ≤ σβξ(α/σ)) m k=1 Y k s I( Y k ≤ σβξ(α/σ)) and σ r-s Υ r (σβξ(α/σ)) Υ s (σβξ(α/σ))
does not exceed η (a-s) when the sample size or number of observations m is large enough. In many applications, the signals S k , k ∈ N, have "finite energy" in the sense that the second-order moments

E[ S k 2 ], k ∈ N, are finite.
In what follows, this "finite energy" hypothesis is made. Therefore, with respect to assumption (A4), we consider the case ν = 2 and assume that S ∈ ∞ (N, L 2 (Ω, R d )). Since the values r and s must verify the inequalities 0 ≤ s < r ≤ 1, our rather natural choice is r = 1 and s = 0.

Given m observations Y 1 , . . . , Y m for which the minimum amplitude a(S) is known to be larger than or equal to some known α, if we set L ∈ N and β = /L for every ∈ {1, . . . , L}, the foregoing suggests estimating σ 0 by a possibly local minimum of

sup ∈{1,...,L} m k=1 Y k I( Y k ≤ β σξ(α/σ)) m k=1 I( Y k ≤ β σξ(α/σ)) -σΘ(β ξ(α/σ)) (6) where Θ is defined for every x ∈ [0, ∞) by Θ(x) = Υ 1 (x) Υ 0 (x) = x 0 t d exp(-t 2 /2)dt x 0 t d-1 exp(-t 2 /2)dt . ( 7
)
Following the terminology proposed in [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF], the estimate obtained by minimizing Eq. ( 6) is an Essential Supremum Estimate (ESE) of the noise standard deviation. This name follows from the fact that the essential supremum norm plays an important role in proposition 3.1 and its generalization stated in [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF]Theorem 1].

In [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF], experimental results are given for the case of independent signals that are two-dimensional random vectors uniformly distributed on a circle centred at the origin with known radius α and that have their probabilities of presence less than or equal to one half. Such random signals model modulated sinusoidal carriers whose amplitudes equal α and whose phases are uniformly distributed in [0, 2π]. The empirical bias and the empirical Mean-Square Error (MSE) of the ESE remain reasonably good for α ∈ {0, 0.25, 0.5, 0.75, . . ., 5} and any probability of presence in {0.1, 0.2, 0.3, 0.4, 0.5}. One conclusion of [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF] is therefore that the asymptotic conditions on which [3, Theorem 1] and proposition 3.1 rely, are not so constraining and can certainly be relaxed in practice.

B. Modified and Complex Essential Supremum Estimates

The computation of the ESE requires prior knowledge of a lower bound for the amplitudes of the signals. This lower bound is not always known. Since [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF] suggests that the convergence stated by ( 5) is quite fast, we hereafter set α to 0 in (6) and, thus, since ξ(0) = √ d (see [START_REF] Pastor | A sharp upperbound for the probability of error of the likelihood ratio test for detecting signals in white gaussian noise[END_REF]), estimate σ 0 by a possibly local minimum σ 0 of sup

∈{1,...,L} ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ m k=1 Y k I( Y k ≤ β σ √ d) m k=1 I( Y k ≤ β σ √ d) -σΘ(β √ d) ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ .
Any minimization routine for scalar bounded non-linear functions is suitable. For instance, the experimental results presented in the next sections were obtained with the MATLAB routine fminbnd.m based on parabolic interpolation (see [START_REF] Press | Numerical recipes in C, The Art of Scientific Computing[END_REF]). This routine is also that employed in [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF]

m 4(m/2-1) 2 . A typical choice is Q = 0.95, provided that m ≥ 24. Set h = 1/ 4m(1 -Q) and k min = m/2 -hm. The left endpoint is then σ min = Y [kmin] / √ d.
The reader is asked to refer to [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF] for justifications regarding the construction of this search interval.

Regarding the applications considered in the sequel, simulations of the same type as those proposed below show that σ 0 is a reasonably good estimate of σ 0 . However, we can propose a better estimate of σ 0 . This new estimate, henceforth called the Modified ESE and denoted by σ 0 , is computed on the basis of σ 0 by setting

σ 0 = γ m k=1 Y k 2 I( Y k ≤ σ 0 √ d) m k=1 I( Y k ≤ σ 0 √ d) (8)
where γ is some constant chosen empirically with respect to the application. For instance, this constant is different for the two applications addressed in the sequel. The rationale is the following. Under the assumptions of proposition 3.1 and in the particular case where β = 1, Eq. ( 5) means that, when the amplitudes of the signals are larger than or equal to some sufficiently large value α and the sample size is large enough, the random variable

m k=1 Y k r I( Y k ≤ σξ(α/σ 0 )) m k=1 Y k s I( Y k ≤ σξ(α/σ 0 )) approximates σ r-s Υ r (ξ(α/σ 0 ))/Υ s (ξ(α/σ 0 )).
For the computation of σ 0 , we have chosen r = 1, s = 0, and considered that the convergence stated by Eq. ( 5) holds true even for small values of α.

Assume now that S ∈ ∞ (N, L 4 (Ω, R d )), which remains reasonable although stronger than the "finite energy" assumption since

∞ (N, L 4 (Ω, R d )) ⊂ ∞ (N, L 2 (Ω, R d ))
. Then, we can choose r = 2 and s = 0. Still assuming that Eq. ( 5) is valid whatever the value of α, we consider that σ 2 0 can be approximated, in a certain mathematical sense that remains to be specified, by the random variable

Υ 0 (ξ(α/σ 0 )) Υ 2 (ξ(α/σ 0 )) × m k=1 Y k 2 I( Y k ≤ σξ(α/σ 0 )) m k=1 I( Y k ≤ σξ(α/σ 0 ))
By setting α = 0 and taking into account that ξ(0) = √ d, the foregoing leads to approximate σ 0 as proposed in Eq. ( 8)

with γ = Υ 0 ( √ d)/Υ 2 ( √ d). When d = 2
, the case addressed in the sequel, this constant equals 1.0937. According to the experimental results presented below, it may be necessary to adjust this value with respect to the application.

In practice, observations are often complex random variables, or equivalently, two-dimensional real random vectors. For instance, such observations can be the complex values provided by the standard I and Q decomposition encountered in most receivers in radar, sonar and telecommunication systems. Complex observations can also be simply the outcome of a Discrete Fourier Transform (DFT) as in some radar and sonar systems but also in a great variety of other applications. In section VI, we will consider an application to speech processing where the DFT plays a crucial role.

Because of the important role played by twodimensional real random vectors and complex random variables in many applications, the Modified ESE when d = 2 will be called the Complex ESE (C-ESE). In the two-dimensional case, note that the expression of Θ simplifies so that

Θ(x) = x 0 t 2 exp(-t 2 /2)dt 1 -exp(-x 2 /2) . (9)
V. APPLICATION TO THE NON-COHERENT DETECTION OF MODULATED SINUSOIDAL CARRIERS IN AWGN The bias, the consistency and the minimum attainable mean square estimation error of the Modified and Complex ESEs are still unknown. However, we can undertake some experiments making it possible to assess the performance of such estimates. The experiments presented in this section are motivated by the following facts.

We keep the notations used so far. Given any real number x, by thresholding test with threshold height h ∈ R, we hereafter mean the test T h = I( • ≥ h) that assigns 1 to y ∈ R d if y ≥ h and 0, otherwise.

Given any non-negative real number α, let T be the statistical test

I( • ≥ σ 0 ξ(α/ σ 0 )). Given k ∈ N, the decision of this test is that ε k is 1 if Y k ≥ σ 0 ξ(α/ σ 0 ) and that ε k is 0, otherwise.
If the Modified ESE σ 0 is a reasonably good estimate of the noise standard deviation, the performance of test T can be expected to approach that of the thresholding test T σ0ξ(α/σ0) with threshold height σ 0 ξ(α/σ 0 ).

To detect the presence of any signal with norm larger than or equal to α and prior less than or equal to one half, it follows from [5, Theorem VII.1] that the probability of error P e {T σ0ξ(α/σ0) } of T σ0ξ(α/σ0) satisfies the following inequalities

P e {L} ≤ P e {T σ0ξ(α/σ0) } ≤ V (α/σ 0 ). (10)
In these inequalities, P e {L} stands for the probability of error of the Minimum-Probability-of-Error (MPE) decision scheme L, that is the likelihood ratio test with the smallest possible probability of error among all possible hypothesis binary tests. In Eq. ( 10), the map V is defined for every x ∈ [0, ∞) and the reader can refer to [START_REF] Pastor | A sharp upperbound for the probability of error of the likelihood ratio test for detecting signals in white gaussian noise[END_REF] for the general expression of this map.

The inequalities in [START_REF] Quackenbush | Objective Measures of Speech Quality[END_REF] become equalities in the least favourable situation where the signal is uniformly distributed on the sphere centred at the origin with radius α and has prior equal to one half (see [START_REF] Pastor | A sharp upperbound for the probability of error of the likelihood ratio test for detecting signals in white gaussian noise[END_REF]Theorem VII.1]). Therefore, in this least favourable case, if the ESE σ 0 is a good estimate of σ 0 , the probability of error P e { T } of T should not significantly exceed V (α/σ 0 ). If the probability of presence of every S k equals one half, P e { T } should even be close to V (α/σ 0 ).

We do not know the theoretical value of P e { T }. Hence, we estimate this probability of error by the Binary Error Rate (BER) obtained by a Monte-Carlo simulation and compare this BER to V (α/σ 0 ). This Monte-Carlo simulation is carried out in the case of two-dimensional real random observations (d = 2). We then choose some p ∈ (0, 1/2] and every S k has a probability of presence equal to p and is uniformly distributed on the circle centred at the origin with radius α. Therefore, the two components of every S k can be regarded as the in-phase and quadrature components of a sinusoidal carrier. In other words, we consider the "non-coherent detection of modulated sinusoidal carriers", a problem particularly relevant for telecommunication and radar processing (see [7, p. 65]).

In the two-dimensional case, the mathematical expressions of ξ and V simplify. According to [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF], for d = 2, we have

ξ(ρ) = (1/ρ)I -1 0 (exp(ρ 2 /2)) (11)
and

V (ρ) = 1 2 exp(-ρ 2 /2) ξ(ρ) 0 t exp(-t 2 /2)I 0 (ρt)dt + 1 2 exp(-ξ(ρ) 2 /2), (12)
where I 0 is the zeroth-order modified Bessel function of the first kind (see [8, Eq. 9.6.47, p. 377]).

For the computation of the C-ESE and on the basis of some preliminary tests such as those described below, the constant γ in Eq. ( 8) is set to 1, the integer part of the constant suggested in the preceding section.

The BER of T is then computed as follows. Independent trials of m observations each are carried out until two conditions are fulfilled. First, at least M trials must be performed. Inasmuch as the decision about the presence or the absence of signals is made on the observations used for estimating σ 0 , the accuracy of the estimate affects m decisions at one go. This effect is then reduced by fixing a minimum number of trials. Second, trials are performed until the total number N e of errors made by test T for detecting the presence or the absence of signals is above or equal to some specified number N . If j is the first trial number larger than or equal to M for which the total number of errors N e becomes larger than or equal to N , the BER of test T is then defined as the ratio N e /(j ×m).

The simulation is achieved with σ 0 = 1. The prespecified number of errors is N = 1000 and the minimum number of trials is M = 1000. We choose L = m and Q = 0.95 on the basis of preliminary trials. The comparison between the BER of T and V (α/σ 0 ) is achieved for p ∈ {0.1, 0.2, 0.3, 0.4, 0.5} and A ∈ {0, 0.25, 0.5, 0.75, . . ., 5}. The results are those of figures 1, 2, 3 and 4 for different values of m. As expected, the performance of T yields performance close to that of T σ0ξ(α/σ0) . Again for signals that have the same probability of presence p less than or equal to one half and that are uniformly distributed on the circle centred at the origin with radius α, we calculated the absolute value |BIAS| and the empirical MSE of the C-ESE computed over M = 1000 trials. Tables I, II and III below display some results obtained for different values of α and p.

The experimental results presented above were achieved without prior knowledge on the amplitude of the signals. Nevertheless, they are significantly close to those given in [START_REF] Pastor | A theoretical result for processing signals that have unknown distributions and priors in white gaussian noise[END_REF] where the estimation of the noise standard deviation is performed via the ESE, which does require such prior knowledge.

VI. APPLICATION TO SPEECH DENOISING

In this section, we describe how the C-ESE can be used by the standard Wiener filter aimed at denoising speech signals corrupted by independent AWGN with unknown standard deviation. We briefly remind the reader of the basics concerning the standard Wiener filtering of noisy speech signals (see subsection VI-A). Then, in subsection VI-B, under the assumption that noise is white and Gaussian, we explain how to use the C-ESE to estimate the noise standard deviation and adjust the filtering. We complete this section by experimental results (see subsection VI-C).

A. Wiener filtering

Let s(t), t = 0, . . . , T -1, be the samples of some speech signal and suppose that these T samples are corrupted by some noise x(t), t = 0, 1, . . . , T -1, so that the samples of the observed signal are

y(t) = s(t) + x(t), t = 0, . . . , T -1. ( 13
)
It is usual to split the T available samples y(t), t = 0, 1, . . . , T -1, into frames of N = 2 q samples each where q is some integer such that NF s ≈ 20ms, F s being the sampling frequency. Frames are generally constructed so that two consecutive ones overlap by one half. The samples of each frame can be weighted. In what follows, for the sake of simplifying the notations, the description of the filtering is made without taking such weighting into account.

For the k th frame, let s k (t), x k (t) and y k (t), t = 0, 1, . . . , N -1, stand for the N samples of the speech signal, noise and the observed noisy speech signal, respectively. We thus have

y k (t) = s k (t) + x k (t).
Now, let Y k (ν), S k (ν) and X k (ν), ν = 0, . . . , N -1, denote the DFT coefficients of y k (t), s k (t) and x k (t), t = 0, 1, . . . , N -1, respectively. For every ν = 0, 1, . . . , N -

1, we have Y k (ν) = S k (ν) + X k (ν).
The standard Wiener filtering of y is based on Malah's decision-directed approach (see [START_REF] Ephraim | Speech enhancement using a minimum mean square error short-time spectral amplitude estimator[END_REF]). Each frequency com-

ponent S k (ν) is then estimated by S k (ν) = W k (ν)Y k (ν)
where W k (ν) is the so-called Wiener gain function. This gain is given by

W k (ν) = R k (ν)/(1 + R k (ν)) (14)
where

R k (ν) = (1 -w)h (ζ k (ν) -1) + w | S k-1 (ν)| 2 E [|X k (ν)| 2 ] ( 15
)
is the so-called decision-directed estimate of the a priori Signal to Noise Ratio (SNR)

E |S k (ν)| 2 /E |X k (ν)| 2 . In Eq. (15), S k-1 (ν) = W k-1 (ν)Y k-1 (ν) is the ν th spectral component of the Wiener denoised speech signal in frame k-1; h(x) = x if x ≥ 0 and h(x) = 0 otherwise; ζ k (ν) = |Y k (ν)| 2 /E |X k (ν)| 2
is the a posteriori SNR; the weighting factor w will be set to 0.98 as recommended in [START_REF] Ephraim | Speech enhancement using a minimum mean square error short-time spectral amplitude estimator[END_REF]. The decision-directed approach described by (15) takes into account the current frame, with weight (1w), and the result of the processing of the previous frame, with weight w.

B. Estimation of the noise standard deviation via the C-ESE when noise is white and Gaussian

When noise is white and Gaussian with standard deviation σ 0 , the noise frequency components are independent and identically distributed so that each X k (ν) is complex Gaussian distributed with

X k (ν) ∼ CN (0, 2NA 2 σ 2 0 ) ( 16
)
where A is some constant depending on the normalization chosen with regard to the Parseval identity. More specifically, if the sequence U (ν), ν = 0, 1, . . . , N -1, is the DFT of the discrete sequence of (possibly complex) values u(t), t = 0, 1, . . . , N -1, this constant is such that

N -1 ν=0 |U (ν)| 2 = NA 2 N -1 t=0 |u(t)| 2 and generally chosen in {1/ √ N, 1, 1/N }.
To perform the standard Wiener filtering according to Eq. ( 15), an estimate of the noise standard deviation must be calculated so as to estimate

E[|X k (ν)| 2 ] = 2NA 2 σ 2 0 .
A basic and popular solution to perform this estimate involves using a Voice Activity Detector (VAD). The estimate of σ 0 is then computed on the basis of the samples of the time frames that the VAD has detected as noise alone. Subspace approaches can also be used to estimate σ 0 by computing the smallest eigenvalues of the noisy speech autocorrelation matrix. In this case, the model order is difficult to choose and the computation of the eigenvalues may prove unstable.

We now present how to use the C-ESE for the estimation of σ 0 . We begin by splitting the T available samples y(t), t = 0, 1, . . . , T -1, into frames of N = 2 q successive samples each. In contrast to the previous subsection, the frames are here constructed so that they do not intersect; it is even conceivable to choose frames separated by some time lapse. We proceed thus for reasons given below.

Let K stand for the number of frames so constructed. The notations used above to denote the observed signal, the speech signal and noise in the k th frame are again y k (t), s k (t) and x k (t), t = 0, 1, . . . , N -1, respectively. Consider now the set of observations Y k (ν), k ∈ {1, . . . , K}, ν ∈ {0, . . . , N -1}, where, again as above, Y k (ν), S k (ν) and X k (ν), ν = 0, 1, . . . , N -1, stand for the respective DFTs of the signals y k (t), s k (t) and x k (t), t = 0, 1, . . . , N -1. Because of the Hermitian symmetry of the DFT, we can restrict our attention to half of this set, namely the complex random variables Y k (ν), k ∈ {1, . . . , K}, ν ∈ {0, . . . , N/2 -1} where N is assumed to be even.

Of course, we still have Y k (ν) = S k (ν) + X k (ν) for every k = 1, 2, . . . , K and every ν = 0, 1, . . . , N -1. However, depending on the Signal to Noise Ratio (SNR) and the type of speech signal present during frame k, some speech time-frequency components can be neglected in comparison to noise and other speech time-frequency components. For instance, high frequency components of voiced speech signals can become negligible with regard to noise and the low-frequency components of the same speech signals. For unvoiced fricative speech signals, lowfrequency components can also be significantly smaller than high-frequency ones and smaller than those due to noise.

According to these remarks, we propound to model the presence and the absence of the speech timefrequency component S k (ν) by a discrete random variable ε k (ν) valued in {0, 1} and write that the observation is

Y k (ν) = ε k (ν)S k (ν) + X k (ν).
With respect to this model, P({ε k (ν) = 1}) is the probability that some speech component is present at frequency ν during frame k. This probability of presence may be larger than one half for low frequency components. However, for high frequency components, this probability of presence becomes less than or equal to 1/2 and even relatively small. Summarizing, we can consider that the speech timefrequency components are less present than absent. If we now split the observation set Y k (ν), k ∈ {1, . . . , K}, ν ∈ {0, . . . , N/2-1}, into subsets of m observations each and if the elements of each subset are chosen randomly amongst the available observations, the observations of each subset can reasonably be expected to be independent all the more since these observations are computed on the basis of non-intersecting frames. We thus assume that the model underlying proposition 3.1 and the construction of the C-ESE are satisfied for each subset of m observations. According to Eq. ( 16), for each subset of m observations, the C-ESE will return an estimate of Aσ 0 √ N since each complex random vector X k (ν) can be regarded as a centred two-dimensional random vector whose covariance matrix equals Aσ 0 √ N I 2 . Therefore, our final estimate of Aσ 0 √ N is simply obtained by averaging all these estimates returned by the C-ESE. Dividing this average by A √ N yields an estimate of σ 0 .

C. Performance evaluation

We randomly select twenty-five sentences in the TIMIT database and twenty-five sentences in the TIDIG-ITS database. These two databases are standard in speech processing. The TIMIT and the TIDIGITS databases are composed of read speech and connected digit sequences, respectively. The sentences selected in these databases are downsampled to 8 kHz before adding white Gaussian noise. We estimate the noise standard deviation as described above via the C-ESE on the basis of frames with N = 256 samples each. Such a frame thus corresponds to 32ms of noisy speech signals. As specified above, we use non-overlapping frames for the estimation of the noise standard deviation. They are not weighted either. On the basis of preliminary tests on signals other than the sentences used to achieve the experimental results presented below, we set the constant γ in (8) to √ 2. As far as the Wiener filtering is concerned, there is a 50% overlap between two adjacent frames and each frame is weighted by a Hanning window.

We evaluate the quality of the filtered speech signals by calculating the standard Segmental Signal to Noise Ratio (SSNR) (see [START_REF] Quackenbush | Objective Measures of Speech Quality[END_REF]) and the Modified Bark Spectral Distortion (MBSD) (see [START_REF] Yang | Modified bark spectral distortion measure which uses noise masking threshold[END_REF]). The SSNR is the average of the SNR values on short segments. The SSNR is not relevant enough to measure the distortion of the denoised speech signals. This is the reason why we use the MBSD. The MBSD is an improved version of the Bark Spectral Distortion (BSD) proposed in [START_REF] Wang | An objective measure for predicting subjective quality of speech coders[END_REF]. Basically, it extends the BSD by incorporating the speech masking threshold so as to perform the measure without taking into account inaudible distortions. The MBSD proves to be highly correlated with subjective speech quality assessment [START_REF] Yang | Modified bark spectral distortion measure which uses noise masking threshold[END_REF].

The average SSNR and MBSD obtained over the twenty-five sentences of the TIMIT database are presented in figures 5(a) and 5(b). Figures 6(a) and 6(b) are the average SSNR and MBSD calculated over the twentyfive sentences of the TIDIGITS database. In all these figures, the solid curves are the performance measurements achieved with the filtering defined by equations ( 14) and (15) when the C-ESE is used instead of the true value of the noise standard deviation. The dashed curves are the results obtained when the filtering is achieved with the exact value of the noise standard deviation. The results are similar for the two sets of sentences considered during these experiments. Clearly, the Wiener filtering adjusted with the noise standard deviation estimate yields results that are significantly close to those obtained when the noise standard deviation is known. In addition, figure 7(a) (resp. figure 7(b)) compares the average value of the C-ESE obtained by processing the 25 sentences of the TIMIT database (resp. of the TIDIGITS database) to the true value of the noise standard deviations at the different SNR's tested during the experiments. These figures illustrate the good behaviour of our algorithm. As mentioned below, they are preliminary to further studies undertaken to get better insight into the behaviour of the C-ESE with respect to applications in speech denoising.

VII. PERSPECTIVES AND EXTENSIONS

We have presented algorithms for estimating the standard deviation of some AWGN when observations derive from signals less present than absent in this background. According to experimental results, this algorithm is very promising. An application to speech denoising has been described where the estimator avoids the use of a DAV or subspace approaches that can be difficult to adjust.

New theoretical developments should be achieved in forthcoming work. In particular, further study of the Modified and Complex ESEs are required to get better insight into the behaviour of these estimates, in particular with regard to the constant γ in Eq. ( 8) since this constant is chosen empirically depending on the application.

More exhaustive experiments are in progress so as to better characterize the statistical behaviour of the C-ESE in the case of noisy speech signals. The analysis of the bias, the consistency and the mean square estimation error of the estimators proposed in this paper should also be addressed from a theoretical point of view.

As an extension of the application to speech denoising described above, on-going work concerns the application of the C-ESE to perceptual speech denoising in AWGN and coloured noise. As far as coloured noise is concerned, the use of the C-ESE could be extended along the following lines. The processing would involve assuming that the coloured noise power spectral density (psd) is a step function. Given a set of DFT bins where the noise psd is assumed to be constant, this constant value could then be estimated by computing the C-ESE on the basis of the values returned in these frequency bins by the DFTs of non-intersecting time-frames. By using nonintersecting frames, we still expect that the time-frequency components can reasonably be considered as independent.

The design of CFAR systems for the detection of radar targets also seems to be a rather natural field of application. Application in ESM for the interception of noncooperant communications could also be investigated. Two other areas of applications are proximity sensing and distributed detection systems (see [START_REF] Varshney | Distributed Detection and Data Fusion[END_REF]). Proximity sensing aims at informing a robot that it is approaching an object. For the design of distributed detection systems, the absence of prior knowledge about the statistics of the observations is still an issue.

  This paper is an extended version of "Estimating the Standard Deviation of Some Additive White Gaussian Noise on the Basis of Non Signal-Free Observations" by D. Pastor, which appeared in the Proceedings of the 2006 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP 2006, Toulouse, France, May, 2006. c 2006 IEEE.

Figure 1 .

 1 Figure 1. V (A/σ 0 ) vs BERs of b T for m = 100

Figure 2 .Figure 3 .

 23 Figure 2. V (A/σ 0 ) vs BERs of b T for m = 200

Figure 4 .

 4 Figure 4. V (A/σ 0 ) vs BERs of b T for m = 800

Figure 5 .

 5 Figure 5. (a) SSNR and (b) MBSD improvement for 25 speech signals randomly chosen in the TIMIT database and additively corrupted by independent AWGN with various SNRs.

Figure 6 .

 6 Figure 6. (a) SSNR and (b) MBSD improvement for 25 speech signals randomly chosen in the TIDIGITS database and additively corrupted by independent AWGN with various SNRs.

  of the noise standard deviation True value of the noise standard deviation (a) Average over 25 sentences of the TIMIT database

Figure 7 .

 7 Figure 7. For each tested SNR, comparison between the true value of the noise standard deviation and the average value of the C-ESE computed over (a) the 25 noisy sentences of the TIMIT database and (b) the 25 noisy sentences of the TIDIGITS database.

  1 , x 2 , . . . , x m be m real numbers. If there exists at least one natural number k ∈ {1, . . . , m} such that x k = 0, the finiteness of the ratio m k=1 |x k | r / m k=1 |x k | s is trivial. Since r > s, the definition of this ratio is then extended by continuity by setting

  to compute the ESE. The search interval [σ min , σ max ] is computed as follows. Sort the observations Y 1 , . . . , Y m , k = 1, . . . , m, by increasing norm. Let Y [k] , k = 1, . . . , m, be the resulting sequence. The right endpoint of the search interval is thenσ max = U [m] / √ d.As far as the left endpoint is concerned, choose a real number Q close to 1 but less than or equal to 1 -

TABLE I .

 I ABSOLUTE VALUE OF THE EMPIRICAL BIAS AND EMPIRICAL MSE OF THE C-ESE WHEN SIGNALS ARE UNIFORMLY DISTRIBUTED ON THE CIRCLE CENTRED AT THE ORIGIN WITH RADIUS α AND HAVE A PROBABILITY OF PRESENCE p EQUAL TO 0.1.

	α	0.5	1	1.5	2	2.5	3	3.5	4
	m = 100	|BIAS| 0.077 MSE 0.036	0.07 0.039	0.053 0.037	0.039 0.032	0.04 0.035	0.049 0.036	0.048 0.034	0.064 0.04
	m = 200	|BIAS| 0.065 MSE 0.023	0.051 0.023	0.03 0.02	0.041 0.022	0.054 0.026	0.052 0.024	0.054 0.022	0.056 0.023
	m = 400	|BIAS| 0.069 MSE 0.016	0.058 0.017	0.05 0.016	0.045 0.015	0.054 0.015	0.061 0.016	0.068 0.017	0.062 0.016
	m = 800	|BIAS| 0.082 MSE 0.014	0.07 0.013	0.062 0.012	0.054 0.011	0.059 0.011	0.067 0.012	0.07 0.012	0.074 0.012

TABLE II .

 II ABSOLUTE VALUE OF THE EMPIRICAL BIAS AND EMPIRICAL MSE OF THE C-ESE WHEN SIGNALS ARE UNIFORMLY DISTRIBUTED ON THE CIRCLE CENTRED AT THE ORIGIN WITH RADIUS α AND HAVE A PROBABILITY OF PRESENCE p EQUAL TO 0.3.

	α		0.5	1	1.5	2	2.5	3	3.5	4
	m = 100	|BIAS| 0.071 MSE 0.037	0.024 0.038	0.02 0.042	0.041 0.046	0.039 0.05	0.014 0.041	0.005 0.038	0.01 0.035
	m = 200	|BIAS| 0.059 MSE 0.024	0.005 0.023	0.028 0.024	0.043 0.028	0.033 0.028	0.013 0.021	0.002 0.018	0.017 0.018
	m = 400	|BIAS| 0.065 MSE 0.017	0.017 0.015	0.026 0.016	0.026 0.017	0.02 0.014	0.002 0.012	0.012 0.01	0.014 0.01
	m = 800	|BIAS| MSE	0.07 0.013	0.033 0.01	0.009 0.01	0.024 0.01	0.002 0.008	0.009 0.006	0.017 0.005	0.023 0.006

TABLE III .

 III ABSOLUTE VALUE OF THE EMPIRICAL BIAS AND EMPIRICAL MSE OF THE C-ESE WHEN SIGNALS ARE UNIFORMLY DISTRIBUTED ON THE CIRCLE CENTRED AT THE ORIGIN WITH RADIUS α AND HAVE A PROBABILITY OF PRESENCE p EQUAL TO 0.5.

	α	0.5	1	1.5	2	2.5	3	3.5	4
	m = 100	|BIAS| 0.056 MSE 0.035	0.02 0.04	0.124 0.068	0.193 0.098	0.223 0.128	0.161 0.099	0.101 0.059	0.07 0.041
	m = 200	|BIAS| 0.038 MSE 0.021	0.042 0.027	0.145 0.052	0.207 0.083	0.209 0.088	0.15 0.057	0.095 0.029	0.08 0.021
	m = 400	|BIAS| 0.045 MSE 0.014	0.028 0.016	0.129 0.036	0.197 0.063	0.177 0.054	0.135 0.031	0.113 0.02	0.105 0.017
	m = 800	|BIAS| 0.059 MSE 0.012	0.012 0.01	0.112 0.024	0.17 0.043	0.16 0.037	0.141 0.026	0.137 0.021	0.145 0.024


L Éo Girardin
email: leo.girardin@math.u-psud.fr.

TWO COMPONENTS IS TOO SIMPLE: AN EXAMPLE OF OSCILLATORY FISHER-KPP SYSTEM WITH THREE COMPONENTS

Keywords: 2000 Mathematics Subject Classification. 35K40, 35K57, 37G10, 92D25 KPP nonlinearity, reaction-diffusion system, Hopf bifurcation, spreading phenomena

In a recent paper by Cantrell, Cosner and Yu [9], two-component KPP systems with competition of Lotka-Volterra type were analyzed and their long-time behavior largely settled. In particular, the authors established that any constant positive steady state, if unique, is necessarily globally attractive. In the present paper, we give an explicit and biologically very natural example of oscillatory three-component system. Using elementary techniques or preestablished theorems, we show that it has a unique constant positive steady state with two-dimensional unstable manifold, a stable limit cycle, a predatorprey structure near the steady state, periodic wave trains and point-to-periodic rapid traveling waves. Numerically, we also show the existence of pulsating fronts and propagating terraces.

Introduction

This paper is concerned with the long-time properties of the following reactiondiffusion system:

(KPP µ )    ∂ t u 1 -∆u 1 = u 1 + µ (-2u 1 + u 2 + u 3 ) -1
10 (u 1 + 8u 2 + u 3 ) u 1 ∂ t u 2 -∆u 2 = u 2 + µ (+u 1 -2u 2 + u 3 ) - 1 10 (u 1 + u 2 + 8u 3 ) u 2 ∂ t u 3 -∆u 3 = u 3 + µ (+u 1 + u 2 -2u 3 ) - 1 10 (8u 1 + u 2 + u 3 ) u 3 , written in vector form as:

∂ t u -∆u = u + µMu -(Cu) • u,
where µ is a positive constant,

M =   -2 1 1 1 -2 1 1 1 -2   , C = 1 10   1 8 1 1 1 8 8 1 1   ,
and • denotes the Hadamard product (component-by-component product) of two vectors. This system (KPP µ ) is a particular case of KPP system, as defined by the author in [START_REF] Girardin | Non-cooperative Fisher-KPP systems: traveling waves and long-time behavior[END_REF] and subsequently analyzed in [START_REF] Girardin | Non-cooperative Fisher-KPP systems: Asymptotic behavior of traveling waves[END_REF]. Let us remind briefly here that this name comes from the fact that the above reaction term is strongly reminiscent of the scalar logistic term u (1 -u) and leads to very similar conclusions regarding extinction, persistence, traveling waves and spreading speed.

From the biological point of view, (KPP µ ) can for instance model a structured population with three coexisting phenotypes subjected to spatial dispersal, phenotypical changes and competition for resources. As explained by Cantrell, Cosner and Yu [START_REF] Stephen Cantrell | Dynamics of populations with individual variation in dispersal on bounded domains[END_REF], the phenotypical changes can come from behavioral switching, phenotypic plasticity or Darwinian evolution, for instance.

The paper is organized as follows. In the remaining of this introductive section, we present and comment our results. In Section 2, we prove our main analytical result. In Section 3, we present our numerical findings. In the appendix, we give elementary proofs of related new results on periodic wave trains (Appendix A) or KPP systems (Appendix B and Appendix C).

1.1. Main result. In their paper, Cantrell, Cosner and Yu [START_REF] Stephen Cantrell | Dynamics of populations with individual variation in dispersal on bounded domains[END_REF] studied the KPP system with competition of Lotka-Volterra type

∂ t u -diag (d) ∆u = Lu -Cu • u
with only two components but in full generality with respect to the parameters. Here we recall that the minimal KPP assumptions are the positivity 1 of d and C as well as the essential nonnegativity 2 of L, its irreducibility and the positivity of its Perron-Frobenius eigenvalue. We also recall that a competition term c (u) • u is referred to as a Lotka-Volterra competition term if the vector field c is linear, namely c (u) = Cu.

Cantrell, Cosner and Yu obtained an almost complete characterization of the long-time asymptotics in bounded domains. In particular, they proved that any constant positive steady state, if unique, is globally attractive under Neumann boundary conditions. Hence it is an important step forward regarding the general study of two-component KPP systems with Lotka-Volterra competition, which have been studied by several authors in the past few years [START_REF] Elliott | Dispersal polymorphism and the speed of biological invasions[END_REF][START_REF] Girardin | Non-cooperative Fisher-KPP systems: Asymptotic behavior of traveling waves[END_REF][START_REF] Griette | Existence and qualitative properties of travelling waves for an epidemiological model with mutations[END_REF][START_REF] Morris | Individual variability in dispersal and invasion speed[END_REF]. We will show briefly in Proposition B.1 how their results and arguments of proof can be applied to the unbounded setting.

On the contrary, for the three-component system (KPP µ ), a qualitatively completely different result will be proved in the forthcoming pages. Before giving the statement, we point out that the parameters of the system are obviously normalized in such a way that, for any value of µ, 1 = (1, 1, 1)

T is a positive constant steady state. Additionally, we define µ H = 7 60 , µ -= 1 10 and µ + = 8 10 . These values satisfy 0 < µ -< µ H < µ + < 1.

In the statement below, (e i ) i∈{1,2,3} denotes the canonical basis of R 3 . For notation convenience, indices are defined modulo 3 (i.e. e 4 = e 1 , e 5 = e 2 , etc.).

Theorem 1.1. The diffusionless system

(KPP 0 µ ) u = u + µMu -(Cu) • u
satisfies the following properties.

(1) 1 is the unique positive steady state.

(2) If µ > µ H , 1 is locally asymptotically stable, but at µ = µ H , it undergoes a supercritical Hopf bifurcation leading to the birth of a unique and locally asymptotically stable limit cycle C µ . Using µ as a parameter, there exists a family of positive limit cycles (C µ ) µ∈(0,µH) and any such family converges, 1 In the whole paper, positive vectors are vectors with positive components and nonnegative, negative and nonpositive vectors are defined analogously. 2 Off-diagonal nonnegativity.

in the sense of the Hausdorff distance, as µ → 0, to

C 0 = i∈{1,2,3} {10e i } ∪ H i ,
where H i is, for (KPP 0 µ ) with µ = 0, the unique heteroclinic connection between 10e i and 10e i+1 , which lies in e ⊥ i+2 . Furthermore, any limit cycle

C µ satisfies C µ ⊂ v ∈ (0, +∞) 3 | 1 ≤ v 1 + v 2 + v 3 3 ≤ 10 3
and is rotating clockwise around span (1) if seen from 0.

(3) Let v ≥ 0. The reaction term at v is cooperative if and only if v ∈ 0, µ µ + 3
and competitive if and only if

v ∈ µ µ - , +∞ 3 . If v ∈ µ µ+ , µ µ-

3

, the off-diagonal entries of the linearized reaction term at v have the following signs:

  • -+ + • - -+ •   .
This is in particular the case if v = 1 and µ ∈ (µ -, µ + ).

We will illustrate numerically that the limit cycle C µ seems to be in fact globally attractive (with respect to initial conditions that are not in the basin of attraction of 0 or 1, namely almost all of them), and therefore also unique. However, we did not manage to prove the global attractivity or the uniqueness.

1.2. Discussion on Theorem 1.1. Although the third property follows from a direct differentiation, it is qualitatively very meaningful. On one hand, in the cube 10µ 8 , 10µ 3 , the system has the structure of a cyclic predator-prey system (rockpaper-scissor-like). On the other hand, a consequence of Cantrell-Cosner-Yu [9, Propositions 2.5 and 3.1] is that two-component KPP systems with Lotka-Volterra competition are competitive in the neighborhood of any constant positive saddle. The case µ ∈ (µ -, µ H ) of Theorem 1.1 above proves that this property fails with three components. For the sake of completeness, we will prove in Proposition C.1 what seems to be the optimal result for an arbitrary number of components: at an unstable constant positive steady state, the reaction term is not cooperative. Notice that, changing in (KPP 0 µ ) µ into 1 µ and normalizing appropriately the time variable, we obtain the system

∂ t u -Mu = µ (u -(Cu) • u) ,
which is obviously similar (by finite difference approximations and Riemann sum approximations, see [START_REF] Girardin | Non-cooperative Fisher-KPP systems: traveling waves and long-time behavior[END_REF]Section 1.5]) to the following nonlocal KPP equation:

∂ t u -∂ yy u = µu (1 -φ y u) .
For this equation, it is now well-known that as µ increases, the steady state 1 is dynamically destabilized (we refer for instance to Berestycki-Nadin-Perthame-Ryzhik [START_REF] Berestycki | The non-local Fisher-KPP equation: travelling waves and steady states[END_REF], Faye-Holzer [START_REF] Faye | Modulated traveling fronts for a nonlocal fisher-kpp equation: a dynamical systems approach[END_REF] and Nadin-Perthame-Tang [START_REF] Nadin | Can a traveling wave connect two unstable states?: The case of the nonlocal Fisher equation[END_REF]). In this context, this property is usually understood as a form of Turing instability. In particular, Nadin-Perthame-Tang [START_REF] Nadin | Can a traveling wave connect two unstable states?: The case of the nonlocal Fisher equation[END_REF] showed numerically how this Turing instability can lead to interesting spreading phenomena, where the classical traveling waves connecting 0 to 1 are replaced by more sophisticated solutions. 1.3. Discussion and numerical results on the spatial structure. Of course, the presence in (KPP µ ) of a third variable x makes the system (KPP µ ) qualitatively different from the nonlocal KPP equation. In fact, as explained by the author in [START_REF] Girardin | Non-cooperative Fisher-KPP systems: traveling waves and long-time behavior[END_REF], (KPP µ ) is more reminiscent of the cane toad equation with nonlocal competition and local or nonlocal mutations [START_REF] Alfaro | Travelling waves in a nonlocal reactiondiffusion equation as a model for a population structured by a space variable and a phenotypic trait[END_REF][START_REF] Arnold | Existence of nontrivial steady states for populations structured with respect to space and a continuous trait[END_REF][START_REF] Bénichou | Front acceleration by dynamic selection in fisher population waves[END_REF][START_REF] Bouin | Travelling waves for the cane toads equation with bounded traits[END_REF][START_REF] Bouin | Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration[END_REF][START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF][START_REF] Prévost | Applications of partial differential equations and their numerical simulations of population dynamics[END_REF]. For these equations, what happens in the wake of an invasion front is still poorly understood. Our results and numerical findings are, in this regard, quite interesting.

Taking profit of the Hopf bifurcation at µ = µ H , we can apply a theorem of Kopell-Howard [START_REF] Kopell | Plane wave solutions to reaction-diffusion equations[END_REF] (we also refer to Murray [

where γ ∈ [0, γ s ] ∪ [γ l , 1] is the parameter, 0 < γ s ≤ γ l < 1, κ γ ∈ R, σ γ ∈
R and p γ is positive and periodic. Without loss of generality, γ can be understood as an amplitude parameter:

• the image of p 0 is 1;

• κ 1 = 0 and the image and period σ -1

1
of p 1 are respectively the limit cycle C µ and its associated period; • γ → p γ is increasing in the sense that the image of

(γ, ξ) ∈ [0, γ 1 ] × R → p γ (ξ) is strictly included in that of (γ, ξ) ∈ [0, γ 2 ] × R → p γ (ξ) provided γ 1 < γ 2 .
Furthermore, there exists γ ∈ (0, γ s ] such that all wave trains of amplitude γ ∈ [0, γ) are unstable with respect to compactly supported, bounded perturbations.

Beware that the Kopell-Howard theorem only shows that there are wave trains close to 1 (small amplitude 0 ≤ γ ≤ γ s ) and close to C µ (large amplitude γ l ≤ γ ≤ 1). It is unclear, even numerically, whether a continuum of wave trains exists (i.e., the equality γ s = γ l is unclear).

The nonlinear stability of a wave train of amplitude close to 1 is a delicate question, as established by Kopell-Howard [START_REF] Kopell | Plane wave solutions to reaction-diffusion equations[END_REF] and subsequently confirmed by Maginu [START_REF] Maginu | Stability of spatially homogeneous periodic solutions of reaction-diffusion equations[END_REF][START_REF] Maginu | Stability of periodic travelling wave solutions with large spatial periods in reaction-diffusion systems[END_REF]. Nevertheless, simply thanks to the fact that the diffusion matrix is the identity, the stability of the limit cycle extends in the following way. Proposition 1.3. If µ < µ H and C µ is locally asymptotically stable with respect to (KPP 0 µ ), then there exists γ ∈ [γ l , 1) such that all wave trains of amplitude γ ∈ γ, 1 are marginally stable in linearized criterion.

The proof of Proposition 1.3 is very simple but is actually not provided in [START_REF] Maginu | Stability of periodic travelling wave solutions with large spatial periods in reaction-diffusion systems[END_REF]. For the sake of completeness, it will be detailed in Appendix A. Note that the notion of stability in the statement above is the marginal stability in linearized criterion [START_REF] Paul | Mathematical aspects of reacting and diffusing systems[END_REF] and not the asymptotic waveform stability [START_REF] Maginu | Stability of periodic travelling wave solutions with large spatial periods in reaction-diffusion systems[END_REF], which might fail in general and remains a difficult and open question.

Numerically, we will observe propagating terraces (succession of compatible waves with decreasingly ordered speeds, as defined by Ducrot, Giletti and Matano [START_REF] Ducrot | Existence and convergence to a propagating terrace in one-dimensional reaction-diffusion equations[END_REF]) where 0 is invaded by 1 and then 1 is slowly invaded by a stable wave train of amplitude γ close to 1. The former invasion takes the form of a traveling wave (here defined as an entire solution with constant profile and speed) whereas the latter takes the form of a pulsating front (more general entire solution connecting at some constant speed two periodic, possibly homogeneous, solutions, defined for instance by Nadin in [START_REF] Nadin | Traveling fronts in space-time periodic media[END_REF] and also known as pulsating traveling wave).

Regarding traveling waves, the following proposition can be straightforwardly established by looking for wave profiles of the form ξ → p (ξ) 1 (we refer to [START_REF] Girardin | Non-cooperative Fisher-KPP systems: Asymptotic behavior of traveling waves[END_REF] for a similar construction).

Proposition 1.4. The system (KPP µ ) admits a family of monotonic traveling plane waves connecting 0 to 1 at speed c ≥ 2. This is indeed such a monotonic profile we observe numerically. By direct application of theorems due to Fraile and Sabina [START_REF] José | Kinetic conditions for the existence of wave fronts in reactiondiffusion systems[END_REF][START_REF] José | General conditions for the existence of a "critical pointperiodic wave front" connection for reaction-diffusion systems[END_REF], there exist also point-to-periodic rapid traveling waves connecting 1 to wave trains of large amplitude; however, these are not the pulsating fronts that we observe numerically, which do not have a constant profile and have two distinct speeds, the one of the invasion front and the one of the wave train.

More precisely, close to the bifurcation value µ H , the speed of the pulsating front of the terrace connecting 1 to the wave train is 2 3 (µ H -µ). On the contrary, the intrinsic speed of the wave train, c γ = σγ κγ , is negative and of very large absolute value (consistently with |c γ | → ∞ as γ → 1). We emphasize that the preceding formula for the invasion speed is linearly determinate (in some sense precised below in Subsection 3.3) and was first predicted heuristically by Sherratt [START_REF] Sherratt | Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave[END_REF]. Interestingly, in Sherratt's predictions, both the invasion speed and the wave train speed do not depend on the initial condition or even on the speed of the first invasion (1 into 0). This is confirmed by our numerical experiments.

Since this pulsating front is parametrized by two distinct speeds (that of the invasion and that of the wave train), the interesting problem of its existence is in fact very difficult. Seemingly similar results on periodic wave trains [START_REF] Kopell | Plane wave solutions to reaction-diffusion equations[END_REF][START_REF] Murray | Spatial models and biomedical[END_REF], pointto-periodic rapid traveling waves [START_REF] José | Kinetic conditions for the existence of wave fronts in reactiondiffusion systems[END_REF][START_REF] José | General conditions for the existence of a "critical pointperiodic wave front" connection for reaction-diffusion systems[END_REF] or even space-periodic pulsating fronts (recently studied by Faye and Holzer [START_REF] Faye | Modulated traveling fronts for a nonlocal fisher-kpp equation: a dynamical systems approach[END_REF]) are proved by means of codimension 1 bifurcation arguments. The space-time periodic pulsating front at hand is a codimension 2 bifurcation problem. Its resolution is definitely outside the scope of this paper and we leave it for future work.

Another prediction of Sherratt [START_REF] Sherratt | Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave[END_REF] is the possible nonexistence of such propagating terraces when the speed c 1 of the first invasion satisfies

c 1 > 7 √ 3 20 3 (µ H -µ) = 7 20 √ µ H -µ .
More precisely, when this condition holds, a periodic wave train of speed c 1 and small amplitude exists and therefore there is the possibility of a point-to-periodic rapid traveling wave connecting directly 0 to this unstable wave train. Recall however that solutions that are initially compactly supported asymptotically spread at speed 2 [START_REF] Girardin | Non-cooperative Fisher-KPP systems: traveling waves and long-time behavior[END_REF], which is clearly smaller than the above threshold close to the bifurcation value, whence these traveling waves are irrelevant regarding biological applications. Anyway, even with initial conditions with appropriate exponential decay [START_REF] Girardin | Non-cooperative Fisher-KPP systems: Asymptotic behavior of traveling waves[END_REF], we numerically obtained propagating terraces and did not manage to catch these non-monotonic rapid traveling waves. Thus their existence remains a completely open problem. We point out here that this existence would be in sharp contrast with a nonexistence result by Alfaro and Coville for the nonlocal Fisher-KPP equation [START_REF] Alfaro | Rapid traveling waves in the nonlocal fisher equation connect two unstable states[END_REF]. Thanks to the λ-ω normal form, Sherratt manages also to find formulas for the amplitude γ and the speed c γ of the wave train. Nevertheless, it is quite tedious to reduce our three-component system to the appropriate two-component λ-ω system, as its phase "plane" is the unstable manifold of the steady state 1 (which is definitely not a Euclidean plane). For the sake of brevity, we choose to omit here this reduction and the precise predictions on the wave train.

The contrast with the very simple dynamics exhibited by Cantrell, Cosner and Yu for the two-component system is striking. This is of course reminiscent of the contrast between two-component and three-component competitive Lotka-Volterra systems: the two-component ones always have a simple monostable or bistable structure, devoid of periodic orbits, whereas some three-component ones have stable limit cycles (as established by Zeeman [START_REF] Zeeman | Hopf bifurcations in competitive three-dimensional Lotka-Volterra systems[END_REF] in her classification of the 33 stable nullcline equivalence classes). However let us emphasize once more that our system is not competitive near 1, so that qualitatively similar observations for three-species competitive systems (for instance, those of Petrovskii-Kawasaki-Takasu-Shigesada [START_REF] Petrovskii | Diffusive waves, dynamical stabilization and spatio-temporal chaos in a community of three competitive species[END_REF]) are actually unrelated to our results.

We emphasize that here all diffusion rates are equal (and normalized), whence there is no Turing instability with respect to the space variable. Obviously, if the phenotypes differ also in diffusion rate, then even more complicated dynamics are to be expected -and can be observed numerically. On this vast topic, we refer for instance to Smith-Sherratt [START_REF] Smith | The effects of unequal diffusion coefficients on periodic travelling waves in oscillatory reaction-diffusion systems[END_REF].

This collection of results confirms that the traveling waves constructed in [START_REF] Girardin | Non-cooperative Fisher-KPP systems: traveling waves and long-time behavior[END_REF] are definitely not the end of the story from the viewpoint of the asymptotic spreading for the Cauchy problem. 1.4. What about more general systems? Theorem 1.1,1-2 (and its various consequences) can be easily extended to general KPP systems with Lotka-Volterra competition, equal diffusion rates and any number of components provided the matrices L and C remain circulant matrices and a Hopf bifurcation does occur (thus some asymmetry is required). Theorem 1.1,3 needs a bit more care and appropriately chosen coefficients but should still hold true in a much more general framework.

Here, we choose to focus on a particular three-dimensional example, mainly because our point is to confirm the existence of oscillatory KPP systems. There are secondary reasons worth mentioning: first, the particular choice we make simplifies a lot the notations and calculations; second, Hopf bifurcations with a hyperbolic transverse component are at their core a three-dimensional phenomenon and taking into account more dimensions is just cumbersome. Recall that j + j = -1, z • z = 1 √ 3 z and z • z = 1 3 1. Recall also that the set of all n × n circulant matrices forms a commutative algebra and that the matrix

U = 1 √ 3   1 1 1 1 j j 1 j j   is a unitary matrix such that   a b c c a b b c a   = U   a + b + c 0 0 0 a + bj + cj 0 0 0 a + bj + cj   U T .
These basic properties bring forth a very convenient spectral decomposition for the problem.

Uniqueness of the constant positive steady state.

Proof. Let µ > 0 and v ≥ 0 be a solution of (

I + µM) v -Cv • v = 0.
Writing v = α1+βz+γz with α, β, γ ∈ C and identifying the real and imaginary parts, we straightforwardly verify that α ∈ [0, +∞) (by nonnegativity of v) and β = γ (by the fact that v is a real vector).

Then the equality (

I + µM) v = Cv • v reads α1 + (1 -3µ) βz + (1 -3µ) βz = α1 + 7j 10 βz + 7j 10 βz • α1 + βz + βz .
After a few algebraic manipulations, this is equivalent to the system

α 2 -α -7 30 |β| 2 = 0 7 √ 3 30 jβ 2 + α -(1 -3µ) + 7 10 -1+i √ 3 2 α β = 0.
On one hand, assuming by contradiction the existence of a solution such that |β| = 0 and taking the square of the modulus of the second line multiplied by 20, we deduce 196 3

|β| 2 = (13α -20 (1 -3µ)) 2 + 147α 2 = 316α 2 -520 (1 -3µ) α + 400 (1 -3µ) 2 , that is 49 3 |β| 2 = 79α 2 -130 (1 -3µ) α + 100 (1 -3µ) 2 .
On the other hand, from the first line, we deduce 49 3 |β| 2 = 70 α 2 -α .

Equalizing the two expressions of 49 3 |β| 2 , we obtain

9α 2 + 10 (7 -13 (1 -3µ)) α + 100 (1 -3µ) 2 = 0.
On one hand, the discriminant of this equation is 100

(7 -13 (1 -3µ)) 2 -36 (1 -3µ) 2 ,
which is itself nonnegative if and only if

7 -26 (1 -3µ) + 19 (1 -3µ) 2 ≥ 0,
that is if and only if 1 -3µ / ∈ 7 19 , 1 , that is if and only if µ / ∈ 0, 4 19 . Therefore, α being real, v cannot possibly exist if µ ∈ 0, 4 19 . Hence necessarily µ ≥ 4 19 . On the other hand, at α = 0, the polynomial α → 9α 2 + 10 (7 -13 (1 -3µ)) α + 100 (1 -3µ)

2 is nonnegative and with derivative 30 (13µ -2), which is positive since we now assume the necessary condition µ ≥ 4 19 > 2 13 . Therefore the polynomial has actually no zero in (0, +∞). Since v ≥ 0 and |β| = 0 imply together α > 0, we find a contradiction.

This exactly means that all solutions satisfy β = 0 and α 2 -α = 0, so that 0 and 1 are indeed the only solutions as soon as µ > 0. (This is of course consistent with a direct differentiation.)

Since µM-C is a circulant matrix, three complex eigenpairs are (-1, 1), (λ µ , z), λ µ , z , where

λ µ = -2µ - 1 10 + µ - 8 10 j + µ - 1 10 j = 3 7 60 -µ + i 7 √ 3 20
.

This proves indeed the local asymptotic stability when µ > 7 60 , the Hopf bifurcation at µ = 7 60 and, as the transverse component is hyperbolic, the uniqueness of the limit cycle C µ close to the bifurcation value.

Thereafter we will also need adjoint eigenvectors satisfying (µM -C)

T z = λz, that is (µM -C) T z = λz.
Using this time the fact that (µM -C) T is circulant, eigenpairs of it are (-1, 1), (λ µ , z), λ µ , z .

2.4.

The first Lyapunov coefficient. First, we recall a well-known statement (we refer, for instance, to Kuznetsov [22, Formula 5.39, p. 180]).

Theorem 2.1. Let N ∈ N such that N ≥ 2, I ⊂ R be an interval containing 0 and

f ∈ C 3 R N × R, R N such that f (0, η) = 0 for all η ∈ I.
Assume that, if η ∈ I ∩ (-∞, 0), 0 is a locally asymptotically stable steady state for the dynamical system u = f (u, η) , and that at η = 0 it undergoes a Hopf bifurcation (with a center subspace of dimension 2).

Let A ∈ M N (R), b : R N × R N → R N and c : R N × R N × R N → R N such that the Taylor expansion of u → f (u, 0) at 0 has the form f (u, 0) = Au + 1 2 b (u, u) + 1 6 c (u, u, u) + O |u| 4 .
Let q ∈ C N be an eigenvector of A associated with the purely imaginary eigenvalue λ ∈ iR + and p ∈ C N be an eigenvector of A T associated with -λ, normalized so that p T q = q T q = 1.

Then the Hopf bifurcation is supercritical, respectively subcritical, if the first Lyapunov coefficient

l 1 (0) = 1 2 |λ|
Re p T c (q, q, q) -2p T b q, A -1 B (q, q) + p T c q, (2iω

0 I N -A) -1 B (q, q)
is negative, respectively positive.

We are now in position to apply this theorem to our case.

Proof. Performing the changes of variable v = 1+w and µ = 7 60 -η and identifying the Taylor expansion of w → 7 60 M -C w -Cw • w at w = 0, which is actually an exact expansion, we find T with respect to the eigenvalue -λ, so that in the preceding statement we have p = q = z.

The most convenient way to identify one by one the terms involved in the expression of the first Lyapunov coefficient is to use again the properties of circulant matrices. Doing so, we find:

b (z, z) = -2 Re (z • Cz) = -2 Re 7j 10 z • z = 7 30 1, A -1 =   U   -1 0 0 0 λ 0 0 0 -λ   U T   -1 = U      (-1) -1 0 0 0 i 7 √ 3 20 -1 0 0 0 -i 7 √ 3 20 -1      U T = U    -1 0 0 0 -i 20 √ 3 21 0 0 0 i 20 √ 3 21    U T , A -1 b (z, z) = - 7 30 1, b z, A -1 b (z, z) = 7 30 (Cz + z) = 7 300 (10 + 7j) z, z T b z, A -1 b (z, z) = 7 300 (10 + 7j) , b (z, z) = -2z • Cz = - 14j 10 z • z = - 14 √ 3 30 jz, i 7 √ 3 10 I -A -1 =   U    1 + i 7 √ 3 10 0 0 0 -λ + i 7 √ 3 10 0 0 0 λ + i 7 √ 3 10    U T    -1 = U       1 + i 7 √ 3 10 -1 0 0 0 -i 7 √ 3 20 + i 7 √ 3 10 -1 0 0 0 i 7 √ 3 20 + i 7 √ 3 10 -1       U T = U    100-i70 √ 3 247 0 0 0 -i 20 √ 3 21 0 0 0 -i 20 √ 3 63    U T , i 7 √ 3 10 I -A -1 b (z, z) = - 14 √ 3 30 j i 7 √ 3 10 I -A -1 z = - 14 √ 3 30 j -i 20 √ 3 63 z = 4 9 ijz, b   z, i 7 √ 3 10 I -A -1 b (z, z)   = 4 9 ijb (z, z) = 4 9 ijb (z, z) = - 56 √ 3 270 iz = -i 28 √ 3 135 z, z T b   z, i 7 √ 3 10 I -A -1 b (z, z)   = -i 28 √ 3 135
Finally, the first Lyapunov coefficient of (KPP 0 µ ) is

l 1 (0) = 10 √ 3 21
Re -14 300 (10 + 7j) -i 28

√ 3 135 = - 10 √ 3 21 × 14 300 × 13 2
= -13 √ 3 90 and, consequently, the limit cycle C µ is indeed locally asymptotically stable close to the bifurcation value.

2.5. Continuation of the limit cycle when the Hopf bifurcation theorem does not apply anymore.

Proof. First, we show that, for any µ ≥ 0, the ω-limit set of (KPP 0 µ ) is, apart from 0, contained in the compact set

I = v ≥ 0 | 1 ≤ v 1 + v 2 + v 3 3 ≤ 10 3 .
Using again the decomposition v = α1 + βz + βz, this simply amounts to verifying that, for any α ∈ [0, 1],

(v + µMv -Cv • v) • 1 ≥ 0 and, for any α ≥ 10 3 , (v + µMv -Cv • v) • 1 ≤ 0.
Using again previous calculations, we end up with

(v + µMv -Cv • v) • 1 = α -α 2 + 7 30 |β| 2 ,
which is obviously nonnegative if α ∈ [0, 1]. Noticing the simple geometric fact that

T α = {v ≥ 0 | v 1 + v 2 + v 3 = 3α}
is an equilateral triangle of perimeter 9 √ 2α whose circumscribed circle is the boundary of the closed two-dimensional ball

B α = v ∈ B α1, √ 6α | v 1 + v 2 + v 3 = 3α , we deduce √ 6α = max v∈Tα |v -α1| = max v∈Bα |v -α1| = max v∈Bα |2 Re (βz)| = 2 max v∈Bα |β| Re e i arg(β) z = 2 √ 3 max v∈Bα |β| max θ∈[0,2π] cos (θ) 2 + cos θ + 2π 3 2 + cos θ + 4π 3 2 = √ 2 max v∈Bα |β| ≥ √ 2 max v∈Tα |β| whence α -α 2 + 7 30 |β| 2 ≤ α -3 10 α 2
, which is indeed negative if α > 10 3 (and, having in mind that 10e 1 = 10 3 1 + 2 Re 10 3 z is a steady state of the particular case µ = 0, this constant is optimal). Considering the dynamical system defined by (KPP 0 µ ) with initial conditions in the unstable manifold of 1, we can reduce it to a two-dimensional flow whose ω-limit set is also included in I. Applying the Poincaré-Bendixson theorem, we deduce for any value of µ ∈ (0, µ H ) the necessary existence of a positive limit cycle C µ in I.

Notice that although a limit cycle that is locally asymptotically stable for the flow embedded in the unstable manifold of 1 does exist, we do not, at this point, have any information on the stability of this limit cycle in the three-dimensional flow.

Next, using the relative compactness (in the topology induced by the Hausdorff distance) of any family (C µ ) µ∈(0,µH) , we can extract a limit point of it as µ → 0, say C. Fixing an appropriate family of initial conditions, we easily derive the existence of a solution of (KPP 0 µ ) with µ = 0 whose full trajectory is contained in C. The corresponding orbit is a fixed point, a limit cycle, a heteroclinic connection or a homoclinic connection. Since 1 does not bifurcate again at µ = 0, the case C = {1} is discarded. The well-known characterization of the ω-limit set of the three-component Lotka-Volterra competitive system corresponding to the case µ = 0 (see Zeeman [35, equivalence class n • 27, p. 22], Uno-Odani [START_REF] Uno | On a Lotka-Volterra model which can be projected to a sphere[END_REF], May-Leonard [START_REF] May | Nonlinear aspects of competition between three species[END_REF], Petrovskii-Kawasaki-Takasu-Shigesada [START_REF] Petrovskii | Diffusive waves, dynamical stabilization and spatio-temporal chaos in a community of three competitive species[END_REF], etc.) shows then that C is indeed a reunion of elements among {10e 1 }, {10e 2 }, {10e 3 }, H 1 , H 2 , H 3 (where we recall that H i is the heteroclinic orbit connecting 10e i and 10e i+1 ). In particular, the limiting system does not admit any periodic limit cycle.

In order to conclude, it only remains to prove that any limit cycle C µ encloses span [START_REF] Alfaro | Rapid traveling waves in the nonlocal fisher equation connect two unstable states[END_REF], so that in the end C = C 0 . To do so, we are going to show that the flow always crosses a plane containing the straight line span (1) in the same direction, that is we are going to show that, for any v = α1 + βz + βz with β = 0, The estimate |β| ≤ √ 3α is this time not precise enough; we truly need to relate the modulus of β and its argument.

(v + µMv -Cv • v) • e i π 2 βz + e i π
By periodicity and invariance by rotation around the axis span (1), it suffices to consider an interval of length 2π 3 for the parameter θ = arg (β). For instance, we take the interval 2π 3 , 4π 3 . In this interval, T α is characterized by the inequality v 1 ≥ 0, which reads α+2 |β| cos θ ≥ 0. Consequently, the studied sign is nonnegative provided

3 |cos θ| ≥ cos 3θ + π 6 for all θ ∈ 2π 3 , 4π 3
,
which is obviously true. Therefore the flow is rotating clockwise around span (1) if seen from 0 (consistently with Figure 3.1), and so is any periodic orbit. Thus any limit point C satisfies indeed C = C 0 , whence any full family (C µ ) µ∈(0,µH) converges as µ → 0 to C 0 .

It might be tempting to use the same ideas to localize more efficiently, and maybe even count, the limit cycles. However, the sign of

(v + µMv -Cv • v) • βz +

Numerical findings

In this section, µ = 13 120 ∈ (µ -, µ H ) is fixed. 3.1. The numerical scheme. All the forthcoming plots are obtained thanks to a simple finite difference scheme, explicit in time and with Neumann boundary conditions on the boundary of a very large spatial interval. It is well-known that such a spatial domain approximates correctly R, at least regarding spreading properties of reaction-diffusion systems and equations. Indeed, the forthcoming results are consistent with previously known theoretical results (such as, for instance, the fact that initially compactly supported solutions for (KPP µ ) invade 0 at speed 2 or the exponential decay of traveling wave solutions [START_REF] Girardin | Non-cooperative Fisher-KPP systems: traveling waves and long-time behavior[END_REF][START_REF] Girardin | Non-cooperative Fisher-KPP systems: Asymptotic behavior of traveling waves[END_REF]).

Source codes are run in Octave [START_REF] Eaton | GNU Octave version 5.1.0 manual: a high-level interactive language for numerical computations[END_REF].

The findings seem to be robust with respect to the numerical parameters.

3.2.

The limit cycle for (KPP 0 µ ). Although we do not know how to prove analytically the global attractivity or the uniqueness of the limit cycle C µ , numerically it seems indeed to be true, as illustrated by Figure 3.1.

The Cauchy problem with diffusion .

The following findings seem to be robust with respect to the initial condition u 0 , as soon as it is compactly supported, nonzero and not in span(1) (stable manifold of 1). For instance, we fix u 0 = (1.01, 1.01, 0.99) T in a small interval in the center of the domain and u 0 = 0 elsewhere.

Once the existence of a pulsating front connecting 1 to a wave train (t, x) → p γ (κ γ x -σ γ t) is observed (see Figure 3.2), we use the phase space to estimate the amplitude γ of the wave train. In order to do so, we plot in Figure 3.3 the trajectory of t → u(t, x) with x appropriatly chosen (say, away from the initial support of the solution but within the final support of the wave train) together with C (obtained by truncating any trajectory of (KPP 0 µ ), see Figure 3.1). This confirms that γ is smaller than, but close to, 1. As a side note, this also confirms that the selected wave train is a stable one (in the sense of Proposition 1.3).

To evaluate the speed of the pulsating front, the most convenient way is to plot an appropriate level set. Since the three components of u always spread together, it is sufficient to plot the level set of only one component, say u 1 . Of course, the value U of u 1 at this level set must satisfy 0 < U < 1, so for instance we fix U = 0.9. We obtain then Figure 3.4.

With Figure 3.4, we can verify that the invasion 1 → 0 occurs at speed 2 and we can evaluate graphically that the invasion p γ → 1 occurs at speed c 1 3 , which corresponds to the linear prediction of Sherratt [START_REF] Sherratt | Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave[END_REF]:

c lin = 2 Re λ 13 120 = 2 3 120 = 1 √ 10 0.3162.
Let us point out that Sherratt's prediction uses the parameter λ 0 of the λ-ω form of the system instead of the real part of the bifurcating eigenvalues and that the equality is perhaps not obvious. As explained earlier, the λ-ω reduction is not We also see on Figure 3.4 that the intrinsic speed c γ of the wave train is negative and, as expected, of large absolute value. (2) Same as before.

(3) Same as before. and then, from the convexity of V , it follows again that u has to be the unique steady state v , which is constant. (5) Same as before except we use as Lyapunov function the integral of V over a spatial period.

Recalling that the profile of a traveling wave connects, in some sense, 0 to an entire solution of the system satisfying the positivity condition above, we deduce directly various sufficient conditions for the convergence of the profile. In particular, all profiles of a monostable system with weak mutations [START_REF] Griette | Existence and qualitative properties of travelling waves for an epidemiological model with mutations[END_REF][START_REF] Morris | Individual variability in dispersal and invasion speed[END_REF] converge to the stable state indeed. This is a new step toward the resolution of a conjecture presented in an earlier work [17, Conjecture 1.1] (still, let us emphasize that the bistable case remains largely open).

Notice that the same conditions also guarantee that the nonnegative nonzero solutions of the Cauchy problem converge locally uniformly to a constant steady state. In particular, two-component systems with weak mutation rates, equal diffusion rates and a unique positive constant steady state satisfy this convergence property. This is of course in striking contrast with the three-component counterexample that is the main point of the present paper.

2 . 1 2. 1 .

 211 Proof of Theorem 1.Well-known facts on circulant matrices. The matrices I, M and C are all circulant matrices. Recall that the 3 × 3 circulant matrix   (a + b + c, 1), a + bj + cj, z and a + bj + cj,

2. 3 .

 3 The linearization at 1: eigenelements and Hopf bifurcation.Proof. The change of variable v = 1 + w leads to(I + µM) v -Cv • v = (µM -C) w -Cw • w.Hence the linearization of the reaction term at v = 1 is exactly µM -

7 60 M 7 √ 3 20.

 76073 -C w -Cw • w = Aw + 1 2 b (w, w) , where A = 7 60 M -C and b : (v, w) → -w • Cvv • Cw. Let λ = λ 7/60 = -λ 7/60 = i The vector z is an eigenvector of 7 60 M -C with respect to the eigenvalue λ and an eigenvector of 7 60 M -C

Figure 3 . 1 . 3 ( 1 ,

 3131 Figure 3.1. Seven trajectories of (KPP 0 µ ) with random initial conditions.

Figure 3 . 3 .Figure 3 . 4 .

 3334 Figure 3.3. In blue, the trajectory at x. In red, the limit cycle C.

( 4 )

 4 If the system is of mixed type and D = I, then by using the Lyapunov functionV = c 1 F 1 (u 1 ) + c 2 F 2 (u 2 ) of [9, Lemma 3.2], we find ∂ t V -∆V = -∇u T .D 2 V.∇u + ∇V • (∂ t u -∆u) = -∇u T .D 2 V.∇u + ∇V • (Lu -Cu • u)

  Assume that (KPP µ ) is set in the spatial domain R n with n ∈ N. If µ < µ H and C µ is locally asymptotically stable, then for any e ∈ S n-1 , (KPP µ ) admits a continuous one-parameter family of traveling plane wave train solutions of the form u : (t, x) → p γ (κ γ x • e -σ γ t)

	27, Chapter 1, Section 1.7]) and
	immediately obtain the following result.
	Proposition 1.2.

  Given a fixed angle arg β ∈ [0, 2π], the nullcline is a straight line in the plane of coordinates (α -1, |β|) whose slope is -14

			√ 39 cos 3 arg β + 2π 3 3 . Unfortunately the
	sign of this slope varies as arg β varies. Hence the best, and really unsatisfying,
	result we can deduce from this is that any limit cycle is in the region of the phase
	space where
	60 13	(µ H -µ) -(α -1) ∈ -	14 39 √	3	|β| ,	14 √ 39	3	|β| .
								βz
	is the same as the sign of 60 13 (µ H -µ) -(α -1) -	14 39 √	3	cos 3 arg (β) +	2π 3	|β| .

Acknowledgments

The author wishes to thank two anonymous referees whose comments helped to improve the clarity of the manuscript.

This work was supported by a public grant as part of the Investissement d'avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH. This work has been carried out in the framework of the NONLOCAL project (ANR-14-CE25-0013) funded by the French National Research Agency (ANR).

Proof. As explained by Maginu [START_REF] Maginu | Stability of periodic travelling wave solutions with large spatial periods in reaction-diffusion systems[END_REF], in order to establish the marginal stability in linearized criterion of wave trains p γ with an amplitude γ sufficiently close to 1, it suffices to prove the strong stability of the spatially homogeneous limit cycle p 1 . Roughly speaking, the strong stability in the sense of Maginu is the linear stability with respect to spatio-temporal perturbations of the form sin (ωx) u(t).

More precisely, the strong stability of the spatially homogeneous limit cycle C µ for (KPP µ ) is defined by Maginu [START_REF] Maginu | Stability of periodic travelling wave solutions with large spatial periods in reaction-diffusion systems[END_REF] as the negativity of all Floquet exponents of all systems

• v evalued at p 1 (which is the periodic profile corresponding to the limit cycle). Let U ω be the fundamental solution associated with (A.1), namely the solution of

It is easily verified that t → e ω 2 t U ω (t) is exactly U 0 . Therefore the Floquet exponents (η i (ω)) i∈{1,2,3} of (A.1) satisfy exactly

. The negativity of the family (η i (0)) i∈{1,2,3} leads to the conclusion.

Appendix B. A remark on the entire solutions of two-component KPP systems with Lotka-Volterra competition set in a Euclidean space

In this section we will use the terminology "eventually cooperative", "eventually competitive" and "mixed type". It refers to the trichotomy of Cantrell-Cosner-Yu [9, Figure 1, Proposition 2.5].

As a preliminary, we point out a result that was just hinted in Cantrell-Cosner-Yu [START_REF] Stephen Cantrell | Dynamics of populations with individual variation in dispersal on bounded domains[END_REF]: in the eventually competitive, bistable case, we can use classical arguments (unstable manifold theorem, Bendixson-Dulac theorem, Poincaré-Bendixson theorem) to show that, exactly as in the Lotka-Volterra case, there exists a partition (B 1 , S, B 2 ) of [0, +∞) 2 such that each B i is the basin of attraction of a stable steady state whereas the separatrix S is the basin of attraction of the nonzero unstable steady state. The separatrix is smooth, contains 0 and is, in the competitive rectangle, the graph of a nondecreasing function.

Proposition B.1. Let n ∈ N, D be a diagonal 2 × 2 matrix with positive diagonal entries, L be a 2 × 2 essentially nonnegative and irreducible matrix satisfying λ PF (L) > 0, C be a 2 × 2 positive matrix and u be an entire solution of

Then u is a constant steady state provided one of the following conditions holds true:

(1) the system is eventually cooperative;

(2) the system is eventually competitive and monostable;

(3) the system is eventually competitive, bistable and there exists t ∈ R such that the image of x → u(t, x) does not intersect the separatrix; (4) the system is of mixed type and D = I;

(5) the system is of mixed type and u is spatially periodic; Proof.

(1) If the system is eventually cooperative, then necessarily u is valued in the cooperative rectangle and, by comparison with a solution that does not depend on x and the fact that the unique steady state v is globally attractive for the diffusionless system, u is exactly v , which is constant. Assume that v is unstable, in the sense that at least one eigenvalue of the linearized operator L v = L -diag(v)C -diag(Cv) has a nonnegative real part.

Then L v is not essentially nonnegative.

Proof. By definition of v, L v v = -Cv • v. This vector is obviously negative. Assuming by contradiction that L v is essentially nonnegative, we deduce by standard properties of essentially nonnegative matrices (e.g., [START_REF] Horn | Topics in matrix analysis[END_REF]) that the Perron-Frobenius eigenvalue of L v , whose real part is maximal among the eigenvalues, is negative. The instability of v is contradicted.

With more general competition terms c(v) [START_REF] Girardin | Non-cooperative Fisher-KPP systems: Asymptotic behavior of traveling waves[END_REF][START_REF] Girardin | Non-cooperative Fisher-KPP systems: traveling waves and long-time behavior[END_REF], the same proof will work provided Dc(v).v is nonnegative. This is a fairly general assumption, reminiscent of the known condition for the existence of traveling waves [START_REF] Girardin | Non-cooperative Fisher-KPP systems: traveling waves and long-time behavior[END_REF]Theorem 1.5].

We also point out that the same simple observation (L v v is negative) yields other interesting properties, for instance the existence of an eigenvalue with negative real part in the case where L v is symmetric.