Sylvain Meignen

Thomas Oberlin

On the Study of Interference and their Localization in the Time-Frequency Plane

Keywords: Time-frequency analysis, short-time Fourier transform, spectrogram zeros, Delaunay triangulation

In this paper, our goal is first to investigate the conditions for the separation of the modes making up multicomponent signals based on the analysis of the short-time Fourier transform. More precisely, we put forward necessary and sufficient conditions for the existence of ridges associated with the modes in the time-frequency plane. The focus is put on signals either made of purely harmonic modes or parallel linear chirps, for which we show that when the modes are strongly interfering in the time-frequency plane, the ridges no longer exist and are replaced by some structures called time-frequency bubbles. Based on a careful study of interference patterns, we show that the zeros of the spectrogram involved in these share very specific features, on which we found a new algorithm to determine interference locations based on the Delaunay triangulation of spectrogram zeros.

I. INTRODUCTION

The analysis of multicomponent signals (MCSs) using timefrequency representations (TFRs) has been the subject of intense research in the last few decades, since the modes making up a MCS are associated with specific regions around TF curves, called ridges [START_REF] Delprat | Asymptotic wavelet and gabor analysis: Extraction of instantaneous frequencies[END_REF], [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF]. The locations of these ridges lead to estimates the instantaneous frequencies (IFs) of the modes [START_REF] Stankovic | Performance of spectrogram as IF estimator[END_REF], [START_REF] Carmona | Multiridge detection and time-frequency reconstruction[END_REF]. Among linear TFRs the most popular ones are the short-time Fourier Transform (STFT) and the continuous wavelet transform (CWT) while the most commonly used quadratic TFR is the spectrogram.

The quality of IF estimation from ridge extraction is tightly related to how well the modes are separated in the TF plane, namely to the choice of the analysis window in the case of spectrogram. Indeed, an inappropriate window choice results in mode interference. In that framework, optimal window length can be found by minimizing the Rényi entropy which is proved to be a good measure of the level of interference in the TF plane [START_REF] Baraniuk | Measuring time-frequency information content using the Rényi entropies[END_REF], [START_REF] Meignen | On the use of Rényi entropy for optimal window size computation in the short-time Fourier transform[END_REF]. However, even when one considers the TFR minimizing the Rényi entropy, there is no guarantee that the modes are not interfering in the TF plane. To improve mode separation for IF estimation, the socalled adaptive short-time Fourier transform (ASTFT) was introduced [START_REF] Kwok | Improved instantaneous frequency estimation using an adaptive short-time fourier transform[END_REF], [START_REF] Zhong | Time-frequency representation based on an adaptive short-time fourier transform[END_REF] and recently adapted to the context of synchrosqueezing transform [START_REF] Li | Adaptive synchrosqueezing transform with a time-varying parameter for non-stationary signal separation[END_REF], [START_REF] Li | Adaptive short-time Fourier transform and synchrosqueezing transform for non-stationary signal separation[END_REF]. Alternative techniques aiming at improving IF estimation from reassigned transform S. Meignen is with the Jean Kuntzmann Laboratory, University Grenoble Alpes and CNRS 5225, Grenoble 38401, France (email:sylvain.meignen@univ-grenoble-alpes.fr), while T. Oberlin is with ISAE-SUPAERO, Université de Toulouse, 31055 Toulouse CEDEX 4, France. This work is supported by the ANR-ASCETE project with grant number ANR-19-CE48-0001-01. [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF] were also developed in [START_REF] Bruni | On the time-frequency reassignment of interfering modes in multicomponent FM signals[END_REF], [START_REF]A fast and robust spectrogram reassignment method[END_REF], the basic principle of these approaches being to disentangle the components at interference locations by modifying the reassignment process. Note that reassignment approaches are also very common to deal with overlapping transients, a situation frequently encountered in echolocation calls [START_REF] Reinhold | Objective detection and time-frequency localization of components within transient signals[END_REF], [START_REF] Brynolfsson | The matched reassignment applied to echolocation data[END_REF]. In this regard, multi-taper approaches [START_REF] Xiao | Multitaper time-frequency reassignment for nonstationary spectrum estimation and chirp enhancement[END_REF] result in significant improvement in terms of transient separation on that type of signals [START_REF] Sandsten | A multitaper reassigned spectrogram for increased time-frequency localization precision[END_REF].

To adapt locally the window length to reduce interference, one often uses the argument that the sparsest representations should lead to the best IF estimation. Nevertheless, such an adaption is not necessary at each time instant to obtain good mode separation. To have an insight into where to modify the window length, an analysis of the TF interference pattern can be carried out. It was remarked in [START_REF] Delprat | Global frequency modulation laws extraction from the gabor transform of a signal: A first study of the interacting components case[END_REF] that the TFR associated with the Gabor transform of the sum of two pure harmonic modes could exhibit so-called time-frequency bubbles (TFBs), consisting of circles of local maxima at the TF location of the two modes. Such structures were also observed in more complex situations [START_REF] Delprat | Global frequency modulation laws extraction from the gabor transform of a signal: A first study of the interacting components case[END_REF]. In practice, TFBs arise when, at some time instant, two modes generate only a single maximum of the spectrogram. When TFBs are present, the estimation of the IF of the modes based on ridge extraction is no longer possible. Though in [START_REF] Delprat | Global frequency modulation laws extraction from the gabor transform of a signal: A first study of the interacting components case[END_REF] some elements were given to explain when such TFBs occured, a clear mathematical analysis was lacking. In the present paper, our first goal is thus to bridge that gap, putting forward a necessary and sufficient condition for the existence of TFBs when the signal is made of two pure harmonic modes and when the TFR is the spectrogram. Then, we show how to extend these conditions to the case of the sum of two parallel linear chirps.

More generally, the spectrogram of close modes can be viewed as a continuum from separate ridges to TFBs, the latter corresponding to the highest level of interference. For IF estimation based on ridge detection, it is crucial to determine where interference occur in the TF plane, because even a low interference level can reduce drastically the quality of IF estimation from TF ridges. For that purpose, we analyze interference patterns, and show that these involve local maxima, zeros and saddle points of the spectrogram. Based on this study, we propose a novel algorithm to localize interference in the TF plane based on Delaunay triangulation of the zeros of the spectrogram involved in interference patterns.

The motivation for using the zeros of the spectrogram can be found in the pioneering work by Flandrin [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF]. In that paper, the zeros of the spectrogram of Gaussian white noise were studied, and further investigated in [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF]. The conclusions of these papers were that, on the one hand, the spectogram zeros tend not to clutter and, on the other hand, they cannot be too far one from another. Based on this last remark, a new technique using Delaunay triangulation of spectrogram zeros was proposed to separate the noise from the signal in the spectrogram [START_REF] Meignen | Adaptive multimode signal reconstruction from time-frequency representations[END_REF]. The rationale supporting the algorithm proposed in that paper was basically that the signal corresponded to zero-free regions, and thus the surrounding zeros generated triangles with long edges in the Delaunay triangulation of the spectrogram zeros. Thus, the criterion chosen to separate the triangles related to noise from those corresponding to the signal was based of the analysis of the edge length in the triangulation. We here propose to develop a variant to this approach to localize interference in the TF plane. More precisely, we study in details the nature of the triangles surrounding the zeros involved in the Delaunay triangulation and then found our algorithm to localize interference on that analysis.

The paper is organized as follows. In the following section, we introduce the notation that we will use throughout the paper. Then, we state the main theoretical contributions regarding the existence of TFBs, studying successively the case of a signal made either of two pure harmonic modes or of parallel linear chirps. Finally, we propose a novel technique to localize mode interference in the TF plane based on the Delaunay triangulation of spectrogram zeros.

II. NOTATION

Considering a signal f ∈ L 1 (R)∩L 2 (R) and a real window h ∈ L ∞ (R) ∩ L 2 (R), the (modified) short-time Fourier transform (STFT) is defined as:

V h f (t, η) = R f (τ )h(τ -t)e -2iπ(τ -t)η dτ, (1)
the spectrogram being the square modulus of V h f . In the sequel, we are going to study the interference between the modes of a multicomponent signal (MCS), defined as the superposition of P AM-FM components:

f (t) = P p=1 f p (t) with f p (t) = A p (t)e 2iπφp(t) , (2)
in which the instantaneous amplitudes (IAs) A p (t) are positive, as well as the instantaneous frequencies (IFs) φ p (t). We further assume that the IFs are sorted in decreasing order, namely φ p+1 (t) > φ p (t) for all t.

III. FREQUENCY DOMAIN INTERFERENCE

In this section, we first investigate the existence of TFBs when the signal is the sum of two pure harmonic modes. We then explain how the obtained results can be extended to the case of two parallel linear chirps. More precisely, we show that the existence of these TFBs is related to that of extrema of the spectrogram. For our study, we will use the Gaussian window h(t) = e -π t 2 σ 2 to define the STFT, because with that choice analytical results can be derived.

A. Existence of TFBs in the Case of the Sum of the Two Pure Harmonic Signals

Let f (t) = f 1 (t)+f 2 (t) with f 1 (t) = Ae i2πξ1t and f 2 (t) = e i2πξ2t , where ξ 1 < ξ 2 . For such a signal, one has V h f1 (t, η) = ĥ(η -ξ 1 )Ae i2πξ1t = σAe i2πξ1t e -π(η-ξ1) 2 σ 2 and V h f2 (t, η) = σe i2πξ2t e -π(η-ξ2) 2 σ 2 and thus:

|V h f (t, η)| 2 = σ 2 (A 2 e -2πσ 2 (η-ξ1) 2 + e -2πσ 2 (η-ξ2) 2 +2Ae -πσ 2 [(η-ξ1) 2 +(η-ξ2) 2 ] cos(2π(ξ 2 -ξ 1 )t)) (3)
which attains, considering η fixed, its maximum at t k = k ξ2-ξ1 , k ∈ Z for which:

|V h f (t k , η)| 2 = σ 2 (Ae -πσ 2 (η-ξ1) 2 + e -πσ 2 (η-ξ2) 2 ) 2 . (4) Studying |V h f (t k , .)|
| log(A)| < 2α α 2 -1 -log α + √ α 2 -1 α - √ α 2 -1 = 2α α 2 -1 -2 arccosh(α). (5)
In any other case, |V h f (t k , .)| 2 has a unique extremum which is a maximum.

The proof is given in Appendix A. Note that one can easily show that the extrema are all in the interval ]ξ 1 , ξ 2 [.

If A = 1, remarking that log α+ √ α 2 -1 α- √ α 2 -1 = 2
arccosh(α) the condition to have three extrema reads:

2α α 2 -1 -log α + √ α 2 -1 α - √ α 2 -1 > 0 ⇔ -arccosh(α) + α α 2 -1 > 0. (6)
As the functionarccosh(x)+x √ x 2 -1 is strictly increasing on [1, +∞[ and is zero at zero, the condition boils down to α > 1.

Now, η being fixed, let us consider the time instants where |V h f (t, η)| 2 attains its minima which correspond to time t = tk = 1/2+k ξ2-ξ1 , k ∈ Z, at which: and one otherwise.

|V h f ( tk , η)| 2 = σ 2 (e -πσ 2 (η-ξ1) 2 -Ae -πσ 2 (η-ξ2) 2 ) 2 . ( 7
The proof is available in Appendix C. This can be generalized to any time t and any amplitude A through the following:

Proposition III.
| log(A)| < -arccosh(X 2 ) + 2α 2 X 2 2 -1 X 2 + γ , (9)
with X 2 = γ(α 2 -1)+α γ 2 (α 2 -2) + 2. In any other case, |V h f (t, .)| 2 has a unique extremum.

The proof is available in Appendix D. Note that when γ = 1, we get X 2 = 2α 2 -1, and then

-arccosh(X 2 ) + 2α 2 X 2 2 -1 X 2 + 1 = -arccosh(X 2 ) + X 2 2 -1 = -arccosh(2α 2 -1) + 2α α 2 -1 (10) = -2 arccosh(α) + 2α α 2 -1
which is exactly the condition of Proposition III.1. Similarly, when A = 1 the right hand side of ( 9) is positive and the condition to have three extrema boils down to α > 1+γ 2 , which is the condition enounced in Proposition III.3. Now, to finish this study we state a necessary and sufficient condition to have three extrema:

Proposition III.5. |V h f (t, .)| 2 whatever t has three extrema if and only if α := π 2 σ(ξ 2 -ξ 1 ) > 1 and | log(A)| < -2 arccosh(α) + 2α α 2 -1. (11)
In any other case there exists some time t where |V h f (t, .)| 2 has a unique extremum.

The proof is given in Appendix E. This last proposition means that if the spectrogram has three extrema at t k then it has three extrema at each time t. In that context, the existence of TFBs is related to the conditions of Proposition III.1 being not fulfilled. In such a case, each mode cannot be associated with a ridge. To illustrate this, we consider the spectrogram of the sum of two pure harmonic modes with A = 1 such that the conditions of Proposition III.1 are fulfilled in Fig. 1 (a) and not in Fig. 1 (b), in which TFBs appear. Finally, in Fig. 1 (c), we consider the case A = 1 when the conditions of Proposition III.1 are fulfilled. We notice that the locations of the extrema of the spectrogram at time instants tk , and in particular the locations of the zeros, are not much changed. But, this is not the case of the extrema locations at time t k .

B. Existence of TFBs in the Case of two Parallel Linear Chirps

Our goal is now to extend the previous study to parallel linear chirps. We thus consider f (t) = f 1 (t) + f 2 (t) with f 1 (t) = Ae 2iπφ1(t) and f 2 (t) = e 2iπφ2(t) , with φ 1 (t) = φ 2 (t) = C constant. It can be easily shown that in such a case one has:

|V h f (t, η)| 2 = r -1 A 2 e - 2πσ 2 (η-φ 1 (t)) 2 1+C 2 σ 4 + e - 2πσ 2 (η-φ 2 (t)) 2 1+C 2 σ 4 +2Ae -π σ 2 [ (η-φ 1 (t)) 2 +(η-φ 2 (t)) 2 ] 1+C 2 σ 4 cos(2πφ(t, η)) , (12)
with r =

1 σ 4 + C 2 and φ(t, η) = φ 2 (t) -φ 1 (t) + Cσ 4 1 + C 2 σ 4 (φ 2 (t) -φ 1 (t)) η - φ 2 (t) + φ 1 (t) 2 . ( 13
)
Now, given a real number γ in [-1, 1], let us consider the set of TF points (t, η) such that cos(2πφ(t, η)) = γ. These points are such that φ(t, η) = arccos(γ)

2π + k = β k . As φ 1 (t) = a 1 t + C 2 t 2 and φ 2 (t) = a 2 t + C 2 t 2
, we have that: From this, we deduce that

(a 2 -a 1 ) t + Cσ 4 1 + C 2 σ 4 η - a 1 + a 2 + 2Ct 2 = β k ⇔ t = -Cσ 4 η + β k (1 + C 2 σ 4 ) + (a 2 2 -a 2 1 )Cσ 4 2 a 2 -a 1 ⇔ t = Dη + B k . (14)
|V h f (Dη + B k , η)| 2 = r -1 A 2 e - 2πσ 2 (η-a 1 -C(Dη+B k )) 2 1+C 2 σ 4 + e - 2πσ 2 (η-a 2 -C(Dη+B k ) 2 1+C 2 σ 4 +2Ae -π σ 2 [ (η-a 1 -C(Dη+B k )) 2 +(η-a 2 -C(Dη+B k )) 2 ] 1+C 2 σ 4 γ = r -1 A 2 e -2π σ2 (η-a 1 +CB k 1-CD ) 2 + e -2π σ2 (η-a 2 +CB k 1-CD ) 2 + 2Ae -π σ2 (η- a 1 +CB k 1-CD ) 2 +(η- a 2 +CB k 1-CD ) 2 γ , (15)
with σ = σ(1-CD) √ 1+C 2 σ 4 = σ √ 1 + C 2 σ 4 .
As the expression (15) is the same as (3) replacing σ by σ and ξ 1 and ξ 2 respectively by a1+CB k 1-CD and by a2+CB k 1-CD . This leads to the following properties Proposition III.6. Assume (t, η) is evolving on the line t = Dη + B k such that cos(2πφ(t, η) = 1, then |V h f (t, η)| 2 has three extrema (two maxima and one minimum) along the line t = Dη + B k if and only if:

π 2 σ(a 2 -a 1 ) > 1,
and a unique maximum in that interval otherwise.

The proof is the same as that of Proposition III.1. If we now consider the case γ = -1, we have the following property:

Proposition III.7. Assume (t, η) is on the line t = Dη + B k , such that cos(2πφ(t, η) = -1, |V h f (t, η)| 2 always has, along the line t = Dη + B k , a unique minimum for η ∈ [ a1+CB k 1-CD , a2+CB k 1-CD ], a unique maximum for η ∈ ] a2+CB k 1-CD , +∞[, and a unique maximum in the interval ] - ∞, a1+CB k 1-CD [
The proof is the same as that of Proposition III.2. Going further we can prove that for any γ ∈ [-1, 1]:

Proposition III.8. If A = 1, |V h f (t, η)| 2 has three extrema along the line t = Dη+B k if and only if α C = π 2 σ(a 2 -a 1 )
satisfies:

α C > 1 + γ 2 ,
and one otherwise.

The proof of Proposition III.8 follows the same line as that of Proposition III.3. The arguments can be generalized to any A by adapting Proposition III.4 as follows:

Proposition III.9. Assuming α C is defined as in Proposition III.8, and that γ = cos(2πφ(t)), |V h f (t, η)| 2 has three extrema along the line t = Dη + B k if and only if α C > 1+γ 2 and

| log(A)| < -arccosh(X 2 ) + 2α 2 C X 2 2 -1 X 2 + γ , ( 16
)
with X 2 = γ(α 2 C -1) + α C γ 2 (α 2 C -2) + 2. In any other case, |V h f (t, η)| 2 has a unique extremum along the line t = Dη + B k .
We finally state a necessary and sufficient condition to have three extrema along the direction D:

Proposition III.10. |V h f (t, η)| 2 has three extrema along the line t = Dη + B k if and only if α C := π 2 σ(a 2 -a 1 ) > 1 and | log(A)| < -2 arccosh(α C ) + 2α C α 2 C -1. (17)
In any other case, where |V h f (t, η)| 2 has a unique extremum along the line t = Dη + B k .

An illustration of this is given in Fig. 2 (a) where we display the spectrogram of two parallel linear chirps for which there always exist 3 extrema in the direction D, making ridge extraction tractable, then, in Fig. 2 (b), we display the case where there exists only one extremum in the direction D for some value of γ, and thus TFBs are present. Finally, in Fig. 2 (c) we display a case where there are still three extrema in the direction D but when A is no longer equal to 1. Note that, as the computation of saddle points involves the determination of the direction D, we do not display extrema in that case, but the conclusion regarding the stability of the zeros observed in the case of two pure harmonic modes remains valid here.

IV. DETERMINATION OF THE INTERFERENCE REGIONS USING SPECTROGRAM ZEROS A. Motivation for Using Spectrogram Zeros

The previous study tells us that when two modes are interfering, there always exist zeros at some time instants located between the IF of the modes, even when TFBs are present. On the contrary, as shown in the previous section, the number of local maxima of the spectrogram involved in interference patterns may vary both in terms of their number and locations. Furthermore, we numerically notice that the saddle points are very unstable in noisy situations. For these reasons, to localize interference in the TF plane we are not going to directly use the saddle points or the maxima of the spectrogram, but rather the zeros generated by mode interference. Another important reason for using spectrogram zeros is that when the window is the unit Gaussian window, as remarked in [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF], STFT admits a Hadamard-Weirstrass decomposition, that can be fully determined by the locations of the spectrogram zeros [START_REF] Boas | Entire functions[END_REF]- [START_REF] Korsch | On the zeros of the husimi distribution[END_REF].

To build our algorithm to localize mode interference in the TF plane, our source of inspiration is the paper by Flandrin [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF] in which the properties of the Delaunay triangulation of the zeros of the spectrogram of Gaussian white noises were investigated. Indeed, it was shown in that paper that the edge length in such a triangulation is bounded above by a maximum value L e (depending on the normalization of the analysis window), and also that the zeros cannot clutter. Note that a deeper mathematical analysis of the properties of these zeros from the perspective of point processes is available in [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF].

B. Localizing Interference Using Delaunay Triangulation of Spectrogram Zeros

In this section, we are going to explain how to use the results of Section III to build an algorithm to localize interference in the TF plane. First, as remarked in [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF] the presence of signal in the TF plane gives rise to triangles in the Delaunay triangulation of spectrogram zeros (DTSZ) having at least one edge with length exceeding L e . Then, another interesting remark is that when two modes are interfering the zeros involved in the interference pattern belong to more triangles of that kind than other zeros of the spectrogram. To quantify this, we first propose to keep only the triangles of the DTSZ that contain at least one edge larger than a prescribed threshold. To do so, we sort out the edge lengths in ascending order and then consider the triangles that contain an edge among the P e % longest, the choice for that parameter being discussed later. At the end of this procedure, we obtain a set of triangles we denote by T (P e ) from now on. Then we compute the number of such triangles the zeros of the spectrogram belong to by associating with each vertex v (we use bold notation because v has two components) of the triangulation T (P e ) a connectivity index as follows:

I(v, P e ) = #{T ∈ T (P e ), v ∈ T }, (18)
where #X denotes the cardinal of the set X. In this regard, we also define the set of corresponding triangles as follows:

T v (P e ) = {T ∈ T (P e ), s.t. v ∈ T } . (19)
To get an insight into what an appropriate value for I(v, P e ) would be, we consider the two different signals corresponding to the spectrograms depicted in Fig. 3 (a) and (d). In Fig. 3 (a), we superimpose on top of the spectrogram the DTSZ (green edges), the white circles representing the vertices v such that I(v, 30) ≥ 1 (input SNR equals 20 dB). Then, in Fig. 3 (b) and (c), the white circles correspond to I(v, 30) ≥ 4 and I(v, 30) ≥ 5 respectively. On the second row of Fig. 3, we perform the same computation as that on the first row except we consider the spectrogram of Fig. 3 (d). These first simulations show that the connectivity index is a good parameter to discriminate the zeros associated with interference pattern from the others. The value for P e has to be chosen all the smaller that more zeros are associated with the signal and finally, by choosing a connectivity index of 5 the zeros associated with interference are kept while many of the other zeros are left apart. In what follows, we denote by "first step" of the algorithm the just described process to select the zeros of the spectrogram.

However, such a simple criterion is not enough to fully discriminate the zeros associated with interference. To go further, a first naive approach would consist of recalling that, from Section III, the distance between the zeroes located between two interfering pure harmonic modes is monitored by 1 ξ2-ξ1 so when the modes are close, not only the edges of the triangles crossing the modes should be long, but also the one in the direction of interference is not arbitrary and should larger than L e when the modes are close. However, as the instantaneous frequencies of the interfering modes are unknown, it is very difficult to build a criterion on that remark and alternatively, we propose to use some information on the energy on the triangles involved in the computation of the connectivity index.

More precisely, consider Φ a zero mean complex Gaussian white noise, then V g Φ is also Gaussian with zero mean and satisfies:

Var ( {V g Φ }) = Var ( {V g Φ }) = σ 2 g 2 2 , (20)
where {X} (resp. {X}) denotes the real (resp. imaginary) part of complex number X, Var the variance, and g 2 is the l 2 norm of the window g. From this, we deduce that

|V g Φ | 2 σ 2 g 2 2 is
χ 2 distributed with 2 degrees of freedom. Then, using table of χ 2 law the probability of false alarm is less than 1 % if |V g Φ | is larger than 3σ g 2 . So, to separate the noise from the signal in the TF plane, one first computes an estimate of γ = σ g 2 using the robust estimator proposed in [START_REF] Donoho | Ideal spatial adaptation via wavelet shrinkage[END_REF], namely:

γ := median | V g f | 0.6745 , (21)
in which f = f + Φ denotes the noisy signal.

Based on this analysis of the STFT of the noise, for each vertex selected at the end of the first step, we compute the proportion of triangles involved in the computation of the connectivity index that are above the noise level by computing the following index, for each v such that I(v, P e ) ≥ 5:

I signal (v, P e , β) = # T ∈ T v (P e ), ∃(t, η) ∈ T s.t. |V g f (t, η)| ≥ βγ #T v (P e ) . ( 22
)
An illustration of the zeros selected using that procedure is given in Fig. 4, when the noise level varies. On the first row, we display the zeros selected by considering I(v, 30) ≥ 5 (first step), when the input SNR is equal to 20, 10, and 5 dB, respectively in Fig. 4 (a), (b) and (c). We notice that though some zeros are not associated with interference pattern, the selected set of zeros contain those we want to discriminate. On the second row of Fig. 4, we show that the criterion ( 22) is very efficient to discriminate the zeros of the interference region from the others. Indeed, in that case, we plot the mask associated with |V g f (t, η)| ≥ βγ for an appropriate value for β (yellow zones on Fig. 4. (d), (e) and (f)), the value of which will be discussed later. Then, we select the zeros such that I signal (v, P e , β) = 1. These correspond to red circles on Fig. 4 (d),(e) and (f), which exactly correspond to the zeros associated with interference. Note that both the zeros associated with TFBs and those associated with more classical interference patterns are detected.

V. PARAMETER DETERMINATION

The algorithm proposed for the determination of the zeros associated with interference depends on three parameters P e , then the value of I(v, P e ), and finally β. Our strategy in that matter is to try to discriminate the zeros associated with a single mode from those associated with several modes. In this regard, we remark that, independently of the choice for P e , the zeros associated with an isolated ridge seldom belong to 5 or more triangles crossing the ridge (such a situation only happens when the ridge has some curvature), so to consider I(v, P e ) larger or equal to 5 seems a good strategy provided P e is chosen sufficiently large so as to select only the triangles crossing the ridges when computing the connectivity index. When P e is small many triangles are associated with noise and then I(v, P e ) is very often larger or equal to 5, but when P e increases the number of vertices selected using I(v, P e ) ≥ 5 should stabilize and be related to the signal. In Fig. 5 we plot the proportion of zeros kept with respect to P e assuming I(v, P e ) ≥ 5, for the spectrograms of Fig. 4 (a), (b) and (c), corresponding to input SNRs of 20, 10 and 5 dB, respectively. We observe that there exist two different regimes, the first one corresponds to a fast decay of the number of zeros kept and is related to the elimination of the zeros corresponding to noise, while the smaller decay observed for larger P e corresponds to the progressive elimination of the zeros associated with the signal. To better illustrate this aspect, we also plot, in Fig. 5, the linear approximation of the data for the different regimes. In this regard, we observe that P e = 30 is a good choice and that the behavior is similar regardless of the noise level. The last parameter we need to fix is β which is used in [START_REF] Boas | Entire functions[END_REF] to fully characterize the zeros associated with interference pattern. To do so, we still consider the spectrograms of Fig. 4 (a), (b) and (c), and that I(v, P e ) ≥ 5, and set P e = 30. Of course, this last parameter should be fixed differently depending on the signal. Our strategy to determine the right choice for β is then based on the fact that when β increases the zeros related to noise should be progressively eliminated by considering the vertices v such that I signal (v, 30, β) = 1. Indeed, when β = 3, the probability that |V g f (t, η)| ≥ βγ at location where only the noise is present is below 1 %. In Fig. 6 we compute the number of vertices v such that I signal (v, 30, β) = 1 with respect to β and again for the spectrogram of Fig. 4 (a), (b) and (c). We notice that, for β ≥ 3, the number of zeros reaches a plateau, this corresponds to all the zeros corresponding to noise have been eliminated, only remain those of interest. We then notice that these zeros are robust to an increase in the value of β, and all the more so that the noise level is low. In that context, to select the right zeros one should consider a value for β associated with such the plateau regime mentioned above, and thus the possible choices depend on the noise level.

VI. CONCLUSION

In this paper, we first put forward necessary and sufficient conditions for time-frequency bubbles to exist when two parallel pure harmonic modes or linear chirps are interfering. We notice that these conditions are very similar except that the modulation has to be taken into account when considering linear chirps. In the second part of the paper, after noticing that the Delaunay triangulation of spectrogram zeros have very interesting properties, we proposed an algorithm to localize interference in the TF plane based on the analysis of this triangulation. Future work should involve the improvement of the proposed criteria to localize interference regions in heavier noise situations, and it would also certainly be interesting to investigate whether the mathematical characterization of the appearance of TF bubbles could be extended to non parallel modes.

VII. APPENDIX

A. Proof of Proposition III.1

Let us consider the function, (Ae -πσ 2 (η-ξ1) 2 + e -πσ 2 (η-ξ2) 2 ) 2 whose extrema are the same as those of l(η) = Ae -πσ 2 (η-ξ1) 2 + e -πσ 2 (η-ξ2) 2 .

Defining ξ = ξ2-ξ1 2 and making the change of variables η = ν + ξ1+ξ2 2 , we get:

l(ν + ξ 1 + ξ 2 2 ) = Ae -πσ 2 (ν+ξ) 2 + e -πσ 2 (ν-ξ) 2 ,
which we denote by q(ν). Differentiating q(ν), we get after rearranging the terms:

q (ν) = -2πσ 2 e -πσ 2 (ν 2 +ξ 2 ) A(ν + ξ)e -2πσ 2 ξν + (ν -ξ)e 2πσ 2 ξν = 2πσ 2 e -πσ 2 ν 2 e -πσ 2 ξ 2 ξ(Ae -2πσ 2 ξν + e 2πσ 2 ξν ) - ν ξ - Ae -2πσ 2 ξν -e 2πσ 2 ξν Ae -2πσ 2 ξν + e 2πσ 2 ξν = 2πσ 2 e -πσ 2 ν 2 e -πσ 2 ξ 2 ξ(Ae -2πσ 2 ξν + e 2πσ 2 ξν ) - ν ξ + tanh(2πσ 2 ξν - log(A) 2 ) . (23)
So q has the same sign as g(ν) = -ν ξ + tanh(2πσ 2 ξν -

log(A) 2
) whose derivative is:

g (ν) = 1 ξ -1 + 2πσ 2 ξ 2 (1 -tanh 2 (2πσ 2 ξν - log(A) 2 )
which is always negative if 2πσ 2 ξ 2 ≤ 1. In such a case, q is decreasing and then and as lim ν→-∞ q (ν) = +∞ and lim ν→+∞ q (ν) = -∞, it annihilates once and thus l has a unique extremum which is a maximum. Now when 2πσ 2 ξ 2 > 1, g has two zeros, and we have the following table of variations:

ν -∞ ν 1 ν 2 +∞ g (ν) -0 + 0 - +∞ g(ν 2 ) g(ν) g(ν 1 ) -∞
If g(ν 1 ) ≥ 0, then g(ν) and hence q (ν) are positive on ] -∞, ν 2 [ and g(ν) and hence q (ν) annihilate and change signs at some ν 2 in ]ν 2 , +∞[. In this case, q(ν) is strictly increasing on ] -∞, ν 2 ] and strictly decreasing on [ν 2 , +∞[, and thus has a maximum at ν = ν 2 which is the unique extremum. If g(ν 2 ) ≤ 0, g(ν) and hence q (ν) annihilate and change signs once for a certain ν 1 in ] -∞, ν 1 [ and then q has a maximum at ν = ν 1 which will be its unique extremum. If g(ν 1 ) < 0 and g(ν 2 ) > 0, g(ν) and q (ν) annihilate and change sign 3 times, once in ] -∞, ν 1 [, once in ]ν 1 , ν 2 [, and once in [ν 2 , +∞[, and thus q has three extrema. We finally need to compute g(ν 1 ) and g(ν 2 ). Remember that ν 1 and ν 2 are the roots of g thus of

1 - 1 2πσ 2 ξ 2 = tanh 2 (2πσ 2 ξν - log(A) 2 )
and therefore

ν 1 = 1 2πσ 2 ξ log(A) 2 + arctanh(-1 - 1 2πσ 2 ξ 2 ) ν 2 = 1 2πσ 2 ξ log(A) 2 + arctanh( 1 - 1 2πσ 2 ξ 2 )
and thus as α = π 2 σ(ξ 2 -ξ 1 ) = √ 2πσξ, we may write:

g(ν 1 ) = - log(A) 4πσ 2 ξ 2 - arctanh(-1 -1 2πσ 2 ξ 2 ) 2πσ 2 ξ 2 -1 - 1 2πσ 2 ξ 2 = - log(A) 2α 2 - arctanh(- √ α 2 -1 α ) α 2 - √ α 2 -1 α = - log(A) 2α 2 - log( α- √ α 2 -1 α+ √ α 2 -1 ) 2α 2 - √ α 2 -1 α g(ν 2 ) = - log(A) 2α 2 - log( α+ √ α 2 -1 α- √ α 2 -1 ) 2α 2 + √ α 2 -1 α .
From this we deduce that g(ν 1 ) < 0 and g(ν 2 ) > 0, and thus q has three extrema, if and only if:

| log(A)| < 2α α 2 -1 -log( α + √ α 2 -1 α - √ α 2 -1 ).

B. Proof of Proposition III.2

The derivative of |V h f ( tk , η)| 2 is null if η satisfies:

(ξ 1 -η)Ae -πσ 2 (η-ξ1) 2 + (η -ξ 2 )e -πσ 2 (η-ξ2) 2 Ae -πσ 2 (η-ξ1) 2 -e -πσ 2 (η-ξ2) 2 = 0. (24)
Looking for a a solution in [ξ 1 , ξ 2 ], we remark that the derivative is null if Ae -πσ 2 (η-ξ1) 2 -e -πσ 2 (η-ξ2) 2 = 0, which leads to

η = log A 2πσ 2 (ξ 1 -ξ 2 ) + ξ 1 + ξ 2 2 . ( 25
)
This extremum corresponds to a zero of |V h f ( tk , η)| 2 , and is thus a minimum of |V h f ( tk , η)| 2 . Now, considering the same formalism as previously, outside of the interval ]ξ 1 , ξ 2 [, the extrema satisfy

η -ξ 1 η -ξ 2 = Ae -2πσ 2 (ξ1-ξ2)(η-ξ 1 +ξ 2 2
) .

Putting

h 1 (η) = η-ξ1 η-ξ2 and h 2 (η) = Ae -2πσ 2 (ξ1-ξ2)(η-ξ 1 +ξ 2 2
) ,

we remark that h 1 (η) < 0 and h 2 (η

) > 0. On in- terval ] -∞, ξ 1 [, as lim η→-∞ h 1 (η) = 1, h 1 (ξ 1 ) = 0, lim η→-∞ h 2 (η) = 0 and h 2 (ξ 1 ) = Ae -πσ 2 (ξ1-ξ2)

C. Proof of Proposition III.3

Let us consider A = 1, for which we put:

l(η) = e -2πσ 2 (η-ξ1) 2 + e -2πσ 2 (η-ξ2) 2 +2e -πσ 2 [(η-ξ1) 2 +(η-ξ2) 2 ] cos(2π(ξ 2 -ξ 1 )t), (27)
and then define γ = cos(2π(ξ 2 -ξ 1 )t). Considering the same change of variables as previously, namely η = ν + ξ1+ξ2 2 and reclling ξ = ξ2-ξ1 2 , we define:

q(ν) := l(ν + ξ 1 + ξ 2 2 ) = e -2πσ 2 (ν+ξ) 2 + e -2πσ 2 (ν-ξ) 2 +2e -πσ 2 [(ν+ξ) 2 +(ν-ξ) 2 ] γ, (28)
whose derivative reads:

q (ν) = 8πσ 2 e -2πσ 2 ν 2 e -2πσ 2 ξ 2 -ν(cosh(4πσ 2 ξν) + γ) + ξ cosh(4πσ 2 ξν) = 8πσ 2 e -2πσ 2 ν 2 e -2πσ 2 ξ 2 (cosh(4πσ 2 ξν) + γ) -ν + ξ sinh(4πσ 2 ξν) cosh(4πσ 2 ξν) + γ , (29)
which has the sign of the odd function k(ν) = -ν + ξ sinh(4πσ 2 ξν) cosh(4πσ 2 ξν)+γ (so we study it only for positive ν). Computing the derivative of this function we get:

k (ν) = -1 + 4πσ 2 ξ 2 cosh(4πσ 2 ξν)(cosh(4πσ 2 ξν) + γ) -sinh 2 (4πσ 2 ξν) (cosh(4πσ 2 ξν) + γ) 2 = -1 + 4πσ 2 ξ 2 1 + β cosh(4πσ 2 ξν) (cosh(4πσ 2 ξν) + γ) 2 , (30)
which has the same sign as r(ν) = -(cosh(4πσ 2 ξν)+γ) 2 +4πσ 2 ξ 2 (1+γ cosh(4πσ 2 ξν)).

Differentiating again, we get

r (ν) = -2(cosh(4πσ 2 νξ) + γ)4πσ 2 sinh(4πσ 2 νξ) +16π 2 ξ 3 σ 4 ξ 2 sinh(4πσ 2 νξ)) = 8πξσ 2 sinh(4πσ 2 νξ) (-cosh(4πσ 2 νξ) -γ + 2πσ 2 ξ 2 ). (31)
which has the sign of s(ν) = -cosh(4πσ 2 νξ) -γ + 2πσ 2 ξ 2 which is decreasing for ν > 0. As s(0) = -1-γ+2πσ 2 ξ 2 , we have two cases. The first one is when

-1-γ+2πσ 2 ξ 2 ≤ 0. As r(0) = -(1+γ) 2 +4πσ 2 ξ 2 (1+γ) = (1+α)(-1-γ+4πσ 2 ξ 2 ), if -1 -γ + 4πσ 2 ξ 2 ≤ 0, we have the following table of variations: ν 0 +∞ r (ν) - r(0) ≤ 0 r(ν) -∞ k (ν) - 0 k(ν) -∞
Note that from the above table of variations we deduce that the function q is decreasing for ν > 0, and q has a unique maximum at 0. Now if -1 -γ + 4πσ 2 ξ 2 > 0 still with -1γ + 2πσ 2 ξ 2 ≤ 0, we get the table of variations:

ν 0 ν 0 ν 1 ν 2 +∞ r (ν) - r(0) > 0 r(ν) -∞ k (ν) + 0 - k(ν) 0 -∞ q (ν) + 0 - q(ν) q(0) -∞
and so in that case, we obtain that q has three extrema.

The second case corresponds to -1 -γ + 2πσ 2 ξ 2 > 0, then -1-γ + 4πσ 2 ξ 2 > 0 then r(0) > 0, and we get the following table of variations:

ν 0 ν 0 ν 1 ν 2 +∞ r (ν) + 0 - r(ν) r(0) > 0 -∞ k (ν) + 0 - k(ν) 0 -∞ q (ν) + 0 - q(ν) q(0) -∞ Conclusion:
The function q has three extrema if 1 + γ < 4πσ 2 ξ 2 , which corresponds to the sought condition, and 1 otherwise.

D. Proof of Proposition III.4

Let us consider the function:

l(η) = A 2 e -2πσ 2 (η-ξ1) 2 + e -2πσ 2 (η-ξ2) 2 +2Ae -πσ 2 [(η-ξ1) 2 +(η-ξ2) 2 ] cos(2π(ξ 2 -ξ 1 )t), (32)
and then put γ = cos(2π(ξ 2 -ξ 1 )t) Considering the same change of variables as previously, namely η = ν + ξ1+ξ2 2 and putting ξ = ξ2-ξ1 2 , we may write:

f 1 (ν) =l(ν + ξ 1 + ξ 2 2 ) =e -2πσ 2 (ν 2 +ξ 2 ) (A 2 e -4πσ 2 ξν + e 4πσ 2 ξν + 2Aγ). (33)
We then make a new change of variables by putting ν = log(A)+µ 4πσ 2 ξ , which enables us to write:

f 1 (ν) =f 2 (µ) =e -2πσ 2 log(A)+µ 4πσ 2 ξ 2 +ξ 2 (Ae -µ + Ae µ + 2Aγ) =e -2πσ 2 log(A)+µ 4πσ 2 ξ 2 +ξ 2 2As(cosh(µ) + γ). (34)
The derivative of f 2 reads:

f 2 (µ) =e -2πσ 2 log(A)+µ 4πσ 2 ξ 2 +ξ 2 2A - log(A) + µ 4πσ 2 ξ 2 (cosh(µ) + γ) + sinh(µ) (35)
which has the sign of:

g(µ) = - log(A) + µ 4πσ 2 ξ 2 +
sinh(µ) cosh(µ) + γ assuming γ = -1, and recalling that α = π 2 σ(ξ 2 -ξ 1 ) = √ 2πσξ, we rewrite:

g(µ) = - log(A) + µ 2α 2 + sinh(µ) cosh(µ) + γ .
Differentiating g we get

g (µ) = - 1 2α 2 + cosh(µ)(cosh(µ) + γ) -sinh 2 (µ) (cosh(µ) + γ) 2 = - 1 2α 2 + 1 + γ cosh(µ) (cosh(µ) + γ) 2 (36)
which has the same sign as

h(cosh(µ)) = - (cosh(µ) + γ) 2 2α 2 + 1 + γ cosh(µ) = - cosh(µ) 2 2α 2 + γ(1 - 1 α 2 ) cosh(µ) + 1 - γ 2 2α 2 = 1 α 2 (- cosh(µ) 2 2 + γ(α 2 -1) cosh(µ) + α 2 - γ 2 2
).

The term inside the parentheses is a second order polynomial in cosh(µ) whose discriminant reads:

∆ = γ 2 (α 2 -1) 2 + (2α 2 -γ 2 ) = α 2 (γ 2 (α 2 -2) + 2) > 0.
The roots of this polynomial are denoted by X 1 and X 2 with X 1 < X 2 , and we know that h

(X) < 0 if X ∈ ] -∞, X 1 [ ]X 2 , +∞[ and h(X) > 0 if X ∈]X 1 , X 2 [.
We now need to know the location of cosh(µ). To do this, let us first compute h(1) = (γ + 1)(1 -γ+1 2α 2 ) which has the same sign as α 2 -γ+1 2 . Let us assume α ≤ γ+1 2 , as h(1) ≤ 0 then 1 belongs to ] -∞, X 1 ] or [X 2 , ∞[. To determine to which interval 1 belongs to we compare it to X1+X2 2 = γ(α 2 -1) < 1. Thus 1 belongs to ]X 2 , +∞[. Finally, as for all µ ∈ R cosh(µ) ≥ 1, h(µ) and thus g (µ) are negative. From this study we deduce the following table of variations:

µ -∞ µ 1 +∞ +∞ g(µ) -∞ f 2 (µ) + 0 - f 2 (µ)
and that the function l(η) has a unique extremum which is a maximum when α ≤ γ+1 2 . If we now assume that α > γ+1 2 , then h(1) > 0 and thus 1 belongs to ]X 1 , X 2 [. As cosh(µ) ∈ ]X 1 , X 2 [ is equivalent to |µ| < arccosh(X 2 ), this means that h(cosh(µ)) < 0 if µ ∈] -∞,arccosh(X 2 )[ ] arccosh(X 2 ), +∞[, and h(cosh(µ)) is positive in the opposite case. We thus get the following table of variations:

µ -∞ -arccosh(X 2 ) arccosh(X 2 ) +∞ g (µ) - 0 + 0 - +∞ g(µ)
-∞ g(µ) vanishes with a change of sign only once at some µ = µ 1 if and only if g(arccosh(X 2 )) ≥ 0 or g(arccosh(X 2 )) ≤ 0.

In this case, we deduce that µ -∞

µ 1 +∞ f 2 (µ) + 0 - f 2 (µ)
meaning l(η) has a unique maximum in that case (which is a maximum). If g(arccosh(X 2 )) < 0 and g(arccosh(X 2 )) > 0, g(µ) will vanish with a change of signs three times at some µ = µ 1 , µ 2 and µ 3 respectively. We therefore deduce the following table of variations for f 2 : µ -∞ µ 1 µ 2 µ 3 +∞ g(µ)

+ 0 -0 + 0 -

f 2 (µ)
In such case, l(η) has three extrema: 2 maxima and a mimimum. Finally, to specify in which situations l has three extrema, as g(arccosh(X 2 )) = -log(A) + arccosh(X 2 )

2α 2 + X 2 2 -1 X 2 + γ g(-arccosh(X 2 )) = - log(A) -arccosh(X 2 ) 2α 2 - X 2 2 -1 X 2 + γ
we deduce that: g(arccosh(X 2 )) > 0 ⇔ log(A) <arccosh(X 2 )) + 2α 2 X 2 2 -1 X 2 + γ g(-arccosh(X 2 )) < 0 ⇔ -log(A) <arccosh(X 2 )) + 2α 2 X 2 2 -1 X 2 + γ .

(37)

From this we get that l(η) has three extrema if and only if α > with X 2 = γ(α 2 -1) + α γ 2 (α 2 -2) + 2.

E. Proof of Proposition III.5

In this section, we assume that -1 < γ ≤ 1. We know that f 2 (µ) has the same sign as:

g(µ) = - log(A) + µ 4πσ 2 ξ 2 +
sinh(µ) cosh(µ) + γ we denote this function by g(µ, γ). We can rewrite: g(µ, γ) = -log(A) + µ 2α 2 + sinh(µ) cosh(µ) + γ .

We know from the previous proposition that the function f 2 has three extrema if and only if α > 1 + γ 2 and g(arccosh(X 2 ), γ) < 0 and g(arccosh(X 2 ), γ) > 0

Let us show that this is equivalent to: ∃ y 0 < 0 g(y 0 , γ) < 0 and y 0 > 0 g(y 0 , γ) > 0

This condition is a necesssary condition since the values y 0 =arccosh(X 2 ) and y 0 = arccosh(X 2 ), satisfy the property. Conversely, let us assume that there exist y 0 < 0 and y 0 > 0 such that g(y 0 , γ) < 0 and g(y 0 , γ) > 0. As lim µ→+∞ g(µ, γ) = -∞ and lim µ→-∞ g(µ, γ) = +∞, this implies that g(., γ) (and also f 2 (µ) ) passes through 0 and changes signs three times, respectively on the intervals ] -∞, y 0 [, ]y 0 , y 0 [, and on ]y 0 , +∞[, and thus f 2 has three extrema.

Let us now show that, if the function f 2 has three extrema for γ = 1, it also has three extrema for γ ∈] -1, 1]. Based on the above characterization of the existence of three extrema based on the study of the function g, let us assume that there exist y 0 < 0 and y 0 > 0 such that g(y 0 , 1) and g(y 0 , 1) > 0. But whatever γ ∈]-1, 1] one has ∀µ ∈ R 1 cosh(µ) ≥ 1 cosh(µ)+1 , which means that ∀µ ≥ 0, sinh(µ) cosh(µ) + γ ≥ sinh(µ) cosh(µ) + 1 and ∀µ ≤ 0, sinh(µ) cosh(µ) + γ ≤ sinh(µ) cosh(µ) + 1 .

Finally, we deduce that ∀µ ≥ 0, g(µ, γ) ≥ g(µ, 1) and ∀µ ≤ 0, g(µ, γ) ≤ g(µ, 1).

This means in particular that g(y 0 , γ) ≥ g(y 0 , 1), and that g(y 0 , γ) ≤ g(y 0 , 1). So we deduce from this that g(y 0 , γ) > 0 and g(y 0 , γ) < 0, and thus whatever γ ∈] -1, 1], the function g satisfies the condition for f 2 to have three extrema. So finally, we may conclude that the function f 2 has three extrema whatever γ ∈ [-1, 1], if and only if α > 1 and | log(A)| < -2 arccosh(α) + 2α α 2 -1.

Fig. 1 :

 1 Fig. 1: (a): spectrogram of two parallel pure harmonic modes with the same amplitude, when there exists three extrema at each time instant (ridges are thus well defined); (b): spectrogram of two parallel pure harmonic modes exhibiting TFBs (ridges associated with each mode can no longer be defined); (c): same as (a) but when the pure harmonic modes have different amplitude (three extrema are still present but the symmetry no longer exists).

Fig. 2 :

 2 Fig. 2: (a): Spectrogram of two parallel linear chirps with the same amplitude, when there exists three extrema at each time instant in the direction D, the upper and lower ridges are also superimposed; (b): Spectrogram of two parallel linear chirps when TFBs are present (three extrema no longer exist in the direction D for each translation factor B k , and ridge extraction is no longer possible); (c): same as (a) but when the linear chirps have different amplitude (three extrema are still present along the direction D).

Fig. 3 :

 3 Fig. 3: (a): spectrogram of two interfering modes with DTSZ superimposed, the white circles correspond to vertices v such that I(v, 30) ≥ 1 (input SNR = 20 dB); (b): same as (a) except the white circles correspond to vertices v such that I(v, 30) ≥ 4; (c): same as (a) with I(v, 30) ≥ 4; (d): spectrogram of three interfering modes with DTSZ superimposed, the white circles correspond to vertices v such that I(v, 30) ≥ 1 (input SNR = 20 dB); (e): same as (d) with I(v, 30) ≥ 4; (f): same as (d) with I(v, 30) ≥ 5.

Fig. 4 :

 4 Fig. 4: (a): spectrogram with DTSZ superimposed, the white circles correspond to the zeros v satisfying I(v, 30) ≥ 5 (input SNR equals 20 dB); (b): same as (a) except the input SNR equals 10 dB; (c): same as (a) except the input SNR equals 5 dB; (d): DTSZ on which the mask corresponds to |V g f (t, η)| ≥ βγ, with β = 8, the red circles then correspond to I signal (v, 30, 8) = 1 (input SNR equals 20 dB); (e): same as (d) except the input SNR equals 10 dB; (f): same as (d) except the input SNR equals 5 dB and β is set to 4.

Fig. 5 :

 5 Fig. 5: Proportion of zeros corresponding to I(v, P e ) ≥ 5 when P e varies and for the signal whose spectrograms are depicted in Fig. 4 (a), (b) and (c).

Fig. 6 :

 6 Fig. 6: Number of zeros kept corresponding to I signal (v, 30, β) = 1 when β varies and for the signal whose spectrograms are depicted in Fig. 4 (a), (b) and (c).

1+γ 2

 2 and then if | log(A)| <arccosh(X 2 )) + 2α 2

  2 , where . means we consider this variable, we get the following result.

	Proposition III.1. |V h f (t k , .)| 2 has three extrema if and only if α := π 2 σ(ξ 2 -ξ 1 ) > 1 and

  The proof is available in Appendix B. This can be generalized to other time instants than t k and tk in the case A = 1, since we have the following proposition.

		300		STFT			300		STFT			300		STFT
						upper ridge						Minima					upper ridge
		280				lower ridge Minima		280				Maxima Saddle points		280			lower ridge Minima
						Maxima											Maxima
	frequency	220 240 260				Saddle points	frequency	220 240 260					frequency	220 240 260			Saddle points
		200						200						200
		0.4 180	0.42	0.44	0.46	0.48	0.5	0.4 180	0.42	0.44	0.46	0.48	0.5	0.4 180	0.42	0.44	0.46	0.48	0.5
				time						time						time
																			)
										Then, we have the following proposition regarding the number
										of extrema of |V h f ( tk , .)| 2 :
										Proposition III.2. |V h f ( tk , .)| 2 has a unique minimum in-
										side the interval [ξ 1 , ξ 2 ], a unique maximum on the interval
										]ξ 2 , +∞[, and a unique maximum on the interval ] -∞, ξ 1 [.
										Proposition III.3. If A = 1, |V h f (t, .)| 2 has three extrema if and only if α = π 2 σ(ξ 2 -ξ 1 ) satisfies:
													α >	1 + cos(2π(ξ 2 -ξ 1 )t) 2		(8)


Jian Hu
email: hujian@whut.edu.cn

Gangyan Li
email: gangyanli@whut.edu.cn

Yeqiong Song
email: song@loria.fr

Research on Calculation Method of Period and Deadline of Frame in Automotive Electronic and Information Integrated Control System

Keywords: Automotive electronic, In-vehicle network, Networked control, Temporal characteristic, Real-time, GCRT

With the increase of the complexity of automotive control, it is necessary to develop AEIICS (Automotive Electronic and Information Integrated Control System) for exchanging information between ECUs (Electronic Control Units) based on in-vehicle network technology. We named the parameter (or variable) contained in and transmitted by frame as signal. Period and deadline are basic temporal characteristics for signal and frame. Designers mainly focus on directly designing the period and deadline of frame but omit the temporal characteristic relationship between signal and frame. In this paper, we researched on the calculation method of period and deadline of frame according to the temporal characteristic of signal. Based on the period of signal, we proved the period of frame is the minimal period of signals composing this frame. Considering the possible time offset between the production time of data of signal and the first transmission time of frame containing this data, we deduced the calculation formula of deadline of frame by applying GCRT (Generalized Chinese Remainder Theorem).

INTRODUCTION

Nowadays, with the requirements of security, comfort and environment protection, automobiles are more and more equipped with ECUs (Electronic Control Units) to implement different control functions. These functions require the use of real-time and reliable exchange of information between different ECUs, sensors and actuators. But it is impossible to realize this exchange of information through the traditional point-to-point links which would induce a disproportionate length of cable, an increase of cost and production time, reliability problems and other drawbacks. Thus, it is essential to develop AEIICS (Automotive Electronic and Information Integrated Control System) with in-vehicle network for sharing information and implementing real-time and correlative control between ECUs [START_REF] Hu | Design and Application of SAE J1939 Communication Database in City-Bus Information Integrated Control System Development[END_REF]. At the beginning of 1980s the engineers of the automotive manufacturers assessed the existing field bus systems for their using in vehicles [START_REF] Domínguez | Communication Networks for Vehicular Electronic Devices[END_REF]. With the development of in-vehicle network technology, the typical in-vehicle network protocols are LIN [START_REF]LIN specification package, revision 2.0[END_REF], CAN [START_REF] Iso | Road vehicles-low speed serial data communication-part 2: low speed controller area network[END_REF]- [START_REF] Iso | Road vehicles-interchange of digital information-controller area network for high-speed communication[END_REF], TTCAN [START_REF] Müller | Fault tolerant TTCAN networks[END_REF]- [START_REF] Iso | Road vehicles-controller area network (CAN)-part 4: time-triggered communication[END_REF], TTP/C [START_REF]Time-triggered protocol TTP/C, high-level specification document, protocol version 1.1[END_REF], FlexRay [START_REF]FlexRay communication system, protocol specification, version 2.0[END_REF], and MOST [START_REF]MOST specification revision 2.3[END_REF] etc. which have different characteristics and are applied in different automotive domains. All kinds of status information and control information are basic parameters for well-working of AEIICS. In this paper, we named the parameter (or variable) transmitted through in-vehicle network as signal. The instance of signal, termed as the data of this signal, identifies different parameter value in AEIICS. For example, engine speed is a signal. 1000rpm is one data of engine speed and a parameter value of engine control system. In a word, the basic function of AEIICS is to share datum of signals between different ECUs based on in-vehicle network technology. Signal transmission in AEIICS is implemented by frame which is the information unit exchanged in data link layer of ISO/OSI (International Standardization Organization/Open System Interconnect) Reference Model. Before sending signals, data link layer packages signals into frame with specific in-vehicle network protocol format. Every frame includes not only the data of signal but also some additional control information such as synchronization information, address information and error control information etc. Currently, research of AEIICS focuses on network message scheduling, system performance analysis and system development method etc [START_REF] Navet | Trends in automotive communication systems[END_REF]. Period and deadline are primary parameters to describe temporal characteristics for signal and frame in AEIICS. In system design process, period and deadline of signal depend on sampling period of sensor, control algorithm and experience of system engineer etc., the period and deadline of frame depend on the temporal characteristic of signal. However, in current research and design of AEIICS, researchers and engineers mainly focus on directly designing the period and deadline of frame but omit the temporal characteristic relationship between signal and frame (such as SAE J1939 [START_REF] Sae | Surface vehicle recommended practice, SAE J[END_REF]). Because of directly specifying the period and deadline of frame in system design, it couldn't ensure the temporal characteristic of signal and there would be system security problem induced by this kind of design method. In this paper, we mainly research on the calculation method of period and deadline of frame in AEIICS. This paper is organized as follows: In section 2, we described the characteristics of signal and frame of AEIICS with tuples respectively. In section 3, we defined the conditions with which the period of frame should be satisfied. We also presented and proved the calculation formula of period of frame according to the periods of signals. By illustrating the deadline of frame isn't always the minimal deadline of signals composing this frame, in section 4, we deduced the calculation formula of deadline of frame according to the temporal characteristics of signals. In section 5, we calculated the period and deadline of frame containing lighting control signals as a case study. Finally, in section 6 we concluded and presented future work.

FORMAL DESCRIPTION FOR CHARACTERISTICS OF SIGNAL AND FRAME IN AEIICS

Regarding temporal characteristic as the main description object, each signal ( ) in AEIICS is characterized by a tuple ( , , , ):

i s + Ζ ∈ i • i N • i T • i D • i C • i N is the node (i.e. ECU) sending . i s • i
T is the production period of . In practice the clocks of the nodes are not synchronized but it is assumed that the first production time of all signals on the same node takes place at the same time. i s • i D is the deadline of . It is the maximal interval between the production of on the sending node and its reception by all consumers. Usually, it is assumed that the deadline of is equal to its period: .

i s i s i s = • i D • i T • i
C is the size of in bits with the assumption that is always smaller than or equal to the maximal data field size of frame (i.e. there is no segmentation of a signal). For CAN data frame, the maximal signal size is 64 bits because the maximal data field size of CAN data frame is 8 bytes.

CALCULATION METHOD OF PERIOD OF FRAME BASED ON TEMPORAL CHARACTERISTICS OF SIGNALS

Problem description

It is assumed that frame consists of signals { , , …, }( ). The periods of and ( ) are ( ) and (

k f k s 1 k s 2 k n s + Ζ ∈ n k, k f k i s + Ζ ∈ i k T k T + Ζ ∈ • i T • i T + Ζ ∈
) respectively. As shown in Fig. 3.1, the calculation problem of the period of can be described as follows: k f One finds out the maximal value of which ensures that the data of produced at

k i k T k i s * • i T ) ( Ν ∈ i k can be transmitted before ) 1 ( + i k * • i T for ∀ k i s ∈{ , , …, }. k s 1 k s 2 k n s • * + i i T k ) 1 ( • * i i T k t t L L L L k j T k * k f k i s 0 0 Fig.3.1 When is transmitted with , the data of produced at k i k f k T k i s * • i T ) ( Ν ∈ i k can be transmitted before ) 1 ( + i k * • i T .

Calculation method of period of frame

k f k s 1 k s 2 k n s + Ζ ∈ n k, k f k i s + Ζ ∈ i k T k T + Ζ ∈ • i T • i T + Ζ ∈ k s 1 k s 2 k n s k s min k s 1 k s 2 k n s = • min T • i T k i s ∈ k s 1 k s 2 k n s k f k
s min Before deducing the calculation method of period of frame, we firstly introduce one of the divisibility theorems. Theorem 1: Assuming b and b , any integer a can be uniquely formulized as follows:

Ζ ∈ 0 ≠ a = b * q+r where q, r and Ζ ∈ b r < ≤ 0 .
For { , , …, }, applying theorem 1, the duration of between

∀ k i s ∈ k s 1 k s 2 k n s k i s • * i i T k ) ( Ν ∈ i k
and is formulized by:

• * + i i T k ) 1 ( = ∆ k i q * r (3.1) + k T where q, r and Ζ ∈ k T r < ≤ 0 . Setting , k T min • = T because = ∆ k i k i * (3.2) - * + • i i T k ) 1 ( = • i T • i T and , • i T min • ≥ T thus: = ∆ k i q * r (3.3) min • T +
where q, r and q , .

Ζ ∈ 1 ≥ < ≤ r 0 min • T Therefore,
is transmitted for q (q and q ) times during the duration [ , ), i.e. contained by is transmitted for q times. It ensures that the data of produced at can be transmitted by before .

k f Ζ ∈ 1 ≥ • * i i T k • * + i i T k ) 1 ( k i s k f k i s • * i i T k ) ( Ν ∈ i k k f • * + i i T k ) 1 (
It has proved can be the period of . We will prove is the maximal allowed period of by using reduction to absurdity as follows. ) (

+ Ζ ∈ ε k f k T min • = T ε + .
According to formula (3.1), we obtain:

= ∆ k i q * ( min • T ε + ) r +
where q, r Ζ ∈ and . ) ( 0

min ε + < ≤ • T r Setting , • i T min • = T according to formula (3.2): k min ∆ min • = T = q * ( min • T ε + ) + r (3.4)
where q, r Ζ ∈ and q 0 ≤ , . ) ( 0

min ε + < ≤ • T r Therefore, when is transmitted with period ( k f min • T ε + ), it can't ensure that the data of produced at k i s • * i i T k ) ( Ν ∈ i k
will be transmitted by before

during the duration [ , )
.

k f • * + i i T k ) 1 ( • * i i T k • * + i i T k
) 1 ( Thus, the assumption is wrong. Finally, we obtain the maximal allowed period of , i.e. the period of defined by us is:

k f k f = k T min • T (3.5)

CALCULATION METHOD OF DEADLINE OF FRAME BASED ON TEMPORAL CHARACTERISTICS OF SIGNALS

Problem description

When it is assumed that the transmission time of frame containing signals { , , …, }( ) is synchronized with the production time of , we have known that the period of is , the minimal period of signal set { , , …, }. So should we conclude that the deadline of is the minimal deadline of the set { , , …, }?

k f k s 1 k s 2 k n s + Ζ ∈ n k, k s min k f min • T k s 1 k s 2 k n s k f k s 1 k s 2 k n s
Because there is possible offset [START_REF] Saket | Frame packing algorithms for automotive applications[END_REF] between the production time of the data of signal and the first transmission time of the frame containing this data, the deadline of is not the minimal deadline of signal set { , , …, } .

f k i s ∈ k s 1 k s 2 k n s k D k i s k i s k f k f
In order to calculate the deadline of , we firstly present and prove the mathematic formula between the production times of two periodic signals.

k f

This states that the interval between the production times of any two datum of two periodic signals always equals to the integer multiple of the great common division of their periods if these two periodic signals have the same initial time.

Relationship between the production times of two periodic signals GCRT (Generalized Chinese Remainder Theorem):

For any , a common residue equation set )) gcd( (mod

n j i ≤ < ≤ 1 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ≡ ≡ ≡ ) (mod ) (mod ) (mod
j i j i m m a a , ≡ (4.1)
When n=2, the common residue equation set

⎩ ⎨ ⎧ ≡ ≡ ) (mod ) (mod 2 2 1 1 m a x m a x
has a solution if and only if

)) , gcd( (mod

2 1 2 1 m m a a ≡ (4.2)
Given two periodic signals and with periods and , and the production time and respectively, it is assumed that they have the same initial time. We will deduce the relationship between and using GCRT.

1 s 2 s • 1 T • 2 T ) , ( 2 1 + • • Ζ ∈ T T 1 p 2 p ) , ( 2
Ν ∈ ∃ 2 1, k k 1 p 2 p • * = 1 1 1

T k p

Because there is possible time offset between the production time of the data of and the first transmission time of containing this data.

k i s k f • * = 2 2 2 T k p Thus: 0 mod 1 1 = • T p
Thus the deadline of is defined as the minimal value of the difference between and , i.e.:

k f • i D ) ( max • i k T T Offset , 0 mod 2 2 = • T p i.e.: }} { | ) ( min{ 2 1 max k n k k k i i k i k s s s s T T Offset D D , , , , L ∈ - = • • ⎪ ⎩ ⎪ ⎨ ⎧ ≡ ≡ • • ) (mod 0 ) (mod 0 2 2 1 1 T p T p Because = , k T min • T thus: Thus the common residue equation set }} { | ) ( min{ 2 1 min max k n k k k i i i k s s s s T T Offset D D , , , , L ∈ - = • • • . ⎪ ⎩ ⎪ ⎨ ⎧ ≡ ≡ • • ) (mod ) (mod 2 2 1 1 T p x T p x
Thus the key of calculating the deadline of is to calculate .

k f ) ( max • i k T T Offset ,
has solution and the solution is x=0.

In order to calculate , it is assumed that is the production time of the data of , is the first transmission time of containing this data. Then equals to .

) ( max

• i k T T Offset , • i t k i s k t k f ) ( • i k T T Offset , - k t • i t
According to GCRT, we know:

)) , gcd( (mod

2 1 2 1 • • ≡ T T p p i.e. 2 1 2 1 | ) gcd( p p T T - • • , thus q , ∃ Ζ ∈ ) gcd( 2 1 2 1 • • * = - T T q p p , (4.3) t t k f L L L L ) ( • i k T T Offset , L L L L ) ( max • i k T T Offset , ) ( • i k T T Offset , ) ( • i k T T Offset , k i s 0 0 • * i i T k • * + i i T k ) 1 ( • * + i i T k ) 2 ( • * + i i T k ) 3 ( k j T k * + ) 1 ( k j T k * + ) 2 ( k j T k * + ) 3 ( k j T k * + ) 4 ( k j T k * + ) 5 ( k j T k * Fig.4.1
: time offset between the production time of the data of and the first transmission time of containing this data; : the maximal value of .

) , (

• i k T T Offset k i s k f ) , ( max • i k T T Offset ) , ( • i k T T Offset Because , min • = T T k thus: k t = min k * min • T ) ( min Ν ∈ k • i t = i k * • i T ) ( Ν ∈ i k Thus: - k t = • i t min k * - • min T i k * • i T ,
• • * = - T T q p p , , i.e.: ) gcd( 2 1 2 2 1 1 • • • • * = * - * T T q T k T k , Setting , , • • = 1 min T T • • = 2 T T i thus: min k * - • min T i k * • i T ) gcd( min • • * = i T T q , ( Ν ∈ q
) Thus the key of calculating the maximal value of T is to calculate the maximal value of q. We calculate the maximal value of q as follows:

Because ≤ 0 min k * - • min T i k * • i T < , i.e.: min • T ) gcd( 0 min • • * ≤ i T T q , < min • T thus: ) gcd( min min • • • < i T T T q ,
And because ,

min min | ) gcd( • • • T T T i , thus: 1 ) gcd( ) ( max min min - = • • • i T T T q , Thus: ) ( max • i k T T Offset , ) ( max min min • • * - * = i i T k T k )) gcd( ) max( min • • * = i T T q , ) gcd( ) max( min • • * = i T T q , ) gcd( ) 1 ) gcd( ( min min min • • • • • * - = i i T T T T T , , = ) gcd( min min • • • - i T T T , Thus: k D =min{ } - • i D )} ( max • i k T T Offset , )} gcd( { min min min • • • • + - = i i T T T D ,
Therefore, the deadline of :

k f )} gcd( { min min min • • • • + - = i i k T T T D D , (4.4) If = , • i D • i T )} gcd( { min min min • • • • + - = i i k T T T T D , (4.5)

CASE STUDY

There are many lighting control signals in AEIICS. We usually configure all these signals into one frame which is transmitted through in-vehicle network. Because of different importance, there are different sampling period and real-time request for each lighting control signal. As shown in Table 5.1, it is different for the periods and deadlines of these signals. Here we assume that the deadline of lighting control signal is equal to its period. According to the culculaiton results, we know the deadline of frame isn't always equal to the minimal deadline of signals. It is essential to obtain the deadline of frame with formula (4.4) because it will affect on the analysis of system schdulability and the selection of nework transmission rate if we configure the deadline of frame as the the minimal deadline of signals composing this frame.

CONCLUSION

In this paper, we researched on the temporal characteristic relationship between signal and frame in automotive electronic and information integrated information control system. We proved the period of frame is the minimal period of signals composing this frame. By applying GCRT, we put forward the mathematic formula of the difference between the production times of two periodic signals. Based on this formula, we finally deduced the calculation formula of deadline of frame with the temporal characteristic of signals as input parameters. In SAE J1939 protocol, the periods and deadlines of some messages are specified, but others depends on the decision of system designer. The further research will focus on the conditions which make the period and deadline of frame equals to the minimal period and deadline of signals composing this frame and keep the period and deadline of frame when other signals are inserted into this frame.

NTCPf

  is the node (i.e. ECU) sending . is the transmission period of . The value of depends on the temporal characteristics of signals composing . is the size of in bits. It consists of the size of the data field and the size of the frame overhead. is the priority of for the priority MAC protocol.jThis paper focuses on how to obtain the period and deadline of frame according to the temporal characteristics of signals in AEIICS.

  It is assumed that is not the maximal allowed period of . Then there exists ( period of , i.e.

s

  The calculation problem of the deadline of can be described as follows: for ∀ { , , …, }, one finds out the value of which ensures that the response time of satisfies with the requirement of deadline of if the response time of satisfies with the requirement of deadline of .

  k

  transmission time of is synchronized with the production time of . As shown in Fig.4.1, is defined as the time offset between the production time of the data of and the first transmission time of containing this data (because the data of maybe be transmitted more than one time); is defined as the maximal value of . if and only if:

T

  Fig.4.2 The relationship between and ) ( • i k T T Offset ,

  40, 40, 40, 100, 50, 50, 50, 500, 100, 500, 500} = 40(ms)

Table 5 .

 5

		1 Periods and deadlines of lighting control
		signals
		Signal Name			Period (ms)	Deadline (ms)
	Running Light Command			500	500
	Alternate Beam Head Light Command		80	80
	Low Beam Head Light Command			80	80
	High Beam Head Light Command			80	80
	Front Fog Lights Command		100	100
	Right Turn Signal Lights Command			50	50
	Left Turn Signal Lights Command			50	50
	Back Up Light Command			50	50
	Marker Light Command			500	500
	Rear Fog Light Command		100	100
	Lighting Data Request Command			500	500
	Work Light Command			500	500
	min	• i T k i s ∈ k s 1	k s 2	s	k n

• T = min { | { , , …, }} k T = min

ACKNOWLEDGEMENT

We would like to give our acknowledgments to Hubei province, Wuhan city and Wuhan University of Technology, China because the research work reported here has been sponsored by "Hubei province digital manufacturing key laboratory open foundation" (SZ0615), "Wuhan city tackle key science and technology problem project" (200710421118) and "Excellent PhD. Dissertation Sustentation Fund of Wuhan University of Technology" (200704). We also want to thank Mathieu Grenier and Nicolas Navet, LORIA, France, for their helpful advices on this study.


Chung Shue Chen

Yanjun Li
email: yanjun.li@loria.fr

Ye-Qiong Song
email: song@loria.fr

An Exploration of Geographic Routing with k-Hop Based Searching in Wireless Sensor Networks

We explore the asymptotic performance of existing geographic routing with a utilization of k-hop neighborhood information. The reachability from source to sink improves as we use more information for routing decision and the average number of hops required decreases significantly from 1-hop to 2-hop searching, which indicates a potential tradeoff between performance enhancement and system complexity. As simple greedy geographic routing is insufficient in lossy wireless environment, we propose a new metric incorporating both advance in distance and link quality. Through simulations, we show its effectiveness over the conventional simple greedy method. The generalization to k-hop based routing and resulting performance are also presented. Simulation results show that with the multi-hop based searching, there is a good improvement in the number of transmissions required from source to sink, which also indicates potential improvement in routing delay and energy consumption in transmissions.

I. INTRODUCTION

Geographic routing was originally developed for packet radio networks in the 1980s [START_REF] Takagi | Optimal transmission ranges for randomly distributed packet radio terminals[END_REF] but now receives renewed interest as a promising routing candidate for wireless ad hoc and sensor networks due to its attractive properties including stateless nature and low maintenance overhead [START_REF] Stojmenovic | Position based routing in ad hoc networks[END_REF]. In principle, the neighbor node which can forward the packet towards destination in the most progress is selected. Node location is identified through some localization techniques or GPS [START_REF] Doherty | Convex position estimation in wireless sensor networks[END_REF], [START_REF] He | Rangefree localization schemes for large scale sensor networks[END_REF]. It may also involve some information exchange among neighboring nodes. Generally speaking, this greedy method [START_REF] Karp | GPSR: greedy perimeter stateless routing for wireless networks[END_REF] is simple and effective in many cases despite some recovery schemes should be incorporated to handle possible network voids or deadlocks in the delivery [START_REF] Bose | Routing with guaranteed delivery in ad hoc wireless networks[END_REF]. However, it is worth pointing out that most studies in this topic are located on 1-hop neighborhood information based schemes which often address a quite limited optimization. The generalization to 2-hop or more generally k-hop based routing has not been profoundly investigated despite interesting and possible [START_REF] Lee | Efficient geographic routing in multihop wireless networks[END_REF]. An exploration will be helpful to provide a reference for future design and consideration, especially for the purpose in emerging real-time wireless sensor networks (WSN) [START_REF] Stankovic | Real-time communication and coordination in embedded sensor networks[END_REF] and some demanding applications [START_REF] Li | Real-time QoS support in wireless sensor networks: a survey[END_REF].

In the following, we will first conduct a case study of geographic routing in different levels of neighborhood information, says k-hop, where k = 1, 2, 3, • • • , and particularly a performance comparison to understand the potential improvement by the extra searching, while assuming the neighborhood information is available. Note that, when k = ∞, the result will naturally fall into the classical shortest path routing with global search which is to serve as a benchmark for comparison. In the beginning, a performance study is conducted in an ideal and deterministic WSN model with simple geographic greedy routing determined by the progress in distance. This aims to offer a basic understanding on the performance of k-hop based routing while isolated from the impacts of other factors.

In practice, the strategy of simply choosing the neighbor node which is closest to the destination or in greatest distance progress may suffer from low packet reception rate (PRR) [START_REF] Zuniga | Analyzing the transitional region in low power wireless links[END_REF] or poor link quality with highly variant channel conditions particularly during hops of large distance. This can lead to packet delivery failure and high retransmission cost [START_REF] Li | Analyzing and modeling of the wireless link for sensor networks[END_REF]. A routing metric beyond simple greedy method is thus necessary and favorable. To conduct the study in a lossy and dynamic WSN, we adopt the model used in [START_REF] Seada | Energyefficient forwarding strategies for geographic routing in lossy wireless sensor networks[END_REF] which also considers link quality transitional regions [START_REF] Zuniga | Analyzing the transitional region in low power wireless links[END_REF]. Meanwhile, a new routing metric which incorporates the proximity to destination and link reliability is proposed and generalized. Results show that the performance of the new metric outperforms the simple greedy method as expected practically. Furthermore, we investigate the performance in k-hop neighborhood based routing under the lossy model and report comparisons.

The rest of the paper is organized as follows. Some related work is discussed in Section II. In Section III, geographic routing in distance metric with k-hop searching, where k = 1, 2, 3, is studied and compared. The lossy WSN model and new metric are reported in Section IV which includes the corresponding simulation study. Finally, we conclude the paper in Section V.

II. RELATED WORK

For system simplicity, most existing routing protocols are based on 1-hop neighborhood information. However, it is promising that multi-hop information can lead to better performance in many aspects including routing, message broadcasting and channel access scheduling [START_REF] Lee | Efficient geographic routing in multihop wireless networks[END_REF], [START_REF] Stojmenovic | Loop-free hybrid single-path/flooding routing algorithm with guranteed delivery for wireless networks[END_REF]- [START_REF] Lou | On reducing broadcast redundancy in ad hoc wireless networks[END_REF]. Details of several distributed information exchange schemes and efficient algorithms in computing 2-hop neighborhood information for wireless ad hoc and sensor networks can be found in [START_REF] Čapkun | GPS-free positioning in mobile ad-hoc networks[END_REF]- [START_REF] Rajendran | Energyefficient collision-free medium access control for wireless sensor networks[END_REF]. In a network of n nodes, a complexity analysis is reported in [START_REF] Calinescu | Computing 2-hop neighborhoods in ad hoc wireless networks[END_REF] that every node can obtain the knowledge of 2-hop neighborhood by a total of O(n) messages, each of O(log n) bits.

For WSN, there are some geographic distance based routing metrics proposed in the literature. In [START_REF] Lee | Efficient geographic routing in multihop wireless networks[END_REF], a link metric called normalized advance (NADV) is proposed for wireless multihop networks. It selects neighbor nodes in tradeoff between link cost and proximity to destination. In general, the link cost can be defined by the application objectives such as delay, energy consumption or link quality. In [START_REF] Seada | Energyefficient forwarding strategies for geographic routing in lossy wireless sensor networks[END_REF], beyond conventional geographic greedy routing in maximum distance advance which may perform poorly in lossy wireless links, a study of routing strategies in distance-hop tradeoff is conducted. Results show that the product of PRR and the distance traversed towards destination is a routing metric for packet forwarding in general. Our work here is more related to them. The new routing metric proposed in Section IV is inspired by the concept of expected advance per transmission in [START_REF] Lee | Efficient geographic routing in multihop wireless networks[END_REF].

It should be noted that, to the best of our knowledge, this paper is the first attempt in exploring geographic routing with k-hop neighborhood information for WSN, where generally k ≥ 1, and giving a comparative study on asymptotic properties with different levels of hop information. Furthermore, for lossy WSN, new distance-PRR joint metrics are developed in generalization while taking into account multihop information for potential geographic routing schemes.

III. A STUDY OF SIMPLE GEOGRAPHIC GREEDY ROUTING

IN k-HOP NEIGHBORHOOD INFORMATION In this section, we extend the conventional geographic greedy routing to the utilization of k-hop neighborhood information to investigate the potential improvement when more information is allowed. In principle, the neighbor node which can forward the packet by the most progress is selected as the next hop forwarder. Given k-hop neighborhood information, the one in next forwarding which will lead to a node, after khop routing, closest to the destination is chosen. In generalization, it is chosen by its potential advantage in terms of k-hop transmissions identified by the neighborhood information in the corresponding range. The determination follows the same concept of preference in the proximity to destination.

As we know, in geographic routing, a selection of forwarding candidate only by the preference of shortest distance to destination can simply fall into a routing deadlock [START_REF] Stojmenovic | Loop-free hybrid single-path/flooding routing algorithm with guranteed delivery for wireless networks[END_REF]. In other words, the packet may be passed between two nodes at the dead end. This can be resolved by many different techniques. A simple and effective solution is to blacklist the visited end node and eliminate it from the forwarding candidate list since a loop-back occurs. Here, we follow this approach. Besides, in this Section, the link model adopted is simply a boolean disk, which means if the packet transmission is in the effective radio range, R, the delivery is considered successful. A lossy link model will be investigated later in Section IV.

In the performance comparison, the number of hops or transmissions required from source to destination (sink) in the packet forwarding is in our major concern as it reflects both the routing delay and energy consumed during transmissions. Here, we will not take into account MAC and queueing delay in order to focus on the issue of routing. The reachability to be reported soon is defined by the ratio of packets from source which can reach the sink finally. When we find the packet is always buffered at a specific node or possibly looped in an isolated region and thus cannot go into the right direction after a large number of hops, it will be considered as routing failure. This indicator can reflect the routing capability of a path searching scheme. However, significantly, it also depends on the node density and can help to indicate the level of network connectivity.

Nodes in the WSN are uniformly distributed in a geographical area of (200m, 200m), while the source and sink are fixed at the location of (30m, 30m) and (170m, 170m) respectively. Clearly shown in Fig. 1, as the number of nodes increases, the reachability from source to sink increases as well generally due to the resulting better network connectivity. We conduct two sets of simulations with different effective radio range, R, of 26m and 22m respectively, and plot in parallel for reference. As expected, when the number of nodes is insufficient, the reachability will be low even by ∞-hop searching due to disjoint clusters or low connectivity between source and sink. It is therefore unsuitable for network formation.

As shown in Fig. 1, when the number of nodes is above 150 in the case R = 26m, the reachability is at a satisfactory level approximately equal to or above 0.9. Reasonably, in the case R = 22m, a higher node density is required to reach this reachability. Detailed studies of critical density for wireless network connectivity can be found in [START_REF] Gupta | Critical power for asymptotic connectivity in wireless networks[END_REF], [START_REF] Krishnamachari | Critical Density Thresholds in Distributed Wireless Networks[END_REF]. Here, we will focus on the comparison among different k-hop searching results. As observable in Fig. 1, there is an improvement in reachability from 1-hop to 2-hop path searching and also an improvement from 2-hop to 3-hop searching by the fact that some routing difficulties can be earlier identified and better resolved, or even avoided, when more neighborhood and connectivity information is provided. However, generally speaking, the improvement is decreasing as the number of nodes is increasing. For example, the reachability tends to 1 in all the curves when the number of nodes goes to 300.

It should be noted that the number of hops required from source to sink is an important indicator from the point of view of end-to-end packet delay and energy consumption. This is at the core of our study. As shown in Fig. 2, when the number of nodes increases, the number of hops required decreases as there can be more and better forwarding choices to reach the sink with potentially less hops. An important observation here is that there is a significant improvement in the number of hops required from 1-hop searching to 2-hop searching, while the improvement from 2-hop to 3-hop searching and even more is relatively small. For example, we can clearly see this in the case R = 26m and the number of nodes equals to 200, in which the reachability is almost 1 as known from Fig. 1. Generally speaking, in this system, 2-hop searching is a good option in different k values by its attractive gain in the reduction of required number of hops. Meanwhile, the system

1-hop, R = 22m 2-hop, 22m 3-hop, 22m ∞-hop, 22m 1-hop, R = 26m 2-hop, 26m 3-hop, 26m ∞-hop, 26m 26m 22m
Fig. 2. Average number of hops required in k-hop searching complexity is concerned. Note that, as expected, the gain will get decreased when the number of nodes increases. In addition, by comparing the curves of R = 26m and R = 22m, as the effective radio range is reduced, under same node density, the number of hops required increases generally as in Fig. 2 since the network is now in a lower connectivity. In Fig. 1, the reachability in R = 22m drops significantly, especially in the region of relatively low node density, e.g. 125 to 175.

To have a better understanding on what kind of scenarios the 2-hop searching will take a smaller number of hops than the 1hop searching does, we study the node topology and resulting routing paths by details. Fig. 3 and 4 depict two typical examples observed. As shown in Fig. 3, both searching can find a routing path from source to sink eventually. However, the two searching will go into different paths since node X. Node V is chosen as the next forwarder in the 2-hop searching scheme instead of node Y, although Y is in fact closer to the During the possible deadlock, the 1-hop searching scheme will encounter more transmissions although the deadlock can be resolved.

sink. The reason is the 2-hop searching has identified that in X's 2-hop range, via V, node W is even closer to sink than node Z. So, node V is better than node Y from the point of view of 2-hop resulting progress. However, the 1hop searching will simply choose node Y. This effect often helps the 2-hop searching to drive into a shorter path than that of the 1-hop searching.

Fig. 4 shows an example of the 1-hop searching in which a possible routing deadlock occurs and is finally resolved. After forwarded to node Y from node X, the packet will be routed back to X. As now Y will be blacklisted, the packet is thus forwarded to node Z although Y is at a closer position to sink than Z. However, by the 2-hop searching, when the packet arrives X, in the next forwarder selection, the potential progress in 2-hop range will be checked. Clearly, Y is found as a dead end as it will be back to X. The 2-hop searching can thus bypass Y and forward the packet from X to Z directly despite Y is in fact closer to sink than Z. However, comparatively, the 1-hop searching will encounter more hops (transmissions) and consequently a delay in the deadlock scenario due to the difference in capability. It is worth pointing out that although the 1-hop and 2-hop searching may enjoy close reachability, the number of hops required from source to sink can be quite different generally.

IV. DESIGN OF A NEW ROUTING METRIC AND ITS GENERALIZATION IN LOSSY WSN MODEL

As reported in [START_REF] Zuniga | Analyzing the transitional region in low power wireless links[END_REF], a strategy of simply choosing the neighbor node which is closest to sink or in the greatest distance progress may suffer from low PRR or poor link quality simply due to the large distance propagation. This can result in packet delivery failure and high retransmission cost. In this section, we will propose a new routing metric which considers both the distance advance and link quality to be incorporated in the decision at each hop. Motivated by the study and observations in Section III, we generalize the new metric in the scenario of having k-hop neighborhood knowledge and compare the performance. The result can help to understand the potential if we can have more of the information for the routing decision in a more practical link model.

A. PRR Model

Before elaborating on the design, we describe the packet reception rate (PRR) model adopted, which is derived in [START_REF] Zuniga | Analyzing the transitional region in low power wireless links[END_REF] for lossy WSN links. It is built on experimental measures of practical systems with respect to existing statistical model of time-varying channels. The PRR, 0 ≤ p(d) ≤ 1, of a wireless link is expressible as:

p(d) = 1 - 1 2 exp - γ(d) 2 1 0.64 8f (1)
where d is the transmitter-receiver distance, γ(d) is the signalto-noise ratio (SNR), and f is the frame size which equals to 50 bytes including preamble, payload and CRC in following MICA2 motes [START_REF]Crossbow. Mica2 motes datasheet[END_REF] and the standard non-coherent FSK modulation. This model takes into account both distance-dependent path loss and log-normal shadowing in identifying the signal strength and randomness in wireless environment.

For transmitting power P t , the SNR, γ(d), is expressible as:

γ(d) dB = P t dB -P L(d) dB -P n dB (2)
where, by MICA2 radios, P t dB is set at 0 dBm, the noise floor P n dB is at -115 dBm, while the path loss P L(d) dB is modeled as:

P L(d) dB = P L(d 0 ) dB -10n log 10 (d/d 0 ) + X σ dB (3)
in which n is the path loss exponent, d 0 is the reference distance (1 meter), and X σ denotes the log-normal shadowing with zero mean and variance σ 2 . In reference to the above model and the discussions in [START_REF] Seada | Energyefficient forwarding strategies for geographic routing in lossy wireless sensor networks[END_REF], we will adopt the following three settings based on [START_REF] Sohrabi | Near ground wideband channel measurement in 800-1000 MHz[END_REF]. They represent different environments and characteristics.

B. Design of Geographic Routing Metric

As reported in [START_REF] Seada | Energyefficient forwarding strategies for geographic routing in lossy wireless sensor networks[END_REF], the product of PRR and the distance progress towards destination is a good decision metric. It is often more effective than a purely distance based routing metric in systems with time-varying links due to a better capture of the wireless nature. Instead of simply taking the product of PRR and the distance progress as a metric, we consider the expected proximity to the sink in terms of p(d).

The technical details are given below.

Denote the current node by i while its neighbor node by j which also represents the next-hop forwarder candidates. The routing algorithm and simulation for 1-hop based path searching is conducted in the following way:

1) The distance d i,j between pair nodes i and j is assumed known via some localization techniques.

2) The PRR, p(d i,j ), is calculated following (1), ( 2) and

(3). We assume the corresponding PRR statistics are known, i.e. the path loss exponent n, log-normal shadowing variance σ 2 , and path loss reference P L(d 0 ).

3) The progress from node i to forwarder candidate j is identified by its expected proximity to destination node, says D, defined as:

d j,D ∆ =p(d i,j ) × d j,D + (1 -p(d i,j )) × d i,D . (4)
This can also be interpreted as the expected position after 1-hop forwarding despite the fact that a packet cannot be stayed or buffered at that point. However, d j,D can serve as the metric in identifying the most suitable candidate. 4) The node with smallest d j,D is selected as the next-hop forwarder, denoted as j * . 5) The packet is then forwarded to node j * . However, this does not mean the transmission is surely successful since the link reliability can be less than 1, i.e. 0 ≤ p(d i,j * ) ≤ 1. There is a success probability of p(d i,j * ) when sending from node i and node j * in the forwarding.

In simulation, a packet forwarding will have success probability equal to p(d). When a transmission fails, a retransmission from the current node i will be initiated with new decision making. The decision will be made based on the possibly new network topology, channel variation and link updates. One may also consider it as an ARQ (automatic repeat request) mechanism but with dynamic PRR. Most geographic routing protocols use 1-hop information, but generalization to 2-hop searching or more is also possible [START_REF] Lee | Efficient geographic routing in multihop wireless networks[END_REF]. The above metric uses 1-hop information. A generalization to 2-hop based routing is established below. Here, we assume 2-hop neighborhood knowledge is obtained via some localization techniques or information exchange. Let's say, the current node is i, while its neighbor node is j. A neighbor node of j is said k. Similarly to (4), by the same concept, we define the following metric to select the best candidate among next-hop forwarders by their expected proximity to sink in 2 hops:

d (2) j,D ∆ = p(d i,j )p(d j,k )d k,D + (p(d i,j )p(d j,k ) + p(d i,j )p(d i,j ))d j,D + p(d i,j )p(d i,j )d i,D (5)
where p(•) ∆ =1-p(•) for the notational convenience. By [START_REF] Karp | GPSR: greedy perimeter stateless routing for wireless networks[END_REF], the neighbor node with smallest d

(2) j,D is selected as the next-hop forwarder.

Iteratively, a generalization to 3-hop based routing is developed below. Consider the current node is i, its neighbor node is j, and a neighbor node of j is k. The neighbor node of k is said l. Similarly to (4) and ( 5), we select the best candidate in their expected proximity to sink in 3-hop range by the following metric:

d (3) j,D ∆ = p(d i,j )p(d j,k )p(d k,l )d l,D + (p(d i,j )p(d j,k )p(d k,l ) + p(d i,j )p(d j,k )p(d j,k ) + p(d i,j )p(d i,j )p(d j,k ))d k,D + (p(d i,j )(p(d j,k )) 2 + p(d i,j )p(d i,j )p(d j,k ) + (p(d i,j )) 2 p(d i,j ))d j,D + (p(d i,j )) 3 d i,D = p(d i,j )p(d j,k )p(d k,l )d l,D + p(d i,j )p(d j,k )(p(d i,j ) + p(d j,k ) + p(d k,l ))d k,D + p(d i,j )((p(d i,j )) 2 + p(d i,j )p(d j,k ) + (p(d j,k )) 2 )d j,D + (p(d i,j )) 3 d i,D . (6)

C. Simulation Results

Sensor nodes are uniformly distributed in an area of (200m, 200m), while source and sink are located at (30m, 30m) and (170m, 170m) respectively. Following (1)- [START_REF] Doherty | Convex position estimation in wireless sensor networks[END_REF], given distance d between a pair of nodes, the corresponding PRR, p(d), is generated in each time instance. The above model and setting refers to a time-varying lossy WSN with dynamic network topology and connectivity. Data are obtained by 3000 runs.

As aforementioned, the following three environment models are employed: (i) sandy flat beach, (ii) engineering building I, and (iii) dense bamboo. Sandy flat beach is considered as common outdoor application environment, while engineering building I refers to indoor situation. Dense bamboo represents a harsh environment. Fig. 5 shows the reachability, while Fig. 6 reports the average number of transmissions required from end to end. Fig. 7 plots the standard deviation of the number of transmissions required in order to indicate the performance stability.

It is observable in Fig. 5(a) that in the lossy WSN, the new routing metric, given by (4), performs a little bit better than the simple (distance only) greedy method in the reachability. However, more importantly, from Fig. 6(a) and 7(a), the number of transmissions required by the new routing metric is significantly less. This is because the simple greedy forwarding considers the proximity to sink only by the physical distance but neglects the reliability in lossy links. Thus, it often tries to propagate in large distance so as to approach sink as much as possible. However, more likely this leads to transmission failure due to corresponding low link reliability, thus results in more retransmissions. On the other hand, the new routing metric has taken into account both the proximity and link quality, and better capture the expected progress per transmission [START_REF] Lee | Efficient geographic routing in multihop wireless networks[END_REF]. Discussions on this can also be found in [START_REF] Seada | Energyefficient forwarding strategies for geographic routing in lossy wireless sensor networks[END_REF].

Regarding the use of multi-hop neighborhood information, there is a clear increase in reachability and decrease in the average number of transmissions required and standard deviation from 1-hop searching to 2-hop searching, under all the three environments, as shown in Fig. 567respectively. The improvement is particularly significant in relatively low node density or in the dense bamboo harsh environment. Results also indicate that in the WSN with high node density and relatively good link quality, e.g. number of nodes around 200 and in engineering building I or sandy flat beach, the gain from 1-hop searching to 2-hop searching is comparatively small. However, in the dense bamboo environment, the improvement is generally substantial.

It is worth pointing out that, look at the number of transmissions required, the improvement from 2-hop searching to 3-hop searching is quite marginal. Besides, as the number of nodes increases, the improvement from a k-hop searching to (k + 1)-hop gets smaller due to the fact that the distance between nodes is now getting shorter. Thus, the PRR will be high between most neighboring nodes. Even if we just adopt the 1-hop searching in the new metric, it could be sufficient to find a good enough routing path. However, without loss of generality, the performance will always be better if more information is allowed with a proper routing strategy.

Comparing results in the three environments, they have shown similar trends. However, in harsh environment of dense bamboo, the magnitude of improvement from 1-hop searching to 2-hop searching is much larger than those in the other two. It indicates the particularly significant advantage in using multi-hop neighborhood knowledge for routing decision in harsh environment or weak links. Generally speaking, the effectiveness in sandy flat beach is still attractive. However, the improvement gets less in the engineering building model.

Note that the ∞-hop searching is conducted to find the shortest path and show the least number of transmissions required in terms of packet broadcasting in the lossy WSN. It is to provide a hint of the best performance and just serves as a benchmark for reference. The result refers to an ideal and very optimistic performance.

V. CONCLUSIONS

The study towards the asymptotic properties of geographic routing with the utilization of k-hop neighborhood information gives a general reference of how much improvement we can possibly achieve. The new metric proposed in this paper takes in account of both advance in distance and link quality. It clearly performs better than the simple greedy method. We also study the improvement of the new routing scheme when more information is utilized. It is found that there is a significant improvement in the number of transmissions required from source to sink by the multi-hop information and utilization, especially in harsh environment. It is expected that the results presented here will help to provide a reference for future design and consideration in emerging WSN and certain applications.

Fig. 1 .

 1 Fig. 1. End-to-end reachability in k-hop searching.

Fig. 3 .

 3 Fig.3. An example to show the typical scenario in which the 1-hop searching encounters more hops than the 2-hop searching does.

1 ) 2 )

 12 Sandy flat beach: n = 4.2, σ = 4.0, P L(d 0 ) = -40.8 dBm. Engineering building I: n = 1.9, σ = 5.7, P L(d 0 ) = -50.5 dBm. 3) Dense bamboo: n = 5.0, σ = 11.6, P L(d 0 ) = -38.2 dBm.

Fig. 5 .Fig. 6 .

 56 Fig.5. End-to-end reachability in the three environments.

Fig. 7 .

 7 Fig. 7. The standard deviation of the number of transmissions required.