(a)

: computation of the output SNR associated with the estimation of ϕ ′ with s f in (computed for different values of λ) or with ω [2] for the first mode of same as (a) but with input SNR= 15 dB. (d): same as (a) but using the filtered estimator F ( q f +ε ) instead of q f +ε , F ( q f ) is also superimposed ; Dictionary matrix W 0 used with FSST2-based approach. For the sake of legibility, in (c) and (d) we do not normalize the column of the matrices.

(F ( q f +ε (t, R f +ε (t))), ω [2] f +ε (t, R f +ε (t))), (F ( q f +ε (t, R f +ε (t))), ϕ ′ ),
∥.∥ ∞ L ∞ norm x ≈ y x is an approximation of y ⌊x⌋ Largest n ∈ Z such that n ≤ x ⌈x⌉ Smallest n ∈ Z such that n ≥ x ⌊x⌉ Rounding of x supp (f ) Support of the function f E{.} Expectation of a random variable ℜ(z) Real part of the complex number z ℑ(z) Imaginary part of the complex number z z * Complex conjugate of z ∈ C i Complex number that verifies i 2 = -1 δ Dirac distribution Sets N Positive integers Z Integers R Real numbers R + Positive real numbers C Complex numbers Spaces L 1 (R) Lebesgue-integrable functions L 2 (R) Squared Lebesgue-integrable functions S(R) Schwartz space S ′ (R) Tempered distribution xix List of Tables Operators Id Identity f ′ Derivative of f f (n) Derivative of order n ∂ t f Partial derivative of f with respect to t Transforms F f Fourier transform of f V g f Short-time Fourier transform of f with a window g W f Wigner-Ville distribution of f xx 1. Introduction
The mathematical frequency representation of signals was first introduced by Fourier [START_REF] Cohen | Time-frequency Analysis: Theory and Applications[END_REF], and initially used to study the propagation of heat between solids [START_REF] Baron Fourier | The analytical theory of heat[END_REF]. The Fourier transform (FT) was thus defined as a tool to perform frequency analysis, but real signals often have time varying frequencies which FT cannot capture. To analyze such signals, time-frequency representations were considered, and transforms such as the short-time Fourier transform (STFT) [START_REF] Gabor | Theory of communication. part 1: The analysis of information[END_REF], spectrogram and wavelet transform [START_REF] Grossmann | Decomposition of Hardy functions into square integrable wavelets of constant shape[END_REF] were introduced, where the latter is effective for analyzing the information content of images [START_REF] Mallat | A theory for multiresolution signal decomposition: The wavelet representation[END_REF]. Such transforms however suffer from Heisenberg's uncertainty [START_REF] Weyl | The theory of groups and quantum mechanics[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], meaning the signal cannot be arbitrarily well localized both in time and frequency. To circumvent this problem, the Wigner-Ville transform was defined [START_REF] Ville | Theorie et application dela notion de signal analytique[END_REF], but at the cost of interfering artifacts that hampered the readability of this representation [START_REF] Martin | Wigner-ville spectral analysis of nonstationary processes[END_REF]. As a consequence, other transforms came to light, as the reallocation method [START_REF] Kodera | A new method for the numerical analysis of non-stationary signals[END_REF] and the synchrosqueezing transform [START_REF] Daubechies | A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models[END_REF],

detailed later in this chapter.

Short-Time Fourier Transform and Spectrogram

The Fourier transform provides a mathematical representation of the frequencies of a signal. Given a signal f ∈ L 1 (R), it is defined by [START_REF] Hlawatsch | Time-frequency analysis[END_REF]Chapter 3],

F f (η) = R f (τ )e -2iπτ η dτ. (1.1)
The frequency representation given by F f does not depend on time. To study signals with varying frequencies, other transforms are required. In this thesis, we are going to study signals with both varying frequency and amplitude, namely non-stationary signals.

Therefore, in order to perform their time frequency (TF) analysis, we recall the definition of the most commonly used time frequency representations (TFRs).

CHAPTER 1. INTRODUCTION

Continuous Formalism

The modified STFT of a signal f ∈ L 2 (R) using a window g ∈ L 2 (R) is an invertible linear and continuous application, defined by [12, Chapter 3],

V g f (t, η) = R f (u)g(u -t)e -2iπη(u-t) du, (1.2)
and the spectrogram by

|V g f (t, η)| 2 . If f is analytic, i.e.
f only contains positive frequencies, and if f, g ∈ L 1 (R) ∩ L 2 (R), and g(0) ̸ = 0, then by applying Plancherel theorem one obtains,

f (t) = 1 g(0) R + V g f (t, η)dη. (1.3) Or if g ̸ = 0, as {u → g(t -u)e 2iπ(t-u)η } is a frame of L 2 (R), f (t) = 1 ∥g∥ 2 R 2 V g f (u, η)g(t -u)e 2iπ(t-u)η dudη. (1.4) Or if F f ∈ L 1 (R), g ∈ L 1 (R)
and R g(τ )dτ ̸ = 0, using the overlap-add method [START_REF] Portnoff | Time-frequency representation of digital signals and systems based on short-time fourier analysis[END_REF][START_REF] Crochiere | A weighted overlap-add method of short-time fourier analysis/synthesis[END_REF],

f (t) = 1 R g(τ )dτ R 2 V g f (u, η)e 2iπ(t-u)η dudη. (1.5)

Discrete Formalism

In order to proceed with the practical study of techniques expressed with continuous variables, we need to adapt the above formalism to a discrete framework. In this regard, we consider the signal f to be a discrete sequence of length L such that f

[n] = f ( n L ), for n = 0, • • • , L -1,
and that (g[n]) n∈Z are the samples at n L of the window, which is further truncated to be supported on {-M, • • • , M } such that 2M + 1 ≤ N , where N is the number of frequency bins. In that context, the discrete STFT of f is defined by [START_REF] Meignen | Retrieval of the modes of multicomponent signals from downsampled short-time Fourier transform[END_REF] V

g f [m, k] := n∈Z f [n]g[n -m]e -2iπ k(n-m) N = n=M n=-M f [m + n]g[n]e -2iπ kn N = N -1 n=0 f [n + m -M ]g[n -M ]e -2iπ k N (n-M ) , (1.6
k ∈ {0, • • • , N -1} corresponding to the frequency k L N . The discrete transform is also invertible [15], if g[0] ̸ = 0, f [n] = 1 g[0]N N 2 -1 k=-N 2 V g f [n, k]. (1.7)
Or if one assumes that f is L-periodic, g ̸ = 0,

f [n] = 1 N ∥g∥ 2 N -1 k=0 n+M m=n-M V g f [m mod L, k]g[n -m]e i2π k(n-m) N . (1.8) Or, still assuming f L-periodic, if M m=-M g[m] ̸ = 0, f [n] = 1 N M m=-M g[m] N -1 k=0 n+M m=n-M V g f [m mod L, k]e i2π k(n-m) N .
(1.9)

In these last two cases, the STFT of f is also L-periodic with respect to time.

Reassignment and Synchrosqueezing

According to the uncertainty principle, for the STFT or the spectrogram, the energy is spread over the TF plane. To circumvent this effect, considering a signal f ∈ L 2 (R) and a window g ∈ S(R), one may define a TF point ( t f (t, η), ω f (t, η)), called centroïd, on which to reallocate the energy of the STFT, provided V g f (t, η) ̸ = 0, and defined by [START_REF] Auger | Improving the readability of time-frequency and timescale representations by the reassignment method[END_REF] t f (t, η)

:= t - 1 2π ∂ η arg V g f (t, η), (1.10)
and

ω f (t, η) := 1 2π ∂ t arg V g f (t, η). (1.11)
In practice, t f (t, η) and ω f (t, η) are computed through the definition of complex reassignment operators [START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF],

tf (t, η) := t + V tg f (t, η) V g f (t, η) (1.12) ωf (t, η) := η - 1 2iπ V g ′ f (t, η) V g f (t, η) , (1.13)
and then ω f (t, η) = ℜ{ω(t, η)} and t f (t, η) = ℜ{ t(t, η)}, where ℜ{x} denotes the real part of x. When f is a pure tone, namely f = Ae 2iπϕ(t) with ϕ(t) = ωt, then ω f (t, η) = ω for any (t, η) in the TF plane [START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF]. This means that ω f gives the frequency ω of the pure tone f . In the general case, where the amplitude A(t) and the frequency ϕ ′ (t) depends on time, ω f becomes an estimate of ϕ ′ . It is also important to note that the operators ωf and tf only depend on the computation of three STFTs and are therefore easy to compute. Then, assuming f ∈ L 1 (R) ∩ L 2 (R) and R g ̸ = 0, the reassignment method (RM) associated with t f and ω f is defined by [START_REF] Auger | Improving the readability of time-frequency and timescale representations by the reassignment method[END_REF],

RM (t, η) = R 2 |V g f (τ, ξ)| 2 δ(t -t f (τ, ξ))δ(η -ω f (τ, ξ))dτ dξ (1.14)
However, RM is not invertible because the phase is lost in the process, and therefore one may instead consider reallocating the STFT along the frequency axis only, which brings about the Fourier-based synchrosqueezing transform (FSST) whose principle is detailed hereafter.

Definition of the synchrosqueezing

FSST was introduced in [START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF][START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], and using both assumptions (1.5) and g ∈ S(R), it is defined by,

T g f (t, η) := R V g f (t, ξ)δ(η -ω f (t, ξ))dξ. (1.15)
Provided g(0) ̸ = 0 one has, when f is analytic, .16) This transform leads to a sparser TFR than the STFT but is very sensitive to frequency modulation. In that context, to improve ω f , the following complex modulation operator was introduced [START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF],

f (t) = 1 g(0) R + T g f (t, η)dη. ( 1
qf (t, η) := ∂ t ωf (t, η) ∂ t tf (t, η) = 1 2iπ V g ′′ f (t, η)V g f (t, η) -V g ′ f (t, η) 2 V tg f (t, η)V g ′ f (t, η) -V tg ′ f (t, η)V g f (t, η)
.

(1.17

)
Local modulation is given by q f (t, η) := ℜ{q f }. An illustration of q f is given in Figure 1.1. Furthermore, for a linear chirp with Gaussian modulated amplitude f , namely t) where log(A(t)) and ϕ(t) are second order polynomials, q f (t, η) is the chirp rate (CR), namely q f = ϕ ′′ (t) [START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF][START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF]. Using qf , the second order frequency reallocation operator was introduced,

f (t) = A(t)e 2iπϕ(
ω [2] f (t, η) = ℜ ωf (t, η) -qf (t, η)(t -tf (t, η)) , (1.18)
which verifies ω [START_REF] Baron Fourier | The analytical theory of heat[END_REF] f (t, η) = ϕ ′ (t). The second order Fourier-based synchrosqueezing transform (FSST2) can then be defined [START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF][START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF],

T g,2 f (t, η) := R V g f (t, ξ)δ(η -ω [2] f (t, ξ))dξ. (1.19)
Note that to evaluate qf and ω [START_REF] Baron Fourier | The analytical theory of heat[END_REF] f only requires to compute five STFTs. The Fourier based synchrosqueezing has been extended to the N -th order [START_REF] Meignen | Synchrosqueezing transforms: From low-to high-frequency modulations and perspectives[END_REF], denoted T g,N f , which we will detail and use only in Chapter 4. Note that the synchrosqueezing transform (SST) was initially introduced in the wavelet framework [START_REF] Daubechies | Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool[END_REF], but in the context of this thesis, we will stick to the STFT formalism.

Multicomponent signals

Many real world non-stationary signals can be described using the formalism of multicomponent signal s (MCSs). For instance, real world audio signals such as piano, clarinet [START_REF] Gribonval | Harmonic decomposition of audio signals with matching pursuit[END_REF], bird songs [START_REF] Sandsten | Classification of bird song syllables using wignerville ambiguity function cross-terms[END_REF], speech [START_REF] Choi | Improved time-frequency representation of multicomponent signals using exponential kernels[END_REF] but also signals from electrocardiograms [START_REF] Herry | Heart beat classification from single-lead ECG using the synchrosqueezing transform[END_REF][START_REF] Acharya | Heart rate variability[END_REF][START_REF] Malik | Dynamic electrocardiography[END_REF], thoracic and abdominal movements [START_REF] Lin | Sleep apnea detection based on thoracic and abdominal movement signals of wearable piezoelectric bands[END_REF], and underwater propagation [START_REF] Bonnel | Physics-based time-frequency representations for underwater acoustics: Power class utilization with waveguideinvariant approximation[END_REF].

A MCS f containing P modes is defined by, Therefore, MCS are commonly analyzed in time and frequency [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF].

f (t) = P p=1 f p (t), (1.20
Then, its ideal TFR is given by,

IT F f (t, η) = P p=1 A p (t)δ(η -ϕ ′ p (t)), (1.21)
where δ is the Dirac distribution. It is said ideal in the sense that it is perfectly localized, both in time and frequency. The SST aims towards this representation and is close to it when the signal meets the assumptions the transform is designed for. An illustration of ideal time frequency representation (ITF) is given in Figure 1.3.

Mode retrieval

When considering a MCS f , there is no general formula to recover each of its modes.

The process of computing an approximation of a mode is called mode retrieval (MR).

Existing techniques use various TFRs such as the continuous wavelet transform, as for instance in [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF], or in the context of synchrosqueezing [START_REF] Daubechies | Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool[END_REF]. When STFT is used, one can take advantage of its linearity, namely,

V g f (t, η) = P p=1 V g fp (t, η), (1.22)
as for instance, in [START_REF] Meignen | Retrieval of the modes of multicomponent signals from downsampled short-time Fourier transform[END_REF] to perform MR based on hard-thresholding.

Then, the retrieval of f p is performed using one of the approximation formulae of

V g fp (t, η) (1.
3), (1.4) or (1.5), depending on what assumptions are made on f and g. FSST was also used for that matter, for instantaneous frequency estimation [START_REF] Thakur | Synchrosqueezing-based recovery of instantaneous frequency from nonuniform samples[END_REF], for non-stationnary modes [START_REF] Thakur | The synchrosqueezing algorithm for time-varying spectral analysis: Robustness properties and new paleoclimate applications[END_REF], a thorough mathematical study being available in [START_REF] Oberlin | The Fourier-based synchrosqueezing transform[END_REF].

When the separation in the TF plane is difficult, adaptive approaches can also be considered, namely the analysis window may vary in time or frequency, as in the case of adaptive wavelets based SST [START_REF] Li | Adaptive synchrosqueezing transform with a timevarying parameter for non-stationary signal separation[END_REF], or in the case of adaptive STFT based SST [START_REF] Li | Adaptive short-time Fourier transform and synchrosqueezing transform for non-stationary signal separation[END_REF] which has recently been used, in [START_REF] Chui | Analysis of an adaptive short-time fourier transform-based multicomponent signal separation method derived from linear chirp local approximation[END_REF] for multicomponent signal separation.

For the sake of completion, we note that other techniques, not based on TFR, are also used, as those based on chirplet [START_REF] Mann | The chirplet transform: Physical considerations[END_REF] or ridglet [START_REF] Candes | Ridglets: The key to high-dimensional intermittency[END_REF] transforms.

Outline of the Following Chapters

In Chapter 2, we will carry out TF analysis of a MCS affected by noise in order to perform MR. In our approach, we are first going to analyze the TF plane to associate TF regions to each mode. This task is called ridge detection (RD) as those regions often look like ridges due to the amount of energy they contain [START_REF] Carmona | Multiridge detection and timefrequency reconstruction[END_REF]. Then, we perform a local analysis in the vicinity of the ridges in order to approximate the corresponding modes using a linear chirp based model. This enables us to express the STFT of the mode in a closed form and then to define a numerical procedure to perform MR.

In Chapter 3, we will perform a TF analysis of the interference of MCSs. We will first characterize an interference pattern called time frequency bubble (TFB) theoretically and subsequently introduce a technique for their identification in the TF plane, which we validate numerically. At end of the chapter, we will present perspective using a different TF approach to identify interference based on the spectrogram zeros. For this purpose, we study the energy distribution related to the Voronoi cell, which will be continued in a near future.

In Chapter 4, we will recall the definition of high order Fourier-based synchrosqueezing transform (FSSTN), and then, we will show how to make an adaptive choice of its order based on an energy based criterion, and validate this technique numerically.

In Chapter 5, considering that FSST2 is relevant in very noisy situations, we will study its associated CR estimator. We illustrate that the errors introduced by the noise have an oscillatory nature and propose a technique to filter out most of them. We will finally validate its accuracy on numerical experiments and its stability on real data.

In Chapter 6, we will study a TF approach for heart rate (HR) estimation on electrocardiogram (ECG) signals. Assuming that an ECG signal can be represented by a harmonic model, we will propose a TF technique for HR estimation, and show that the choice of the TFR has a huge impact on the accuracy. Finally, we will introduce an algorithm based on FSST for the estimation of the HR and evaluate it on synthetic and real ECG signals. [START_REF] Stanković | A measure of some time-frequency distributions concentration[END_REF][START_REF] Stankovic | Performance of spectrogram as IF estimator[END_REF]. The detection of those curves is preformed by a RD.

In a first part of this chapter, we consider the ridge detection in the case of noise. In the literature, there exist RD based on different TFR, for instance the continuous wavelet transform [START_REF] Carmona | Characterization of signals by the ridges of their wavelet transforms[END_REF], the Wigner Ville transform [START_REF] Djurović | An algorithm for the Wigner distribution based instantaneous frequency estimation in a high noise environment[END_REF], or on STFT [START_REF] Colominas | Fully adaptive ridge detection based on STFT phase information[END_REF]. In our approach, we use TF regions with high energy as an indicator of signal information [START_REF] Carmona | Multiridge detection and time-frequency reconstruction[END_REF][START_REF] Zhu | Two robust approaches to multicomponent signal reconstruction from STFT ridges[END_REF]. The technique we detail is designed to be robust to noise, and therefore, it does not assume that the ridge is continuous. Instead of trying to follow a ridge, it uses structures based on local maxima along the frequency axis (LMMFs).

In the second part of this chapter, we design a MR technique based on a RD for noisy signals. We show that even with perfect ridge detection, the quality of the reconstruction cannot be guaranteed. One possible approach is to perform hard thresholding and proceed to the reconstruction of the signal [START_REF] Meignen | Retrieval of the modes of multicomponent signals from downsampled short-time Fourier transform[END_REF]. However, TF regions related to the signal are still affected by noise. Our approach tries to limit the effect of the noise in the signal domain by using a linear chirp model.

At the end of this chapter, we finally perform numerical validation of the techniques and draw conclusions on their behaviors.

CHAPTER 2. MODE RETRIEVAL BASED ON ROBUST RIDGE DETECTION AND LINEAR CHIRP APPROXIMATION

As we will always consider noisy situations, let us introduce the noisy signal,

f (t) := f (t) + ε(t), (2.1)
where f is the MCS and ε is a complex, zero mean, white Gaussian noise with variance σ 2 ε . We assume ε to be independently and identically distributed (i.i.d) with respect to time and that its real and imaginary parts are independent as well. In the remainder of this chapter, we will also use a discrete setting for f , by considering it as a sequence of length L with f [n] = f n L .

Basics of Ridge Detection

Performing ridge detection on f consists of associating TF curves to the modes of f . The objective of this section is to understand basics of what makes a RD robust to noise.

Therefore, we first study the effect of ε on non-robust RDs, the focus being put on a very commonly used one, which we call simple ridge detector (S-RD), in the sequel and that we introduce hereafter.

Simple STFT-Based Ridge Detection

The RD introduced by [START_REF] Carmona | Characterization of signals by the ridges of their wavelet transforms[END_REF] is based on finding the location of the maxima of energy in the spectrogram. It consists of determining the P ridges associated with the modes in the TF plane by computing:

max Γ P p=1 L-1 n=0 |V g f [n, Γ p [n]] | 2 -α(∆ 1 Γ p [n]) 2 -β(∆ 2 Γ p [n]) 2 , (2.2)
with Γ = (Γ p ) p=1,••• ,P , α and β both positive, and in which

∆ 1 Γ p [n] (resp. ∆ 2 Γ p [n]
) is the first order (resp. second) order finite differences approximating ϕ ′′ p [n] (resp. ϕ

p [n]). To consider penalization terms is however not relevant at low noise level [START_REF] Meignen | On demodulation, ridge detection, and synchrosqueezing for multicomponent signals[END_REF], and these are therefore often left apart. Alternatively, one can use a bound B f on the frequency modulation of the modes to extract a first ridge, and then replace (2.2) by a peeling algorithm. In a nutshell, a first ridge is extracted as follows [START_REF] Colominas | Fully adaptive ridge detection based on STFT phase information[END_REF]:

max Γ 1 L-1 n=0 |V g f [n, Γ 1 [n]] | 2 , s.t. |∆ 1 Γ 1 [n]| ≤ B f . (2.3) 2.2. BASICS OF RIDGE DETECTION In practice, to compute Γ 1 , one initializes Γ 1 [n 0 ] to k 0 := argmax 0≤k≤N -1 |V g f [n 0 , k]|.
Then, to define Γ 1 on {n 0 + 1, • • • , L -1}, one uses the following induction principle starting from n = n 0 :

Γ 1 [n + 1] := argmax k |V g f [n + 1, k]|, s.t. |k -Γ 1 [n]| ≤ B f . (2.4)
The same principle is applied on {0, • • • , n 0 -1}, starting from n = n 0 and replacing n + 1 by n -1 in (2.4). Finally, the procedure is run again starting from other initial time indices to define other candidates for Γ 1 , and the ridge finally kept among all the candidates is the one maximizing the energy in the TF plane, i.e.

n |V g f [n, Γ 1 [n]]| 2 .
Then, to compute the next ridges, one introduces the following set, for p > 1:

S p (n) = {0 ≤ k ≤ N -1, s.t |k -Γ q [n]| ≥ ∆, 1 ≤ q < p} , (2.5)
and then the p th ridge is computed using the same principle as for Γ 1 except the initialization index is defined as

k 0 := argmax k∈Sp(n) |V g f [n 0 , k]| and then (2.4
) is replaced by:

Γ p [n + 1] := argmax k∈Sp(n+1) |V g f [n + 1, k]|, s.t.|k -Γ p [n]| ≤ B f . (2.6)
Such a technique is iterated until P ridges (Γ p ) p=1,••• ,P are extracted. There are however two strong limitations to S-RD:

* Each ridge is built assuming that the selection of the maximum around the current position provides with a continuous curve on which the coefficients are associated with signal information. This assumption does not hold when the noise is too strong. Indeed, in such a case, zeros of STFT may appear at TF locations corresponding the IFs of the modes. This is illustrated on a linear chirp in Figure 2.1 in which we display the magnitude of the LMMFs associated with the three largest STFT modulus maxima at each time instant, along with the true IF (signal to noise ratio (SNR) = -10 dB).

* The jumps allowed between two successive time indices depend on the modulation parameter B f which is fixed a priori, and thus this method does not adapt to the local variations of the frequency modulation of the modes.

We next recall the RD proposed in [START_REF] Colominas | Fully adaptive ridge detection based on STFT phase information[END_REF], where the chaining of the coefficient is made more adaptive by removing the modulation parameter B f and exploiting the modulation

Modulation Based Ridge Detection

To circumvent the lack of adaptivity of S-RD to the local variations of the frequency modulation of the modes, a novel approach called modulation based ridge detector (MB-RD) was proposed in [START_REF] Colominas | Fully adaptive ridge detection based on STFT phase information[END_REF]. In a nutshell, this approach extracts the first ridge, MB-RD uses the same induction principle as S-RD, but considers the complex modulation operator q f defined in Section 1.2 instead of B f , namely,

Γ 1 [n + 1] := argmax k |V g f [n + 1, k]|, |k -Γ 1 [n] -q f [n, Γ 1 [n]] N L 2 | ≤ C , (2.7)
the user-defined constant C compensating for potential local frequency modulation estimation errors. MB-RD is proved to be slightly sensitive to C when the noise level is low [START_REF] Colominas | Fully adaptive ridge detection based on STFT phase information[END_REF] (the simulations in that paper only considered a SNR larger than 0 dB), but the estimator q f is not robust enough when the noise is too strong, which leads to bad results.

Furthermore, both techniques are based on the assumption that the IF of a mode at each time index can be associated with a continuous sequence of energetic TF coefficients, which may not be the case in heavy noise situations. Another limitation of S-RD and MB-RD is that they build the ridges one after the other using the peeling algorithm recalled in Section 2.2.1: if the ridge detection fails for one mode, it will also fail for the 2.3. DEFINITION OF A NEW ROBUST RIDGE DETECTOR next ones. To deal with all these issues, the concept of relevant ridge portions (RRPs) is introduced in the following section, and subsequently used to define a new RD, for which the approach significantly differs from S-RD and MB-RD.

Definition of a New Robust Ridge Detector

The objective of this section is to define a RD robust to noise. Its algorithm can be divided in three major steps. The first step consists of identifying energetic chains of coefficients belonging to the signal, called RRPs. The second step associates RRPs that have connected basins of attraction. Such a basin is a TF regions, where all the coefficients points towards the same RRP. The last step consists of two nested iterative procedures. The first one selects a combination of RRPs to associate with each mode.

Each of them represents initial TF coefficients that can be trusted to be associated with one mode of the signal. The second one, uses the initialization given by the first procedure to associate a TF curve with each mode using an iterative procedure based on weighted spline. We name the solution given by this algorithm relevant ridge portion based ridge detector (RRP-RD).

Definition of Relevant Ridge Portions

To have a robust RD, we first define a structure that persists even when noise strongly impacts the signal. For that purpose, we use q f to define a stable direction in the TF plane between two neighboring LMMFs. Denoting a generic LMMF by [n, m[n]], we introduce,

Definition 2.3.1. Let [n, m[n]] and [n + 1, m[n + 1]
] be two LMMFs, they satisfy relation

([n, m[n]] ∼ [n + 1, m[n + 1]]) if and only if, m[n + 1] = argmin k k -m[n] -q f [n, m[n]] N L 2 , s.t. [n + 1, k] LMMF m[n] = argmin k m[n + 1] -k -q f [n + 1, m[n + 1]] N L 2 , s.t. [n, k] LMMF . (2.8) If [n, m[n]] ∼ [n + 1, m[n + 1]], q f is a relevant operator giving the direction from [n, m[n]] to [n + 1, m[n + 1]
] and vice-versa, meaning that q f gives a stable orientation.

A ridge portion (RP) is then defined by extending relation ∼ to any time instants n, n 0 that are not necessarily neighbors. Assuming, without loss of generality, that n > n 0 , we define relation ↔,

Definition 2.3.2. Let [n, m[n]] and [n 0 , m[n 0 ]] be two LMMFs such that n ≥ n 0 + 1,
then define:

([n 0 , m[n 0 ]] ↔ [n, m[n]]) ⇔    ∃ [n 0 + 1, m[n 0 + 1]], • • • , [n -1, m[n -1]] LMMFs ∀i = n 0 , • • • , n -1, [i, m[i]] ∼ [i + 1, m[i + 1]] (2.9) A RP R containing LMMF [n 0 , m[n 0 ]],
is finally defined by:

R[n 0 , m[n 0 ]] = {[n, m[n]], s.t. [n, m[n]] ↔ [n 0 , m[n 0 ]]} . (2.10) Note that, if [n, m[n]] ∈ R[n 0 , m[n 0 ]] then R[n, m[n]] = R[n 0 , m[n 0 ]].
The definition of RPs allows considering the discontinuities in the ridges associated with a mode. However, to separate the noise from the signal, a criterion needs to be defined. Using the assumptions made on the noise ε [53], we have that,

Var (ℜ{V g ε [n, k]}) = Var (ℑ{V g ε [n, k]}) = σ 2 ε ∥g∥ 2 2 .
Then, remarking that

|V g ε | 2 σ 2 ε ∥g∥ 2 2
is χ 2 distributed with two degrees of freedom and assuming the variance of the noise σ 2 ε is known, the probability that

|V g ε [n, k]| ≥ βσ ε ∥g∥ 2 is smaller than 10% if β = 2, and 1% if β = 3. So, by considering only the LMMFs [n, m[n]] such that |V g f [n, m[n]]| ≥ βσ ε ∥g∥ 2
with β ≥ 2 one removes many LMMFs corresponding to noise. In practice, σ ε is unknown, but can be approximated using the estimate γ = σ ε ∥g∥ 2 defined by [START_REF] Donoho | Ideal spatial adaptation via wavelet shrinkage[END_REF][START_REF] Sachs | Applied Statistics: A Handbook of Techniques[END_REF],

γ = median ℜ V g f [n, k] n,k 0.6745 ,
in which median represents the median of the coefficients. One can then set a meaningful threshold on the STFT , taking into account the above probabilities, and considering

S(β) = [n, k], |V g f [n, k]| ≥ βγ . ( 2
[n 0 , m[n 0 ]] ↭ [n, m[n]] ⇔          ∃ [n 0 + 1, m[n 0 + 1]], • • • , [n -1, m[n -1]] LMMFs ∀i = n 0 , • • • , n -1,    [i, m[i]] ∼ [i + 1, m[i + 1]] [i, m[i]] ∈ S(β) (2.12)
Any LMMF belonging to a RRP is called a relevant LMMF. To properly separate signal from noise, a study of the parameter β is required.

Fixing the parameter β for the RRPs

A crucial issue is then how to fix the parameter β. For that purpose, we numerically investigate the existence of relevant LMMFs around the true IF when β varies. To define a meaningful interval around the true IF, we first compute analytical expressions of the STFT assuming a local linear chirp approximation for the modes. To this end, we set g to be a Gaussian window for two reasons:

* First because it simplifies the computation of q f . * Second because the expression of the STFT of a signal that can be locally approximated by a linear chirp is particularly simple [START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF].

Indeed, let g be the discrete Gaussian window g[n] = e -π n 2 σ 2 L 2 , to increase the readability of the TFR, we choose the value of σ whose associated STFT modulus minimizes the Rényi entropy and should therefore provide good localization. Such a choice for σ is proved to minimize interference between the modes in the TF plane [START_REF] Baraniuk | Measuring timefrequency 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]. With such a window, if f can be locally approximated by a linear chirp with constant amplitude A, i.e. f [n] ≈ Ae 2iπϕ [n] , with ϕ a second order polynomial, one has [START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF],

|V g f [n, k]| ≈ ALσ(1 + σ 4 ϕ ′′ [n] 2 ) -1 4 e -π σ 2 (k L N -ϕ ′ [n]) 2 1+σ 4 ϕ ′′ [n] 2 , (2.13)
which is a Gaussian function with standard deviation,

∆ LC [n] = 1 √ 2πσ 1 + σ 4 ϕ ′′ [n] 2 . (2.14)
For V g f , the LMMF associated with the largest STFT modulus maximum at time index n has its ordinate in the interval

I[n] = [(ϕ ′ [n] -∆ LC [n]) N L , (ϕ ′ [n] + ∆ LC [n]) N L ]
. When some noise is added, this is no longer true. In the following experiment, we assume ϕ ′ CHAPTER 2. MODE RETRIEVAL BASED ON ROBUST RIDGE DETECTION AND LINEAR CHIRP APPROXIMATION and ϕ ′′ are known, and assess the proportion corresponding to the number of relevant LMMFs with ordinate in I[n] when n varies, namely:

P(β) = # {[n, m[n]] ∈ S(β) with m[n] ∈ I[n]} L , (2.15)
where #X denotes the cardinal of the set X. Using numerical validation, we show how to set β as high as possible to filter out the noise while keeping P(β) as high as possible.

Such a value for β enables to keep a maximum of relevant LMMFs in I which will help us associate them to the same modes by connecting their corresponding RRPs.

We carry out the validation of the choice of β for the three signals of

Grouping RRPs Based on Basins of Attraction

In this subsection, the objective is to group RRPs using the connection of their basin of attractions (BAs) to finally assign these groups to one mode of the signal. For that purpose, we will use the reallocation vectors

( t f [n, k], ω f [n, k]) defined in 1.2. Indeed, a
BA associated with the RRP R i can be defined by,

B i := [n, k]; argmin [x,y]∈RRPs ∥( t f [n, k], ω f [n, k]) -[x, y]∥ ∈ R i . (2.16)
However, many BAs are associated with noise, and therefore we impose more conditions to keep the most energetic part of B i . To do so, we consider a hard thresholding (HT) approach,

B HT i = B i S(2) if R i S(3) ̸ = ∅ ∅ otherwise , (2.17)
In short, B HT i is more likely related to the signal than B i , since it contains at least one coefficient in S(3).

We then gather together connected B HT i s in the TF plane to obtain a set of larger TF regions which are denoted by {(C HT j ) j } in the sequel.

New RD Definition

In this subsection, we detail how we use the gathered BAs to associate a TF a curve to each mode of the signal. This is done using Algorithm 1 which we detail here. First, in order to avoid mixing ridges associated with different modes, the algorithm looks for P elements in {(C HT j ) j } coexisting in time. Indeed, ridges that coexist in time cannot belong to the same mode. In this approach, we consider the P -tuple coexisting on the longest set of time indices first. Therefore, Algorithm 1 starts with the ranking procedure "sort-bycoexistence", to rank P -tuples of {(C HT j ) j } from the ones coexisting the longest to the shortest, and we denote this ranking by

({C κ p=1,••• ,P }) κ . The first tuple of this ranking, namely (C 0 p ) p=1,••• ,P , should contain one long ridge for each mode p. Because {(C HT j ) j } is finite, ({C κ p=1,••• ,P })
κ is a finite sequence, and we denote the last element index by κ max . Next, we iterate over this ranked set to associate possible solutions for the ridge detection. At the κ-th iteration, where κ ∈ {0, • • • , κ max }, we consider that the coefficients defined by 3) represent robust signal information that should be used as a basis to define the ridge of the mode p. In Algorithm 1, the procedure "fit-spline" defines an initial ridge using the following weighted spline approximation formula,

A κ p = (C 0 p ∪C 1 p ∪• • •∪C κ p )∩S(
s κ p,initial = argmin s   (1 -λ) [n,m[n]]∈A κ p m[n] L N -s( n L ) 2 |V g f [n, m[n]]| + λ 1 0 (s ′′ (t)) 2 dt,   .
(2.18)

Then, "fit-spline" further iterates, by taking into account the elements in {(C HT j )} that s κ p,initial crosses in (2.18). It continues to iterate until the spline stabilizes, meaning no new element are crossed in {(C HT j )}. We denote the stabilized spline by s κ p .

Algorithm 1: New RD Input :

{C HT j },P , S(3) (C κ p ) = sort-by-coexistence({C HT j }) for κ = 0, • • • , κ max do for p = 1, • • • , P do A κ p = C 0 p ∪ • • • ∪ C κ p ∩ S(3) s κ p = fit-spline(A κ p , {C HT j }) Compute energy E κ p using (2.20) end end Compute index κ f in using (2.19) Define the final ridges as s f in p := s κ f in p
Finally, we define the index κ f in of the solution by

κ f in := arg max κ s.t. (s κ p )p are not crossing. P p=1 E κ p , (2.19)
where E κ p corresponds to the energy associated with mode p and is defined by,

E k p := [n,m[n]]∈A k p , [n,⌊s k p [n] N L ⌉]∈C k p |V g f [n, m[n]]|, (2.20)
in which ⌊X⌉ denotes the nearest integer to X. In short, {s f in p } is chosen among the κ max + 1 solutions. It maximizes an energy and contain no crossing splines. We denote this ridge detector RRP-RD.

We display on the first row of

Comparison of RRP-RD, S-RD and MB-RD on Simulated Signals

In this subsection, we draw conclusions on the efficiency of RRP-RD compared with S-RD and MB-RD. In Figure 2.4 we show numerical results on the comparison of these AND LINEAR CHIRP APPROXIMATION algorithms when the input SNR varies between -10 and 0 dB for the three signals of

where y is an approximation of x. The first thing to notice is that RRP-RD always gives a significantly better SNR, from which we deduce that RRP-RD is indeed more robust to noise. Furthermore, the result are better when the smoothing parameter λ is larger.

This can be explained by the fact that for those three signals, the IF is smooth, and a large λ enables to reduce the impact of the irregularities brought by strong noise.

We notice that the results for the two linear chirps of It is worth noting that, because the noise is strong, S-RD behaves better than MB-RD.

Indeed, this can be explained by the fact that MB-RD relies on the unstable frequency modulation estimate q f . In fact, the latter may point very far from the true IF of a mode, even if it is evaluated close to the IF, thus explaining why S-RD leads to better results.

In the sequel, we will study the improvement brought by RRP-RD for MR.

LINEAR CHIRP APPROXIMATION

Linear Chirp Approximation

Once the ridges are detected, another issue is that of MR. For that purpose, we perform local TF analysis around the ridge to get an approximation of its associated mode, based on a simplified linear chirp model. Let us assume that the signal is a linear chirp of the form f (t) = Ae 2iπ(at+b t 2 2 ) , and g(t) is the Gaussian function e -π t 2 σ 2 . The STFT can then be seen as the Fourier transform of a Gaussian, and reads,

V g f (t, η) = f (t) R e -π 1 σ 2 +iϕ ′′ (τ -t) 2 e -2iπ(τ -t)(η-ϕ ′ (t)) dτ. (2.22) CHAPTER 2.

MODE RETRIEVAL BASED ON ROBUST RIDGE DETECTION AND LINEAR CHIRP APPROXIMATION

For a Gaussian function of the form h(t) = e -πzt 2 , where z = re iθ , with r ∈ R + and θ ∈ R, such that ℜ(z) > 0, its Fourier transform admits the following closed form,

F(h)(η) = r -1 2 e -i θ 2 e -π z η 2 . (2.23)
The proof is analogous to the case where z is real. Based on (2.23), one can deduce the analytic expression of V g f [START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF],

V g f (t, η) = Ar -1 2 e -i θ 2 e -πσ 2 (1+ibσ 2 ) 1+(bσ 2 ) 2 (η-a-bt) 2 = V g f (t, a + bt)e -πσ 2 (1+ibσ 2 ) 1+(bσ 2 ) 2 (η-a-bt) 2 = V g f (t, ϕ ′ (t))e -πσ 2 (1+iϕ ′′ (t)σ 2 ) 1+(ϕ ′′ (t)σ 2 ) 2 (η-ϕ ′ (t)) 2 , (2.24)
with r =

1 σ 4 + ϕ ′′ (t) 2 and θ = tan -1 (-ϕ ′′ (t)σ 2 ).
If the monocomponent signal f departs from a linear chirp, the above equalities become approximations, which are furthermore only valid for η in the vicinity of ϕ ′ (t). Applying this approximation to f p , we obtain for η in the vicinity of ϕ ′ p (t),

V g fp (t, η) ≈ V g fp (t, ϕ ′ p (t))e -πσ 2 (1+iϕ ′′ p (t)σ 2 ) 1+(ϕ ′′ p (t)σ 2 ) 2 (η-ϕ ′ p (t)) 2 (2.25)
Then, remarking that at TF location (t, ϕ ′ p (t)), the STFT of the mode f q , for q ̸ = p, is close to 0 we replace in (2.25) f (t, Γ p (t)), and ψ ′′ p (t) = q f (t, Γ p (t)), this enables, from (2.24), to approximate the STFT of f p by

V g fp (t, ϕ ′ p (t)) by V g f (t, ϕ ′ p (t)), V g fp (t, η) ≈ V g f (t, ϕ ′ p (t))e -πσ 2 (1+iϕ ′′ p (t)σ 2 ) 1+(ϕ ′′ p (t)σ 2 ) 2 (η-ϕ ′ p (t)) 2 . ( 2
V g fp (t, η) ≈ V g f t, ψ ′ p (t) e - πσ 2 (1+iψ ′′ p (t)σ 2 ) 1+(ψ ′′ p (t)σ 2 ) 2 ( η-ψ ′ p (t) ) 2 .
(2.27)

Finally, the reconstruction of mode f p is carried out either through equation (1.5, page

2), (1.4, page 2), or (1.3, page 2) replacing V g f by the approximation of V g fp given by (2.27).

It is important to remark here that the approximation of V g fp uses only STFTs of f evaluated in the vicinity of the ridge associated with f p , which will prove to be the very reason why the just described approach for mode retrieval is well adapted to handle the mode-mixing issue. However, up to now, we have only presented the approach for continuous time signals, and it needs to be adapted to discrete time signals, as explained hereafter.

Linear Chirp Based Mode Retrieval Technique

To implement the just described mode retrieval technique, we use the discrete setting of f , recalling that f [n] = f ( n L ). In that context, the discrete mode

(f p [n]) n is associated with a ridge (Γ p [n]) n corresponding to frequency indices in {0, • • • , N -1}. We then set ψ ′ p [n] := ω [2] f [n, Γ p [n]] and ψ ′′ p [n] := q f [n, Γ p [n]] approximating ϕ ′ p ( n L ) and ϕ ′′ p ( n L )
respectively. Remarking, using rectangular integration, that

V g f ( n L , k L N ) ≈ V g f [n,k] L
, using (2.26) and (2.27), we obtain,

V g fp [n, k] ≈ LV g fp ( n L , k L N ) ≈ LV g f n L , ϕ ′ p ( n L ) e - πσ 2 (1+iϕ ′′ p ( n L )σ 2 ) 1+(ϕ ′′ p ( n L )σ 2 ) 2 ( kL N -ϕ ′ p ( n L ) ) 2 ≈ LV g f n L , ψ ′ p [n] e - πσ 2 (1+iψ ′′ p [n]σ 2 ) 1+(ψ ′′ [n]σ 2 ) 2 ( kL N -ψ ′ p [n]) ) 2 , (2.28)
for which an approximation of

V g f n L , ψ ′ p [n]
x. Because V g f n L , ψ ′ p [n] ≈ V g fp n L , ψ ′ p [n]
, using (2.27), we may write,

V g f n L , ψ ′ p [n] ≈ V g f ( n L , k 0 L N )e πσ 2 (1+iψ ′′ p [n]σ 2 ) 1+(ψ ′′ p [n]σ 2 ) 2 (k 0 L N -ψ ′ p [n]) 2 ≈ 1 L V g f [n, k 0 ]e πσ 2 (1+iψ ′′ p [n]σ 2 ) 1+(ψ ′′ p [n]σ 2 ) 2 (k 0 L N -ψ ′ p [n]) 2 .
(2.29)

From (2.28) and (2.29) we get that,

V g fp [n, k] ≈V g f [n, k 0 ]e πσ 2 (1+iψ ′′ p [n]σ 2 ) 1+(ψ ′′ p [n]σ 2 ) 2 [ ( k 0 L N -ψ ′ p [n]) 2 -( kL N -ψ ′ p [n]) 2 ] ≈V g f [n, k 0 ]e πσ 2 (1+iψ ′′ p [n]σ 2 ) 1+(ψ ′′ p [n]σ 2 ) 2 [ L(k 0 -k) N ( L(k 0 +k) N -2ψ ′ p [n]) ]
.

(2.30)

If we denote Ṽ g fp the estimation of V g fp given by (2.30), the retrieval of f p can be carried out either through (1.7, page 3), if g[0] ̸ = 0,

f p [n] ≈ 1 g[0]N N 2 -1 k=-N 2 Ṽ g fp [n, k], (2.31)
or assuming f is L-periodic, through (1.8, page 3), if we assume f is L-periodic, g ̸ = 0,

f p [n] ≈ 1 N ∥g∥ 2 N -1 k=0 n+M m=n-M Ṽ g fp [n, k]g[n -m]e i2π k(n-m) N . (2.32)
In a noisy context, all the above expressions can be used replacing f by its noisy version f . The techniques given by (2.31) and (2.32) are respectively denoted M 1 and M 2 . When making abstraction of the reconstruction given by (1.7, page 3) and (1.8, page 3), this technique is called linear chirp based retrieval (LCR).

Alternative Techniques for Mode Retrieval: HT, SSR-HT

The HT technique when used for the retrieval of the mode f p of a MCS [START_REF] Meignen | Retrieval of the modes of multicomponent signals from downsampled short-time Fourier transform[END_REF] consists of considering in the reconstruction process only some coefficients of V g f in the vicinity of Γp , the ridge associated with fp the noisy version of f p . For that purpose, one defines 2.7. ALTERNATIVE TECHNIQUES FOR MODE RETRIEVAL: HT, SSR-HT for each time indexed by n, an interval

J p [n] = η - p [n], η + p [n] such that: η - p [n] := argmax k kL N < Γp [n], |V g f [n, k]| < 3γ η + p [n] := argmin k kL N > Γp [n], |V g f [n, k]| < 3γ .
(2.33)

The reconstruction of the modes is carried out summing the coefficients over J p using

(1.7, page 3), if g[0] ̸ = 0, f p [n] ≈ 1 g[0]N k∈Jp[n] V g f [n, k], (2.34)
or assuming f is L-periodic, through (1.8, page 3), if we assume f is L-periodic, g ̸ = 0,

f p [n] ≈ 1 N ∥g∥ 2 k∈Jp[n] n+M m=n-M V g f [n, k]g[n -m]e i2π k(n-m) N . (2.35)
Similarly to LCR, we respectively name those techniques HT M 1 and HT M 2 .

For the sake of comparison, we also briefly recall the principles of a denoising technique called shifted-symmetrized-regularized hard-thresholding (SSR-HT), and introduced in [START_REF] Pham | A novel thresholding technique for the denoising of multicomponent signals[END_REF], which consists of an improvement of HT technique based on a linear chirp approximation for the mode. It first shifts the real and imaginary parts of V g f such that their respective maximum in J p is centered on Γ p in the frequency axis. We denote the shifted STFT by V S g p, f . Then the symmetrization is defined by,

SY M g f ,p [n, k] = 1 2 V S g p, f [n, k] + V S g p, f [n, 2Γ p [n] -k] × 1 Jp [k].
(2.36)

where 1 Jp [k] = 1 if k ∈ J p and 1 Jp [k] = 0 otherwise. Next, SSR-HT smoothes real and imaginary parts of SY M f ,p using piecewise cubic Hermite interpolation, that we denote by SSR g f ,p . Finally, as for LCR and HT, if g[0] ̸ = 0,

f p [n] ≈ 1 g[0]N N 2 -1 k=-N 2 SSR g f ,p [n, k], (2.37)
or assuming f is L-periodic, through (1.8, page 3), if we assume f is L-periodic, g ̸ = 0,

f p [n] ≈ 1 N ∥g∥ 2 N -1 k=0 n+M m=n-M SSR g f ,p [n, k]g[n -m]e i2π k(n-m) N . ( 2

Numerical Applications

In this section, we first show the benefits of LCR when the input SNR is between 0 and 10 dB. At such noise levels, all RD are working fine and RRP-RD does not bring any significant improvement. Therefore, we assume that there are no problems with RD, and we focus on comparing HT, SSR-HT with LCR. Then, we further study HT and LCR techniques, when the input SNR is between -10 and 0 dB, in which cases, RRP-RD is used to detect the ridges. Finally, we show that LCR can be used to solve the mode mixing issue.

Comparison of Mode Retrieval Techniques on a Monocomponent Signal

We study the efficiency of LCR for the signals of We notice that LCR gives better output SNR than HT and SSR-HT in each case, whatever the noise level. We also note that LCR is only slightly sensitive to the technique used for mode reconstruction (M 2 seems a little better than M 1 , but the latter is much faster than the former). On the contrary, SSR-HT behaves worse when the modulation decreases, since the quality of reconstruction is hampered by frequency resolution when M 2 is used. Finally, comparing the results of Figure 2.6 (b) and (c), we notice that LCR is efficient for a wide range of frequency modulation for the mode.

Study of the Robustness of Mode Retrieval

In this section, we investigate the robustness of mode retrieval techniques MR based on S-RD and LCR (S-LCR-MR) and MR based on RRP-RD and LCR (RRP-LCR-MR), and draw a conclusion of the impact of ridge detection on MR.

First, RRP-LCR-MR provides significant improvement compared with other MR techniques that rely on RRP-RD [START_REF] Laurent | A novel ridge detector for nonstationary multicomponent signals: Development and application to robust mode retrieval[END_REF]. Now, we focus on the impact of RD in heavy noise situations. To do so, we consider the signals of output SNR. Even when the noise has a strong impact on the entire STFT, a good ridge detection is necessary to improve MR. In fact, the accuracy of RRP-RD allows to get more accurate IF and CR estimates, by using s f in p instead of ω [START_REF] Baron Fourier | The analytical theory of heat[END_REF] f and (s f in p ) ′ instead of q f . This can be checked numerically, as depicted in Figure 2.8 and justify their usage as IF and CR estimators.

Mode Mixing Issue

In this section, we study the ability offered by HT, SSR-HT and LCR to accurately separate two close modes in a MCS provided they admit two separate ridges. When the modes are close in the TF plane, the modulus of the STFT between the two modes is oscillating, due to interference. Therefore, some coefficients located between two modes cannot be reasonably associated with one of them. This phenomenon is often called mode mixing in the literature. It was studied in the context MR, by comparing empirical mode decomposition and FSST [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], and also in [START_REF] Zheng | Partly ensemble empirical mode decomposition: An improved noise-assisted method for eliminating mode mixing[END_REF], to improve the empirical mode decomposition. In such a case, one may also consider the adaptive transforms [START_REF] Chui | Analysis of an adaptive short-time fourier transform-based multicomponent signal separation method derived from linear chirp local approximation[END_REF], and perform MR based on the signal sepration operator.

To illustrate how we tackle this issue, we consider the two mode signal whose STFT is depicted in Figure 2.9 (a), for which mode separation is challenging. In particular,

we notice that at a low noise level and when using HT or SSR-HT, the threshold used for the determination of intervals J p is lower, at many time instants, than the minimum value of the modulus of STFT between the two ridges associated with the modes. This results in mode mixing with these two techniques. In Figure 2.9 (b), we display the results corresponding to the retrieval of the modes f 1 and f 2 . With HT and SSR-HT techniques, when the noise level is low, we get that for many time indices n, J 1 [n] = J 2 [n], hence the poor results. On the contrary, the separation of the two modes with LCR remains satisfactory whatever the noise level.

As a remark, the reconstruction formula used with HT, SSR-HT, and LCR is (1.7, page 3), to use (1.8, page 3) instead would not change the conclusions.

Application to Gravitational-Wave Signals

In this section, we investigate the applicability of RRP-RD and RRP-LCR-MR to a Because the signal is real, it cannot be analytic and we therefore use the reconstruction formula associated with real signals,

f [n] = 2 g[0]N ℜ    N 2 -1 k=0 V g f [n, k]    . (2.39)
For the sake of comparison, we consider the fourth-order Fourier-based synchrosqueezing transform (FSST4), an efficient reassignment technique introduced in [START_REF] Pham | High-order synchrosqueezing transform for multicomponent signals analysis-with an application to gravitational-wave signal[END_REF]. Note that, as mentioned in Section 1.2, the generalized FSSTN technique will be defined in Chapter 4. We first display in Figure 2.10 (a), the modulus of the STFT of such a signal, along with the ridge associated with FSST4 and the spline obtained using RRP-RD. We notice that RRP-RD and FSST4 leads to very similar results, and that both techniques enable the detection of the ringdown.

We then investigate the quality of signal reconstruction by comparing it with the one given by the numerical relativity [START_REF] Abbott | GW151226: Observation of gravitational waves from a 22-solar-mass binary black hole coalescence[END_REF], when one uses either MR based on RRP- may not be that interesting to use the latter technique instead of RRP-MR in that context. It is also worth noting here that FSST4-MR is very specific to signals containing very oscillatory phases, which is the case of the gravitational wave when the ringdown occurs. For any other modes like those studied before, a lower order synchrosqueezing transform would lead to better results [START_REF] Meignen | On the use of short-time fourier transform and synchrosqueezing-based demodulation for the retrieval of the modes of multicomponent signals[END_REF]. Another limitation of the techniques based on synchrosqueezing transforms is that there are not adapted to long signals since the reconstruction of the modes from these transforms when the hop-size is larger than one is not tractable [START_REF] Meignen | Retrieval of the modes of multicomponent signals from downsampled short-time Fourier transform[END_REF]. Finally, as FSST4-MR uses S-RD, it depends on the user defined parameter B f which is not the case of RRP-MR or RRP-LCR-MR. These last techniques thus offer a nice alternative to the best state-of-the-art techniques for the reconstruction of very complex signals such as gravitational-wave signals.

Conclusion

In this chapter, we proposed a new ridge detection technique and investigated its potential use for MR. We introduced the combination of the so-called RRP-RD and LCR techniques, which leads to satisfactory MR results. However, when considering heavy noise, it becomes a lot more challenging since the noise affects the coefficients of the entire TF plane. In particular, we showed that even with very accurate IF and CR estimates, MR is still challenging and the proposed combination of RRP-RD and LCR is a good answer to that problem. However, it is important to note that our approach is unable to deal with crossing modes, since, in that case, one cannot associate a ridge with each mode. This requires a better understanding of interference, which we will be the topic of Chapter 3. We will also see in Chapter 4, that the noise still has a huge impact on the amplitude of the modes and, therefore an ideal RD is not sufficient for accurate MR.

Time Frequency Analysis of Interference

Introduction

In this chapter, we study the time and frequency interference of a MCS, which arise when two modes are close in the TF plane, thus creating interference patterns.

In a first part of this chapter, we carry out the mathematical study of the conditions for interference patterns called TFBs to appear. TFBs were originally introduced 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] with mathematical characterizations of this phenomenon in the case of two pure harmonics. We also study this situation and extend it to signals made of two parallel linear chirps. Then, because in practice, the conditions for TFBs to appear are difficult to evaluate, and for the purpose of considering general MCSs, we detail a special RD that is able to detect the merging of ridges that coincides with the presence of TFBs.

We then perform a numerical validation of the technique, on noisy MCSs.

In the second part of this chapter, we consider a different approach to detect interference patterns based on the analysis of the zeros of the spectrogram, and show perspectives on how to perform a separation between different types of zeros.

On the Existence of Time Frequency Bubbles

In this section, we investigate the number of existing local maxima of the spectrogram along the frequency axis (LMFs) when two modes are interfering. Depending on the situation in time and frequency, we show that when modes are merging or are too close in the TF plane, there may only be one LMF associated with both modes.

In what follows, to characterize these situations, we will use the notion of TFB, which corresponds to the following definition:

This chapter is based on the work by S. Meignen, T. Oberlin, and N. Laurent, "On the Study of Interference and their Localization in the Time-Frequency Plane" submitted to IEEE TSP.

Definition 3.2.1. Two modes make a TFB in the TF plane when the LMFs associated with each mode merge.

We now investigate the existence of TFBs when the signal is the sum of two pure tones and then explain how the obtained results can be extended to the case of two parallel linear chirps. As explained in Chapter 2, a Gaussian window is well adapted for such a study as it allows working with closed form expressions thus enabling the mathematical analysis. Therefore, we set the window g(t) = e -π t 2 σ 2 .

Case of the Sum of Two Pure Tones

Let

f (t) = f 1 (t) + f 2 (t) with f 1 (t) = Ae i2πξ 1 t and f 2 (t) = e i2πξ 2 t
, where ξ 1 < ξ 2 . As one

has V g f 1 (t, η) = σAe i2πξ 1 t e -π(η-ξ 1 ) 2 σ 2 and V g f 2 (t, η) = σe i2πξ 2 t e -π(η-ξ 2 ) 2 σ 2 , |V g 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.1)
A first property of |V g f (., η)| 2 , where . means we consider this variable, is that it attains its maximum at t k = k ξ 2 -ξ 1 , k ∈ Z for which:

|V g f (t k , η)| 2 = σ 2 (Ae -πσ 2 (η-ξ 1 ) 2 + e -πσ 2 (η-ξ 2 ) 2 ) 2 . ( 3
| log(A)| < -2 arccosh(α) + 2α α 2 -1.
In any other case, |V g f (t k , .)| 2 has a unique extremum which is a LMF.

The proof is given in Appendix A.

:= π 2 σ(ξ 2 -ξ 1 ) > 1 and | log(A)| < -2 arccosh(α) + 2α α 2 -1.

ON THE EXISTENCE OF TIME FREQUENCY BUBBLES

In any other case, there exists some time t where |V g f (t, .)| 2 has a unique extremum which is a LMF.

The proof is given in Appendix A.2. This proposition means that, if the spectrogram exhibits two LMFs at t k , it is the same for any time t, and it tells us when two pure tones are generating TFBs.

Note that, when A = 1 and since -2 arccosh(α) + 2α √ α 2 -1 is an increasing function which equals 0 at 1, when α > 1 we always have -2 arccosh(α)

+ 2α √ α 2 -1 > 0,

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 ϕ k (t) = a k t + C 2 t 2 with C constant.
As one has a similar formula to (2.24, page 24),

V g f 1 (t, η) = f 1 (t)r -1/2 e - πσ 2 (1+iCσ 2 )(η-ϕ ′ 1 (t)) 2 1+C 2 σ 4 e -i θ 2 (3.3) with r = √ 1+C 2 σ 4 σ 2
and θ = atan(-Cσ 2 ), one obtains: with

|V g f (t, η)| 2 = 1 r A 2 e -2π (η-ϕ ′ 1 (t)) 2 r 2 σ 2 + e -2π (η-ϕ ′ 2 (t)) 2 r 2 σ 2 + 2Ae -π [ (η-ϕ ′ 1 (t)) 2 +(η-ϕ ′ 2 (t)) 2 ]
ϕ(t, η) = (a 2 -a 1 ) t + C r 2 σ 2 η -Ct -a 1 +a 2 2
. Let us then consider the set of TF points (t, η) such that cos(2πϕ(t, η)) = λ. These points correspond to

ϕ(t, η) = arccos(λ) 2π + k for k ∈ Z ⇔ t = Dη + B λ k (3.5) with D = C C 2 -r 2 σ 2 and B λ k = C 2 (a 2 2 -a 2 1 ) + r 2 σ 2 ( arccos(λ) π + 2k) 2(r 2 σ 2 -C 2 )(a 2 -a 1 )
.

We deduce from (3.4) that:

|V g f (Dη + B 1 k , η)| 2 = 1 r Ae -π( 1-CD rσ ) 2 (η- a 1 +CB 1 k 1-CD ) 2 + e -π( 1-CD rσ ) 2 (η- a 2 +CB 1 k 1-CD ) 2 2 , (3.6)
which is the same expression as (3.2), replacing σ, ξ 1 and ξ 2 by 1-CD rσ ,

a 1 +CB 1 k 1-CD and
| log(A)| < -2 arccosh(α C ) + 2α C α 2 C -1,
and a unique maximum along that line otherwise.

The proof is the same as that of Proposition 3.

| log(A)| < -2 arccosh(α C ) + 2α C α 2 C -1,
and a unique maximum along that line otherwise.

So the presence of TFBs is here related to what happens when λ = 1. Let us denote by

(t k , η k ) a local maximum on the line t = Dη + B 1 k and then a neighboring point (t, η) of (t k , η k ) in the direction of vector (1, C). For such a point, one has η -ϕ ′ 1 (t) = η k -ϕ ′ 1 (t k ) and η -ϕ ′ 2 (t) = η k -ϕ ′ 2 (t k
) and, as cos(2πϕ(t, η)) < 1, one gets that (t k , η k ) is also a maximum in the direction of vector (1, C), and thus a local maximum of the bidimensional function |V g f (t, η)| 2 . Note finally that, contrary to what happens in the pure tones case, (1, C) is not orthogonal to (D, 1).

An illustration is given in Figure 3.2 (a) where the spectrogram of two parallel linear chirps is displayed and for which there exist three extrema along direction t = Dη + B λ k . Then, in Figure 3.2 (b), the case where there exists only one maximum for some value of λ is displayed (TFBs are thus present). Finally, in Figure 3.2 (c) one considers a case where there are three extrema along direction t = Dη + B λ k when A is no longer equal to 1. In each case, the zeros and the local maxima of the spectrogram are also plotted.

Localizing TFBs Based on a New RD

In this section, we design a technique for the detection of TFBs. Indeed, the conditions we expressed above require the knowledge on the frequency and modulation of the signal, which is not the case in practice.

From the definition of TFBs, we propose to define a new RD to find the LMFs that are merging. However, the previously introduced RD such as S-RD in Section 2.2.1, MB-RD in Section 2.2.2 and RRP-RD in Section 2.3.3 are not compatible for such a task, because they do not allow ridges to cross. Therefore, in the section below, we explain how to build a new RD enabling mode crossing.

New RD Enabling the Merging of Ridges

In Chapter 2, the RDs we introduced are specifically designed to work for non-crossing modes and cannot be used as is to detect TFBs. To build the new RD enabling the merging of ridges, and because we only consider input SNR above 5dB, we propose to modify S-RD introduced in Section 2.2.1, which gives good results in such situations.

In the case of S-RD, the ridge Γ p defined by means of (2.6, page 13) cannot cross the previously detected ridges, and is not necessarily made of LMFs. To circumvent these two limitations of S-RD, we still use the same initialization process, namely, for an

Γ p [n 0 ] = arg max 0≤k≤N -1 |V g f [n 0 , k]|, (3.7)
but replace (2.6, page 13) by

Γ p [n + 1] := argmin k |k -Γ p [n]| ≤ B f , s.t.|V g f [n + 1, k]| LM F . (3.8)
This has two straightforward consequences. First, it is possible that Γ p merges with some of the previously computed ridges, and second, the ridge points are all LMFs. As a consequence, if no LMF is found in the interval of interest, the construction of the ridge is stopped, which was not the case with S-RD. The second consequence is that, contrary to S-RD, when two LMFs are found in the interval of interest one considers the closest to the current position rather than the global maximum on the interval.

Note that to extract a ridge with S-RD, the procedure is run with different initializations and the ridge associated with the largest energy is kept. When using (3.8) to compute Γ p instead, the latter may have points in common with the previously computed ridges, thus the energy we associate with Γ p corresponds this time to the sum of the spectrogram coefficients that do not belong to the previously computed ridges. Meanwhile we also define the length of Γ p as the length of the ridge portions belonging only to Γ p and not to the previously computed ridges. Thus, in that context, we no longer perform the ridge detection until P ridges are detected, but until a small fraction of the energy is not contained in the vicinity of the detected ridges. In what follows, this RD will be called redundant ridge detector (R-RD), to emphasize the fact that some ridges may have some points in common.

An illustration of the practical impact of using R-RD instead of S-RD is given in In the presence of noise, R-RD is slightly adapted so as to ensure the ridges are made of significant LMFs. For that purpose, let f = f + ε be a noisy signal, with ε a zero mean complex Gaussian white noise. To find out the regions where significant coefficients are located, we consider the hard-thresholded spectrogram S(3) that we defined in 2.11.

This means that we consider coefficients such that |V g f [n, k]| > 3γ where γ is an estimate of σ∥h∥ 2 . In such a noisy context, to build a ridge, one considers the same initialization procedure as previously replacing f by f and maximizing only over TF points in S(3), and then the construction follows (3.8) considering only the TF points that are in S(3). All the ridges are extracted iteratively, and the extraction is stopped when the last computed ridge has a length below some threshold. When some noise is present, it is more efficient to put a bound on the ridge length rather than on the remaining energy away from the ridges since, due to noise, this energy can be large. This RD will be called redundant ridge detector in noise (R-RDN) in the sequel.

Using New RDs to Detect and Localize TFBs

Based on the previous analysis, a simple criterion to investigate whether the proposed new RDs are efficient to detect TFBs is to find out the instances when some ridges have merged. To measure this, let us consider that either with R-RD or R-RDN, depending on whether the signal is noisy or not, K ridges (Γ k ) k=1,••• ,K have been extracted. Then, let us define the following quantity:

R(σ) =    1, if min n,1≤i<j≤K (|Γ i [n] -Γ j [n]|) = 0 0, otherwise.
(3.9) Thus, R(σ) = 1 when some ridges have some points in common, and it is equal to zero otherwise. We will see in the Results section that the values of R computed from the extracted ridges are in accordance with the necessary and sufficient conditions on the existence of TFBs given by the theoretical analysis of Section 3.2.

LOCALIZING TFBS BASED ON A NEW RD

Once we know two ridges have some points in common in the TF plane, we spot where the merging occurs, the position of which being an estimate of the location of a TFB. An illustration of the procedure is given in Figure 3.4, first row and in the absence of noise (so we use R-RD for ridge detection). The STFT in Figure 3.4 (a) to (c) are computed with σ = 0.02, 0.03, and 0.04 respectively. The stars on these figures indicate where some ridges have merged. It is interesting to notice that, as expected, a smaller σ creates more interference between the modes and thus more TFBs. In that case, by increasing σ, it is possible to get three disjoint ridges (see Figure 3.4 (c)), telling us that the separation of the modes will be easier with such a σ. Then, on the second row of Figure 3.4, we give an illustration in the presence of noise, considering the same signal as on the first row, and with σ = 0.04 as we know such a value leads to three different ridges in the absence of noise. In that context, we use R-RDN for ridge detection, and consider different noise levels (input SNR of 20, 10 and 5 dB) corresponding to Figure 3.4 (d) to (f). In these figures, we plot TF regions corresponding to S(3) and then the ridges computed in these regions using R-RDN. We notice that for the studied noise, the merging of the ridges are located in the same regions as in the noiseless case, when the input SNR is larger than 5 dB, which means the detection of TFBs using the proposed approach seems to be robust to noise. Nevertheless, at this stage, no general conclusions can be drawn since only one realization of the noise is considered, therefore the quality of the detection is going to be investigated more in details in the following section.

Results

In this section, we first investigate how the proposed new RDs enable us to detect TFBs.

Then, we study whether the localization of the merging of ridges computed with our algorithm is a robust criterion to localize TFBs. Finally we propose an application of the potential use of the localization of TFBs for the analysis of vibrato signals.

3.4. RESULTS

TFB Detection Based on new RDs

We investigate here how well the merging of ridges with the RDs we propose reflects the presence of TFBs in the TF plane. To do so, we first consider the case of the sum of two pure tones and then the case of two parallel linear chirps. For that purpose we use the signal of Figure 3.1 (a), which corresponds to the sum of two pure tones with the same amplitude. In such a case, we know ξ 2 -ξ 1 and we make σ vary: TFBs are present when α = π 2 σ(ξ 2 -ξ 1 ) ≤ 1 and not otherwise. This defines a function of σ equal to 1 when a TFB is present and to 0 otherwise. This function is denoted by ground truth (GT) hereafter. Then, from (3.9), R(σ) = 1 when some ridges have some points in common and is equal to zero otherwise. So R should be close to GT if it is relevant to detect the presence of TFBs. To check this, In Figure 3.5 (a), we display R along with GT, in noiseless and noisy cases, R-RD and R-RDN being used respectively to compute the different ridges. We notice that the appearance of TFBs corresponds to σ below a certain threshold, which is computed with a great accuracy using (3.9) for input SNR larger than 20 dB. This threshold is slightly increased when an heavier noise is added to the signal, meaning the noise creates some TFBs that are not present in the noiseless case. R, for the just mentioned signal, and we notice that the appearance of TFBs is detected for a very similar value of σ than that given by the theory when little noise is added. As in the previous case, we notice that, for some values of σ, the noise creates some TFBs that are not present in the noiseless case. Note finally that when (3.9) is used to estimate the presence of TFBs when the signal is made of two parallel linear chirps, similar results to those displayed in Figure 3.5 are obtained.

Now, if one considers the signal of

The question we would like to address is whether the TFBs generated by noise is detrimental to the performance of the procedure we propose to localize the TFBs of the noiseless signal.

TFBs Localization Based on Ridge Merging

We hypothesized in the previous section that the locations of the merging of ridges were a good indicator of those of TFBs, even in noisy situations. The results depicted on the second row of Figure 3.4 were encouraging in that matter, but should be generalized to many different realizations of the noise to conclude. For that purpose, we again consider the signal studied in Figure 3.4, and compute the STFT with σ = 0.04, for which, in the absence of noise, the merging of ridges only occurs between the pure tone signal and the linear chirp. To investigate whether this is still the case in the presence of noise, we compute two masks corresponding to where the interference between two modes f 1 and f 2 are significant, namely we define the region of interference of level T as:

I T (f 1 , f 2 ) = [n, k], s.t.|ℜ(V g f 1 [n, k]V g f 2 [n, k] * )| ≥ T . ( 3

Application to Voice Signals

In this section, our goal is to give an illustration of the potential interest of the procedure we propose to localize TFBs in practical situations, and in particular for the study of voice signals. Indeed, for such signals, it is of interest to have access to its different harmonics, which can be done by extracting the ridges associated with them in the TF plane. A natural question one then asks is whether this ridge extraction can be trusted.

We are going to show that the procedure we propose to localize TFBs offers, in some

way, an answer to that question.

To illustrate this, we consider the voice signal of a 22-year-old female taken from the recording 62 of the Saarbruken Voice Database [START_REF] Woldert-Jokisz | Saarbruecken voice database[END_REF]. This signal contains 126902 samples for a length of 2.538 seconds and for efficiency purpose, we resample it by a factor of 10.

It contains sustained /a/ vowels of the type "low-high-low" meaning there is a change on the pitch of the signal. The STFT of such a signal is depicted in

Analysis of Interference Based on the Zeros of the

Spectrogram

In this section, we consider an approach based on the zeros of the spectrogram to localize interference. In some sense, this approach can be seen as dual to the one proposed just above, in that ridges are often associated with high energy TF regions.

The zeros have already been used in TF analysis, for instance, in the context of the auditory system [START_REF] Gardner | Sparse time-frequency representations[END_REF], where, using reassignment, the repulsive effect of the zeros is illustrated, or in the context of [START_REF] Balazs | The pole behavior of the phase derivative of the short-time fourier transform[END_REF] where the relation between the poles of the phase with the location of the zeros is investigated.

Then in [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF], it is explained that the STFT can be viewed as a function of the complex plane, using the Bargmann factorization of the STFT [START_REF] Gröchenig | Foundations of time-frequency analysis[END_REF]Section 3.4]. It can in fact be regarded as an entire function that is completely characterized by its zeros, which, by construction [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF], are the zeros of the STFT.

The Delaunay triangulation based on the spectrogram's zeros is a TFR that can be used for denoising and MR purposes [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF]. Succinctly explained, since the distance between two neighboring zeros should not exceed a certain value [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF], it is proposed in that paper to consider the Delaunay triangulation of the zeros and then to discriminate the zeros related to the signal from those related to noise by considering a threshold on the length of the triangle edges. To build the signal domain, the following simple procedure is considered [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF]: The shape of triangles associated with noise have a specific distribution.

When a triangle is an outlier for this distribution, it is associated with signal. Adjacent outliers are then considered belonging to the same mode, enabling MR.

Further mathematical description is given in [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF], on statistical properties, notably the link between zeros and Gaussian analytic functions. When one considers a complex and i.i.d. white Gaussian noise ε, the zeros of its spectrogram are the set of TF coefficients satisfying |V g ε (t, η)| 2 = 0, and their distribution is uniform on the TF plane [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF].

Perspectives for a Zero Based Technique to Identify Interference

In this section, we study how to identify the signal domain in the TF plane, using the zeros of the spectrogram. As explained in the previous section, this has been done in [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF],

but in our case, we also want to identify interference in the TF plane,meaning we want to split the TF plane into three categories associated with either noise, interference or signal.

We also consider a different approach, where we study the Voronoi cells associated with the zeros instead of their Delaunay triangles. To illustrate the perspectives described in this section, we consider the signal f = f 1 + f 2 + ε, where f 1 and f 2 are two pure harmonics (PHs) that interfere with each other, and ε is a complex and i.i.d. white Gaussian noise. In Figure 3.9 (a), we depict the Voronoi cells of the zeros of the modulus of the STFT, namely |V g f |, where g is a Gaussian window. In order to separate the TF plane, we aim at identifying three types of zeros associated with one category. We propose to discriminate them using the distribution of energy in their Voronoi cells; such a cell is denoted by C i f , and is a polygon composed of Taking into account these observations, our goal is to study the maximum energy on a segment C i,j f , defined by In order to identify the type of cell C i f , we assume that there exist a threshold E such that if M i,j

J i segments C i,j f for j ∈ {1, • • • , J i }.
M i,j f = max (t,η)∈C i,j f |V g f (t, η)| 2 . ( 3
f > E, then C i,j f crosses the signal domain. Then, by defining

E max i = max j∈{1,••• ,J i } M i,j f (3.12)
E min i = min j∈{1,••• ,J i } M i,j f , (3.13)
we propose to separate the cells by utilizing the following scheme.

If

E min i > E, C i f is related to interference, if E max i ≤ E, C i
f is related to noise and otherwise, C i f is related to signal. However, the choice of E is delicate, and to set it in a meaningful way, we propose to study the distribution M i,j ε , meaning we only consider pure noise.

A similar problem was dealt with in [START_REF] Flandrin | On spectrogram local maxima[END_REF], where the distribution of the spectrogram maxima was identified to follow a Gumbel law. In our case, we restrict our selection of the maximum in a segment of a cell. Namely, first look to determine the law associated with M i,j f . Using numerical computation, we depict the histogram of its probability density function in Figure 3.10 and plan identify the law and parameters in a near future.

Conclusion

In this chapter, we studied two different approaches on interference. The first is based on TFBs, for which we carried out the study of the mathematical conditions for their appearance. Then, we introduced a technique to detect TFBs for a general MCS and 52 3.7. CONCLUSION use numerical validation to show its accuracy in noisy conditions. We also performed its validation on a real signal, which illustrates that the technique can be used as an indicator of the ridge detection quality.

In the last part of the chapter, we study a different approach based on the zeros of the spectrogram. We proposed a scheme that uses the distribution of the energy around the Voronoi cells associated with the zeros to determine TF regions associated with noise, interference or signal. Assuming that the maxima of the energy in the border of a Voronoi cells and the maxima of the spectrogram are similarly distributed, a new technique could be designed to separate the TF regions, which we plan to investigate in the future.

Towards Adaptive and Robust Synchrosqueezing

Introduction

In this chapter, we first recall that to get closer to the ideal time frequency representation when the modes of a MCS have strong frequency modulation, FSST2 was generalized to any order N through FSSTN [START_REF] Pham | High-order synchrosqueezing transform for multicomponent signals analysis-with an application to gravitational-wave signal[END_REF]. While to use a N larger than 2 can be profitable in some instances, it is hard to determine which N is the best in a particular situation.

Ideally, one would seek to find N in accordance with the order of the phase [START_REF] Pham | Demodulation algorithm based on higher order synchrosqueezing[END_REF], but the latter is unknown in general. Furthermore, FSSTN, with N large, is sensitive to noise [START_REF] Pham | Demodulation algorithm based on higher order synchrosqueezing[END_REF]. In this chapter, we design a technique that aims to determine the best N locally in the TF plane. Our approach is based on energy concentration, which we interpret as the quality of the reassignment performed by FSSTN. Note that in [START_REF] Daubechies | ConceFT: Concentration of frequency and time via a multitapered synchrosqueezed transform[END_REF], a multitapered wavelet based SST was introduced to obtain better concentrated TFR and represents an alternative to our technique.

Definition of FSSTN

The objective of this section is to recall the definition of FSSTN, introduced in [START_REF] Meignen | Synchrosqueezing transforms: From low-to high-frequency modulations and perspectives[END_REF][START_REF] Pham | High-order synchrosqueezing transform for multicomponent signals analysis-with an application to gravitational-wave signal[END_REF].

To motivate this work, we mention that even though FSST2 leads to a concentrated TFR for linear chirps with Gaussian modulated amplitudes, it is not perfect for a general MCS.

For instance, if a mode has a polynomial phase of order greater than two, then one needs a transform of higher order to obtain an accurate reassigned transform, this is precisely the aim of FSSTN.

The definition of FSSTN is based on a Taylor expansion of the signal's amplitude and phase. Let us consider the single mode signal f (t) := A(t)e 2iπϕ(t) , where we assume

This chapter is based on the work published in [START_REF] Laurent | A new adaptive technique for multicomponent signals reassignment based on synchrosqueezing transform[END_REF] (to appear).

log(A) and ϕ to be polynomial of degree N . Let t ∈ R, the N -th order expansion of A and ϕ at τ is

f (τ ) = exp N k=0 (τ -t) k k! [log(A(t))] (k) + 2iπϕ (k) (t) . (4.1)
Considering g ∈ L 2 (R), if f ∈ L 2 (R), its STFT satisfies, V g f (t, η) = R exp N k=0 τ k k! [log(A(t))] (k) + 2iπϕ (k) (t) g(τ )e -2iπτ η dτ (4.2)
By differentiating with respect to t and then by dividing by 2πV g f one obtains,

ωf (t, η) := ∂ t V g f (t, η) 2iπV g f (t, η) = N k=1 r [N ] k (t) V t k-1 g f (t, η) V g f (t, η) , (4.3)
where

r [N ] k (t) = 1 (k -1)! [log(A(t))] (k) 2iπ + ϕ (k) (t) . (4.4)
In order to define the N -th order local IF estimate, we first set,

ω[N] f (t, η) = ωf (t, η) - N k=2 r [N ] k (t) V t k-1 g f (t, η) V g f (t, η) = r [N ] 1 (t) (4.5)
and then define ω

[N ] f (t, η) := ℜ{ω [N ]
f (t, η)} where in our case, ω

[N ] f (t, η) = ϕ ′ (t). A simple way to compute r [N ] 1
is first, using (4.3), to remark that

1 2iπ ∂ t V g f (t, η) = N k=1 r [N ] k (t)V t k-1 g f (t, η), (4.6)
and that

∂ η V g f (t, η) = -2iπV tg f (t, η).
Then, we deduce the following equality [START_REF] Meignen | Analysis of reassignment operators used in synchrosqueezing transforms: With an application to instantaneous frequency estimation[END_REF],

1 2iπ        ∂ t V g f i 2π ∂ η ∂ t V g f . . . i N -1 (2π) N -1 ∂ N -1 η ∂ t V g f        =        V g f V tg f • • • V t N -1 g f V tg f V t 2 g f • • • V t N g f . . . . . . . . . . . . V t N -1 g f V t N g f • • • V t 2(N -1) g f              r [N ] 1 r [N ] 2 . . . r [N ] N       = DR. (4.7) 4.2

. DEFINITION OF FSSTN

Finally, based on simple properties of the determinant of matrices, setting,

M 1 =        ∂ t V g f V tg f • • • V t N -1 g f i 2π ∂ η ∂ t V g f V t 2 g f • • • V t N g f . . . . . . . . . . . . i N -1 (2π) N -1 ∂ N -1 η ∂ t V g f V t N g f • • • V t 2(N -1) g f        , (4.8)
one obtains

r [N ] 1 = 1 2iπ det(D -1 ) det(M 1 ) = 1 2iπ det(M 1 ) det(D) . (4.9)
Then, as

∂ t V g f = i2πηV g f -V g ′ f ,
one gets, for any l ≥ 1:

∂ l η ∂ t V g f = (-2iπ) l -kV t l-1 g f -V t l g ′ f + 2iπηV t l g f , (4.10)
leading to det(M 1 ) = i2πη det(D) -det(U 1 ) -det(V 1 ) with U 1 =        0 V tg f • • • V t N -1 g f V g f V t 2 g f • • • V t N g f . . . . . . . . . . . . (N -1)V t N -2 g f V t N g f • • • V t 2(N -1) g f        , (4.11)
V 1 =        V g ′ f V tg f • • • V t N -1 g f V tg ′ f V t 2 g f • • • V t N g f . . . . . . . . . . . . V t N -1 g ′ f V t N g f • • • V t 2(N -1) g f        , (4.12)
and thus

ω [N ] f = ℜ r [N ] 1 = η - 1 2π ℑ det(U 1 ) + det(V 1 ) det(D) . (4.13)
Finally, FSSTN is defined by [START_REF] Meignen | Synchrosqueezing transforms: From low-to high-frequency modulations and perspectives[END_REF][START_REF] Pham | High-order synchrosqueezing transform for multicomponent signals analysis-with an application to gravitational-wave signal[END_REF],

T g,N f (t, η) := R + V g f (t, τ )δ η -ω [N ] f (t, τ ) dτ. ( 4
f (t) = 1 g(0) R T g,N f (t, η)dη. (4.15)

Adaptive selection of the order

The objective of this section is to set the order N of FSSTN locally, to get a sharp TFR.

In the general case, one does not know the signal's modulation, and therefore it is not possible to set the order of FSSTN to the degree of the phase. An adaptive way to set N would be to measure the concentration of energy, and find N that maximizes it.

Let us first consider the general case of a noisy MCS affected by the same kind of white Gaussian noise as in Chapter 2.

f = f + ε (4.16)
where f is a MCS in L 1 (R) ∩ L 2 (R) of P modes, and ε a white Gaussian noise. Our interest is to concentrate the energy to ease the reading of the signal in the TF plane, and for that purpose, we consider that for each mode p, its ridge Γ p is given by a RD. We then measure the concentration of energy on the ridge in T g,N f . In the ideal case, considering a monocomponent signal f (t), based on the reconstruction formula of FSSTN, all the information is located on the ridge,

f (t) ≈ T g,N f (t, Γ(t)) g(0) . (4.17)
Based on that, we measure the energy concentration using the following estimate of the instantaneous amplitude, when g(0) > 0,

E p N (t) = T g,N f (t, Γ p (t)) g(0) , (4.18)
where g is a window in

L 1 (R) ∩ L 2 (R) and RD is performed on T g,N f
. With a continuous ridge, where in the ideal case Γ p (t) = ϕ ′ p (t), the selection of N could be made by maximizing (4.18), by introducing

N p (t) = arg max N ∈N E p N (t). (4.19)
For the sake of evaluating the relevance of this method, using the invertibility of FSSTN,

ADAPTIVE SELECTION OF THE ORDER

we define,

R N p fp (t) := T g,N p (t) f (t, Γ p (t)) g(0) , (4.20)
which is an estimate of the p-th mode, namely f p (t) ≈ R N p fp (t). This enables us to make the comparison of this technique with the signal in terms of SNR. Before that, in order to proceed with the practical study of this method, we first detail the corresponding discrete framework.

We consider the sequence f

[n] = f ( n L ), where n ∈ {1, • • • , L}.
The discrete FSSTN definition is given by,

T g,N f [n, w] = J k=1 V g f [n, k]δ{w, ⌊ ω [N ] f [n, k]⌉} (4.21)
where ⌊x⌉ is the closest integer to x ∈ R and δ is the Kronecker symbol. We also recall that the discrete version of ω

[N ] f used in 4.21 only uses the discrete formalism of the STFT.

However, by using this transform directly, instabilities related to the noise and numerical computations may significantly worsen the quality of mode reconstruction with our method. Indeed, when ω

[N ] f [n, k]
is located between two frequency bins, it is likely that the energy will be dispatched in the upper and lower bins. Obviously, this will impact the measure given by (4.18). To avoid this effect, we consider a variant of T g,N f , where we shift the frequency grid by q ∈ {0, 1 2 } of a frequency bin and choose q that maximizes the concentration. For that purpose, we define the shifted variant of FSSTN by,

T g,N f ,q [n, w] := J k=1 V g f [n, k]δ{w, ⌊ ω [N ] f [n, k] + q⌉}, (4.22)
where J is the number of frequency bins. Then a second issue is that we cannot compute FSSTN for all N ∈ N, and therefore we set N 0 ∈ N and consider orders up to N 0 . Taking the shift into account, we define the energy concentration for

N ∈ {1, • • • N 0 } in the discrete framework by, E p N [n] := max q∈{0, 1 2 } T g,N f ,q [n, Γ p [n]] g(0) , (4.23)
which represents an estimate of the instantaneous amplitude in the discrete case. Note that the ridge given by Γ p is also computed on T g,N f ,q , and that as a consequence, the shift is fixed independently for each mode.

CHAPTER 4. TOWARDS ADAPTIVE AND ROBUST SYNCHROSQUEEZING

In that context, we set

N p [n] = arg max N ∈{1,••• ,N 0 } E p N [n]
. This allows us to set the proper choice of q depending on n using,

Q p [n] := arg max q∈{0, 1 2 } T g,N f ,q [n, Γ p [n]] . (4.24)
Finally, to approximate the mode, as in (4.20), we set,

R N p fp [n] = 1 g(0) T g,N p (t) f ,Qp[n] [n, Γ p [n]]. (4.25)
To draw conclusions on the quality of our method, we can now compute the output SNR using this estimate of f p .

Numerical Results for the Monocomponent Case

In this subsection, we consider monocomponent signals with constant amplitude A = 1. We first look at amplitude estimates given by E N , for N ∈ {2, • • • , 4}. Note that first order is not considered, as it does not take any modulation into account, and is significantly worse than other orders.

In Figure 4.1 (a), the signal has a high-order and oscillating frequency. In Figure 4.3, we proceed with the comparison between the adaptive selection of the order for R N f and R N f measured in terms of output SNR associated with mode reconstruction when the level of noise varies. It is clear that the adaptive technique allows good reconstruction quality. We also note that when the noise level is strong, R N f gives similar results to FSST2.

We next focus on the results for MCSs.

Numerical Results for the Multicomponent Case

In this subsection, we illustrate the benefits of adaptive selection of the orders for each mode.

For that purpose, we consider the signal of means that an adaptive selection of the order depending on time is beneficial compared with using FSSTN with a fixed order N . The rationale for such a behavior is that with, in the case of MCSs, the choice of the order can be set independently for each mode.

Finally, we also note that when the input SNR is around 10 dB or below, fixing N = 2, namely considering R 2 fp is competitive, like in the monocomponent case presented in previous the subsection.

Application to a Bat Echolocation

We now investigate how the proposed adaptive technique improves the quality of reconstruction on a bat echolocation call whose STFT magnitude is depicted in Figure 4.5 (a).

In our study, we only look for the three most energetic modes in ridge detection and put

Conclusion

In this chapter, we defined an indicator of energy concentration that we used in order to select the appropriate order of FSSTN to improve TFR readability. To assess this improvement, we measured the quality of mode reconstruction considering a fixed or adaptive order in synchrosqueezing transforms. Numerical results showed that the adaptive technique is of better quality in general, but that FSST2 is also a good choice in very noisy conditions.

Study of the Chirp Rate Estimate in Case of Noise

In this chapter, the objective is to estimate the CR, in a way that is robust to noise.

Such an estimate is of interest, notably in synthetic aperture radar and inverse synthetic aperture radar imaging [START_REF] Xia | Discrete chirp-fourier transform and its application to chirp rate estimation[END_REF], and in voice signals to study the jitter phenomenon [START_REF] Baken | Clinical measurement of speech and voice[END_REF][START_REF] Miramont | Voice jitter estimation using high-order synchrosqueezing operators[END_REF].

We focus on the CR estimator in FSST2 and first study its behavior in noisy conditions.

In the previous chapter, the numerical simulations showed that the second order is relevant to consider when the noise is strong, and enables good energy concentration, and that to consider a higher order for FSSTN does not improve the result when the input SNR is around 10 dB or below. When one performs a study in very noisy conditions, FSST2 is a valid choice of representation. However, in Chapter 2, we highlighted that this estimate is sensitive to noise and required further investigations. With that in mind, we study how the noise damages the CR estimate and show that the estimation error is oscillatory. Then, we detail a technique to filter out most of these oscillations using the power spectral density of the error. We finally evaluate the quality of the filtered estimate and its stability to noise on a real voice signal.

To carry out this study, a general definition of the STFT is required, because the noise is not a Lebesgue-integrable function. With a proper definition, the STFT of the noise can be considered as a random variable. For that purpose, we will consider that ε is a tempered distribution, namely ε ∈ S ′ (R), and define the STFT by [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF],

V g ε (t, η) = R g(τ -t)e -2iπ(τ -t)η dB ε (τ ), (5.1)
where g ∈ S(R). The right-hand side is an Itô integral with respect to the one dimensional Brownian motion B ε (τ ). The derivative of Brownian motion is a white Gaussian noise and therefore one can think of dB ε (τ ) as ε(τ )dτ . It follows that V g ε (t, η) has a covariance given by the 2-by-2 square matrix σ 2 ε ∥g∥ 2 2 I 2 [START_REF] Holden | Stochastic partial differential equations[END_REF].

We now proceed with the analysis of the CR estimate for a monocomponent signal f , setting f = f + ε and g(t) = e -π t 2 σ 2 . Remarking that g ∈ S(R) and using the properties of g, we compute a simpler expression of q f . For the sake of clarity, we omit timefrequency coordinates (t, η) in the most complicated expressions, and first remark that in our context, one has

q f = - 1 2π ℑ    (V g f ) 2 V g f V t 2 g f -(V tg f ) 2    . (5.2)
Indeed, in equation (1.17, page 4) with our choice for g, the STFTs depending on g ′ can be simplified using

V g ′ f (t, η) = -2 π σ 2 V tg f (t, η).
Having a simple expression of q f , our interest is now to denoise it.

Study of a Simplified CR Estimate

Now, we further simplify q f by identifying negligible terms and thus defining new estimators. In Section 2.5, we detailed that the STFT of a linear chirp is a Gaussian centered at ϕ ′ (t) with respect to frequency. This means that its frequency derivative vanishes at

ϕ ′ (t). Recalling that, -∂ η 2iπ V g f (t, η) = V tg f (t, η), (5.3)
and assuming that f can be approximated by a linear chirp, we use that V tg f (t, ϕ ′ (t)) ≈ 0, and simplify the operator by making a Taylor expansion, in the vicinity of the ridge.

q f = - 1 2π ℑ        V g f V t 2 g f 1 1 - (V tg f ) 2 V g f V t 2 g f        ≈ - 1 2π ℑ V g f V t 2 g f . (5.4)
From (5.4), we define the new estimate based on the last approximation by,

qf := - 1 2π ℑ V g f V t 2 g f . (5.5)
Taking the noise into account, one gets,

qf+ε = - 1 2π ℑ V g f + V g ε V t 2 g f + V t 2 g ε . . (5.6)
Numerically, we show that qf+ε varies similarly to q f +ε for the three signals of Figure 5.1 (a), (b) and (c), as illustrated on Figure 5.1 (d), (e) and (f). We also notice that the quality of the estimation changes with the modulation of the signal. In order to denoise qf+ε , we are going to study it as a random variable. We first study the bias of ℜ{q f +ε (t, η)} as an estimate of ϕ ′′ (t).

Proposition 5.0.1. At a STFT ridge point (t, ϕ ′ (t)) of a linear chirp one has:

E {q f +ε } = E - 1 2π ℑ {Z} = ϕ ′′ (t)   1 -e - |V t 2 g f | 2 σ 2 ε ∥t 2 g∥ 2    ,
where σ 2 ε is the variance of the noise and Z(t, η) =

V g f +ε (t,η) V t 2 g f +ε (t,η)
.

The proof is given in Appendix B.1. We remark that the bias depends on V t 2 g f , which we develop below to study its effect. In the general case, one has [77]

V t k g f (t, η) = i 2π k ∂ k η V g f (t, η) = -σ 2 (1 + iϕ ′′ (t)σ 2 ) 1 + (ϕ ′′ (t)σ 2 ) 2 i(η -ϕ ′ (t))V t k-1 g f (t, η) - k -1 2π V t k-2 g f (t, η) , (5.7)
and thus,

V t k g f (t, ϕ ′ (t)) = σ 2 (1 + iϕ ′′ (t)σ 2 ) 1 + (ϕ ′′ (t)σ 2 ) 2 k -1 2π V t k-2 g f (t, ϕ ′ (t)), (5.8)
which implies that [19, proposition 1] (omitting (t, ϕ ′ (t)) for the sake of simplicity),

|V t 2 g f | = σ 2 V g f 2π 1 + (ϕ ′′ (t)σ 2 ) 2 = Aσ 2 2π(1 + (ϕ ′′ (t)σ 2 ) 2 ) 3/4 .
(5.9)

Then since ∥t 2 g∥ 2 = σ 5 τ 4 e -2πτ 2 dτ = 3σ 5 (2π) 2 4 √ 2 , we deduce that,

|V t 2 g f | 2 ∥t 2 g∥ 2 = A 2 4 √ 2 3σ(1 + (ϕ ′′ (t)σ 2 ) 2 ) 3/2 .
(5. [START_REF] Kodera | A new method for the numerical analysis of non-stationary signals[END_REF] This means that the bias is growing with the modulation. In the case of the signals we following first order Taylor expansion,

qf+ε = - 1 2π ℑ V g f + V g ε V t 2 g f + V t 2 g ε = - 1 2π ℑ        V g f V t 2 g f 1 1 + V t 2 g ε V t 2 g f + V g ε V t 2 g f 1 1 + V t 2 g ε V t 2 g f        ≈ - 1 2π ℑ V g f V t 2 g f (1 - V t 2 g ε V t 2 g f ) + V g ε V t 2 g f (1 - V t 2 g ε V t 2 g f ) ≈ qf + 1 2π ℑ    V g f V t 2 g ε (V t 2 g f ) 2 - V g ε V t 2 g f    .
(5.11)

In this last expression, as qf (t, η) is very close to ϕ ′′ (t), this means that the other term is responsible for the estimation error, that we detail hereafter,

ℑ    V g f V t 2 g ε (V t 2 g f ) 2 - V g ε V t 2 g f    = ℑ V g f V t 2 g f V t 2 g ε V t 2 g f - V g ε V g f = ℑ V g f V t 2 g f ℜ V t 2 g ε V t 2 g f - V g ε V g f + ℜ V g f V t 2 g f ℑ V t 2 g ε V t 2 g f - V g ε V g f = -2πϕ ′′ (t)ℜ V t 2 g ε V t 2 g f - V g ε V g f + 2π σ 2 ℑ V t 2 g ε V t 2 g f - V g ε V g f . (5.12)
We therefore define the error function by, G(t) := V t 2 g ε (t,ϕ ′ (t))

V t 2 g f (t,ϕ ′ (t)) -V g ε (t,ϕ ′ (t)) V g f (t,ϕ ′ (t))
. As we saw previously on Figure 5.1 (d), (e) and (f), G is oscillatory, we study its power spectral density (PSD) for the sake of removing these oscillations.

Power Spectral Density of the Error

Case of a Pure Harmonic Signal First, in the particular case where ϕ ′′ (t) is null, namely f (t) = e 2iπη 0 t , we remark that

V g f (t, ϕ ′ (t)) = f (t)
and rewrite:

G(t) = 2π σ 2 V t 2 g ε (t, ϕ ′ (t)) V g f (t, ϕ ′ (t)) - V g ε (t, ϕ ′ (t)) V g f (t, ϕ ′ (t)) = 2π σ 2 V t 2 g ε (t, ϕ ′ (t)) -V g ε (t, ϕ ′ (t)) e -2iπη 0 t , (5.13)
σ ε σ 3 2πη 2 c,0 e -πσ 2 η 2 c,0 = 2σσ ε ce -1 ⇔ η c,0 = -W(-ce -1 ) σ √ π .
(5.17)

In the case of a linear chirp, when ϕ ′′ = b and with the same reasoning, the maximum is

reached at η m = √ 1+b 2 σ 4 σ √ π . It verifies H G (η m ) = 2σ ε σe -1 . ( 5
σ ε σ 3 2πη 2 c,b (1 + b 2 σ 4 ) e - πσ 2 η 2 c,b 1+b 2 σ 4 = 2cσ ε σe -1 ⇔ η c,b = (1 + σ 4 b 2 ) 1/2 -W(-ce -1 ) σ √ π . (5.19)
However, this time, η c,b depends on the chirp rate b. Therefore, in the numerical simulations, we are going to set the low-pass filter at frequency η c,0 . As η c,0 ≤ η c,b , we expect to filter out more energy.

Evaluating the CR Filtering Process

In this section, we evaluate the quality of the filtered CR estimate that we denote by F ( q f ) with q f . We make this comparison in terms of output SNR with the true CR, namely SN R(ϕ ′′ , E -ϕ ′′ ), where E is any of the estimates. We also set c = 1 √ 10 resulting on a cut-off frequency at 10% of the maximum value. We study this choice of c for modulated and non-modulated modes in order to draw conclusions.

We first illustrate the effect of the filter in Figure 5.2 (d), (e) and (f). We remark that most of the oscillations are removed even though some remain but at low amplitude. Then, in Figure 5.2 (a), (b) and (c), we notice that up to a SNR of 40 dB, for the three signals considered, F ( q f ) offers significant improvement. Indeed, a close look at the box-plots tells us that the worst estimation of F ( q f ) is significantly better than the best estimate of q f , assuming that the RD is able to detect the ridge, which may not be the case for the signal of Figure 5.1 (c) at 0 dB. The purpose of Figure 5.3 is to show that the threshold η c,0 is still interesting to consider, even when the signal is modulated. On Figure 5.3 (a), we remark that the same cut-off frequency can be used, resulting only in a small loss of energy. On Figure

(b)

, for the most modulated signal, the loss of energy is even smaller.

Improving Linear Chirp Based Mode Retrieval

In this section, we investigate how to use the filtered CR estimator defined in Section 5.0.3 to improve LCR technique introduced in Section 2.6. We recall that this technique assumes a local linear chirp approximation for a mode, thus its STFT can be approximated by [START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF]:

V g f (t, η) ≈ V g f,approx (t, η) := A(t)e -πσ 2 1+iC(t)σ 2 1+(C(t)σ 2 ) 2 (η-B(t)) 2 , ( 5.20)
where A(t), B(t) and C(t) are respectively estimates of V g f (t, ϕ ′ (t)), ϕ ′ (t) and ϕ ′′ (t). In that context, mode reconstruction is performed by summing the approximation given by (5.20) over frequencies to obtain:

f (t) ≈ 1 g(0) R V g f,approx (t, η)dη, (5.21)
LCR technique is then based on a specific choice for A(t), B(t) and C(t). To compute the estimate B(t) of ϕ ′ (t), the technique uses the local instantaneous frequency estimator

ω [2] f = ℜ ω[2]
f that we introduced with FSST2 in section 1.2, with: where ωf and tf were introduced in Section 1.2. More precisely, denoting by R f (t) the STFT ridge corresponding to f , B(t) is set to ω

ω[2] f = ωf + ∂η ωf ∂η tf × (t -tf ) if ∂ η tf ̸ = 0 ωf otherwise, ( 5
(ϕ ′′ , ϕ ′ ) are (F ( q f +ε (t, R f +ε (t))), ω [2] f +ε (t, R f +ε (t))), (F ( q f +ε (t, R f +ε (t))), ϕ ′ ),
2] f (t, R f (t)), C(t) to q f (t, R f (t)), and finally A(t) to V g f (t, ω [2]
f (t, R f (t))). When one considers the noisy signal f + ε with n a complex Gaussian white noise, similar estimates are derived replacing f by f + ε in the definitions of A, B and C, which are thus all impacted by noise.

We now investigate the improvement brought by setting C(t) to F ( q f +ε (t, R f +ε (t))) instead of q f +ε (t, R f +ε (t)) in LCR technique, A(t) and B(t) remaining unchanged. We are going to compare this technique to original LCR and, to measure the impact of the estimator B on LCR technique, to implementations of LCR in which B(t) is set to ϕ ′ (t) and C(t) either to F ( q f +ε (t, R f +ε (t))) or to ϕ ′′ (t). To clarify this, as LCR can be viewed as a function of 3 parameters, we define:

LCR A,B,C (t, η) := A(t)e -πσ 2 1+iC(t)σ 2 1+(C(t)σ 2 ) 2 (η-B(t)) 2 . ( 5

.23)

In that framework, the original LCR technique corresponds to for 100 realizations of the noise, the range corresponding to the 5% and 95% quantiles of q f +ε evaluated on the STFT ridge corresponding to the fundamental frequency of the noisy signal whose STFT modulus is displayed in Figure 5.5, q f on the ridge corresponding to the fundamental frequency of the noiseless signal being also superimposed (input SNR = 5 dB for the noise); (b): same as (a) but with input SNR= 10 dB; (c): same as (a) but with input SNR= 15 dB. (d): same as (a) but using the filtered estimator F ( q f +ε ) instead of q f +ε , F ( q f ) is also superimposed ; (e): Same as (d) but with input SNR= 10 dB; (f): same as (d) but with input SNR= 15 dB.

(A, B, C) = (A 0 , B 0 , C 0 ) := V g f +ε (t, ω [2] f +ε (t, R f +ε (t))), ω [2] f +ε (t, R f +ε (t)), q f +ε (t, R f +ε (t)) , (5.24

SST Based Heart Rate Estimation

Introduction

The estimation of the HR is still an active topic of research and includes various approaches such as time based R-peak detection [START_REF] Pan | A real-time qrs detection algorithm[END_REF][START_REF] Kaplan | Simultaneous qrs detection and feature extraction using simple matched filter basis functions[END_REF] or TF analysis [START_REF] Huang | Ecg arrhythmia classification using stft-based spectrogram and convolutional neural network[END_REF][START_REF] Shaik | A method for qrs delineation based on stft using adaptive threshold[END_REF]. The latter is relevant for the study of physiological signals, as they can be highly non-stationary.

ECG signals are such cases, in which the shape of heart beat, namely the QRS complex, varies and can make time based analysis difficult. ECGs are actively studied using TF analysis, for instance, it is used to detect sleep apnea [START_REF] Quiceno-Manrique | Detection of obstructive sleep apnea in ecg recordings using time-frequency distributions and dynamic features[END_REF], to study arrhythmia [START_REF] Qurraie | Ecg arrhythmia classification using time frequency distribution techniques[END_REF], to identify coronary artery diseases [START_REF] Sharma | A new method to identify coronary artery disease with ecg signals and time-frequency concentrated antisymmetric biorthogonal wavelet filter bank[END_REF] or to separate fetal ECG from that of the mother on abdominal ECG signal of pregnant women [START_REF] Su | Extract fetal ecg from single-lead abdominal ecg by deshape short time fourier transform and nonlocal median[END_REF].

In this chapter, our goal is to estimate the HR of an ECG signal, which we assume to be harmonic. Thus, we are going to look for the fundamental frequency (FF) as it represents an estimate of the HR. To perform this estimation, some existing TF techniques try to estimate the spectrogram using non-negative matrix factorization, as for instance on PCGs [START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF] and ECGs [START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF] signals. In our case, we propose to define a dictionary matching both the choice of TFR and the harmonicity of the signal, which we only use to estimate the HR.

The layout of the chapter is as follows: first we define the harmonic model of a MCS and then use earth mover distance (EMD) to select a component in a dictionary for HR estimation. Finally, we discuss the influence of the choice of TFR on the quality of HR estimation and perform numerical assessment of the technique on synthetic and real ECG signals.

Sparsity of Harmonic Signals and Application to ECG Signals

A harmonic signal is a specific type of MCS defined by,

f (t) = P p=1
A p (t)e 2iπpϕ(t) . (6.1)

For such a signal, ϕ ′ (t) is the FF. In the case of an ECG signal, the instantaneous HR in beats per minute (bpm) at time t is given by 60ϕ ′ (t).

Dictionary Based Analysis of the Magnitude

We propose to study the magnitude of harmonic signals in a TFR X using a dictionary W containing K e components. A good choice of the number of components K e enables a good frequency resolution. In order to compare X with a component of W, we propose to use the EMD and to collect this information in a matrix H. We finally use the latter to select the best component and estimate the heart rate. The dictionary matrix W is fixed and describes possible structures of harmonicity depending on the variation of the FF. In the case of ECG signals, we consider a variation of the FF between 30 and 180 bpm and thus, with an accuracy of 1 bpm, set K e to 151.

In Figures 6.1 (a) and (b), we respectively display the STFT and FSST2 magnitude of one real ECG signal. We look for dictionaries allowing us to reasonably characterize these TFRs, and we assume FSST and FSST2 share the same dictionary. In the case of STFT, the dictionary takes into account the spreading of the energy in the TF plane.

We denote this dictionary by W 0 and give an illustration in Figure 6.1 (c), where the spreading of energy assumes that the analysis window is Gaussian. In the case of FSST and FSST2, the dictionary is built assuming perfect energy reassignment, meaning a single nonzero coefficient is enough to represent a harmonic. We denote this dictionary by W 1 and give an illustration in Figure 6.1 (d). Note that for these two dictionaries, a high component index takes more harmonics into account that for a low component index.

We assume that, as long that the highest component has enough harmonics and that the frequency range of the dictionary takes the sampling rate into account, it provides a good representation of the signal. We can confirm this in Figures 6.1 Our objective is to compare HR detection when the considered representations are STFT, FSST and FSST2. We do not take into account higher orders of FSSTN as it 6.2. SPARSITY OF HARMONIC SIGNALS AND APPLICATION TO ECG SIGNALS the minimal value of n med using

n med = min   n ≥ 1 s.t. n/2 k=0 C k n P k f (1 -P f ) n-k ≤ p 0   . (6.6)
The estimation of the HR at time indices 1, • • • , n med can then be done through Algorithm 2, where we look for the components in Î to avoid selecting a component related to noise or harmonics. Then, having defined the HR for the first n med indices, the estimation carries on with Algorithm 3 which is used to estimate HR after time instant n med . To make our detection more robust, we propose to take into account the standard deviation of previous detections, denoted by ∆ n , to define a range of relevant components. Such a range is defined using a parameter γ, the choice of which is discussed in Section 6.2.4.

Note that because Algorithm 3 considers only a subset of H, it does not require to compute the EMD for all components K e .

Algorithm 2: HR estimation step 1

Input : W,X, n med for n = 1, • • • , n med do ĩ[n] = argmin i∈[1,Ke] d(W i , X n ) end îM = median ĩ[n], n = 1, • • • , n med for n = 1, • • • , n med do î[n] = argmin i∈[ 3 4 îM , 3 2 îM ] d(W i , X n ) end Algorithm 3: HR estimation step 2 Input : W,X, n med , î[1, • • • , n med ], γ for n = n med + 1, • • • , L do ∆ n = std î[k], k = 1, • • • , n -1 î[n] = argmin [ î[n-1]-γ∆n, î[n-1]+γ∆n] d(W i , X n ) end HR Mean B STFT B FSST B FSST2 60

Determination of the Probability of False Detection

To compare the different probabilities of false detection depending on the type of TFR used to define X, we consider synthetic ECGs generated as in [START_REF] Mcsharry | A dynamical model for generating synthetic electrocardiogram signals[END_REF], compute the corresponding STFT, FSST, and FSST2 and then (6.3) for each time instant when X is associated with one of the just mentioned TFRs. With these synthetic signals, the mean HR component i M is known, and we measure good detection by considering I = [ 3i M 4 , 3i M 2 ] and then

B = #{n, argmin i d(W i , X n ) ∈ I} N , (6.7)
the probability of false detection being 1 -B.

We numerically notice that, with STFT or FSST2 the minimum of the EMD is much more likely to be associated with a sub-harmonic than when FSST is used, resulting in a much higher probability of false detection, and this is true whatever the mean HR ( the results are depicted in Table . 6.1).

An illustration of this is given in Figure 6.2 (a)-(c) when the TFR are either STFT, FSST and FSST2, respectively: the white curves correspond to the indices associated with the minimal EMD for each time index n. We notice that false detections usually occur when HR changes rapidly, and, as the signal is quasi-harmonic, to make a linear chirp approximation on the modes, as is done in FSST2, results in error at these time indices. As our plan is to use EMD for HR computation, we will make use of FSST rather than STFT or FSST2.

Study of Synthetic ECG

The goal of this section is to define an appropriate choice for parameter γ used in Algorithm 3. For that purpose, we investigate how good this algorithm is to estimate the along with the TF region corresponding to standard deviation ∆ n , when n varies.

Conclusion

In this chapter, we studied the particular case of harmonic signals in the context of ECG to estimate the HR. We first showed that the choice of TFR can have a huge impact on the accuracy of estimations, and that even though FSST does not take frequency modulation into account, it provides the best results. Based on a reliable HR estimator and a relevant choice of TFR, we designed an algorithm to track the HR in time. We evaluated this algorithm on synthetic and real ECG signals to show numerically the 7. Conclusion

Contributions

In the context of this thesis, we performed STFT based TF analysis of MCSs.

In Chapter 2, the contribution is two-fold: first, we designed a technique to identify ridges in the TF plane under very noisy conditions. The key point was to take into account that there is no continuous ridge associated with a mode, and to gather ridge portions to identify each mode. The second contribution was to use a linear chirp model to perform MR. Given the ridges of the signals, we locally approximate the linear chirp and recover the TFR for each mode.

In Chapter 3, we defined interference structures called TFBs and performed their complete mathematical characterization for the case of two parallel linear chirps. For practical purposes, we designed a technique to detect the TFBs, based on the detection of merging rigdes. Finally, we considered perspectives, using a different approach based on the spectrogram zeros.

In Chapter 4, we studied FSSTN and remarked that the selection of the order of this technique can be done by maximizing the concentration of energy. Then, since FSST2 performs well in noisy situations, we studied its associated CR estimator and observed that the errors due to noise are oscillatory. Therefore, we analyzed the power spectral density of the latter, to design a filtering technique, improving its accuracy.

In Chapter 6, we performed HR estimation on ECG signals, assuming that they can be modeled as a harmonic signal. We first compared the accuracy of a HR estimator based on the EMD, and demonstrated that the choice of the TFR can be critical. We finally designed an algorithm to detect the HR, taking into account the modulation of the signal.

Future research

To improve the state of the art, related to the techniques introduced in this thesis, we propose several directions we detail below. In what follows, we omit (t, η) in the definition of the STFTs!s (STFTs!s). Let us first remark that:

B. Towards Robust and Adaptive Synchrosqueezing

E{V g f +ε } = V g f , E{V t 2 g f +ε } = V t 2 g
f since V g ε and V t 2 g ε are both with zero mean. Let us denote by m = (V g f , V t 2 g f ) T . Then denoting by v = (V g f +ε , V t 2 g f +ε ) T , and σ 2 ε the variance of the noise (the real and imaginary part of the noise are assumed to be decorrelated), we have:

E{(v -m)(v -m) H } = E{|V g ε | 2 } E{V g ε (V t 2 g ε ) * } E{(V g ε ) * V t 2 g ε } E{|V t 2 g ε | 2 } = σ 2 ε ∥g∥ 2 σ 2 ε ∥tg∥ 2 σ 2 ε ∥tg∥ 2 σ 2 ε ∥t 2 g∥ 2 .
Based on the results in [START_REF] Li | On the ratio of two correlated complex gaussian random variables[END_REF]Eq. (44)], denoting Z =

V g f +ε V t 2 g f +ε
, we have that:

E {Z} = V g f V t 2 g f + ∥tg∥ 2 ∥t 2 g∥ 2 - V g f V t 2 g f e - |V t 2 g f | 2 σ 2 ε ∥t 2 g∥ 2 ,
so that one gets: Remarking that in that case ϕ ′ (t) = η 0 , we may write:

E - 1 2π ℑ {Z} = - 1 2π ℑ V g f V t 2 g f (1 -e - |V t 2
V t 2 g ε (t, η 0 )e -2iπη 0 t = R n(τ )(τ -t) 2 g(τ -t)e -2iπη 0 τ dτ V g ε (t, η 0 )e -2iπη 0 t = R n(τ )g(τ -t)e -2iπη 0 τ dτ, so we get that:

G(t) = 2π
σ 2 V t 2 g ε (t, η 0 )e -2iπη 0 t -V g ε (t, η 0 )e -2iπη 0 t = R n(τ )( 2π σ 2 (τ -t) 2 -1)g(τ -t)e -2iπη 0 τ dτ = σ 2 2π R n(τ )g ′′ (τ -t)e -2iπη 0 τ dτ.

Let us then write the auto-correlation function of this random variable:

E [G(t)G(t -x) * ] = R 2
E n(τ )n(τ ′ ) g ′′ (τ )g ′′ (τ ′ + x)

σ 4 e -2iπη 0 (τ -τ ′ ) 4π 2 dτ dτ ′ = σ 4 4π 2 σ 2 ε R g ′′ (τ )g ′′ (τ + x)dτ = σ 4 4π 2 σ 2 ε g ′′ * g ′′ (-x), using the fact that g ′′ is even. So the studied process is wide-sense stationary and its power spectral density, i.e. the Fourier transform of the auto-correlation function, reads: In such a case, we have ϕ ′ (t) = a + bt, and thus we may write:

F σ 4
V t 2 g ε (t, ϕ ′ (t))e -2iπϕ(t) = e iπbt 2 R ε(τ )(τ -t) 2 g(τ -t)e -2iπϕ ′ (t)τ dτ V g ε (t, ϕ ′ (t))e -2iπϕ(t) = e iπbt 2 R ε(τ )g(τ -t)e -2iπϕ ′ (t)τ dτ, Furthermore, from [19, proposition 1], we may write:

V g f (t, ϕ ′ (t)) = (1 + b 2 σ 4 ) -1 4 e -i atan(-bσ 2 ) 2 e 2iπϕ(t) = Be 2iπϕ(t) , so that one may rewrite G as:

G(t) = 2π σ 2 (1 -ibσ 2 )V t 2 g ε -V g ε B -1 e -2iπϕ(t) = B -1 e iπbt 2 R
ε(τ ) 2π σ 2 (1 -ibσ 2 )(τ -t) 2 -1 g(τ -t)e -2iπϕ ′ (t)τ dτ = B -1 e iπbt 2 σ 2 2π R ε(τ )g ′′ (τ -t)e -2iπϕ ′ (t)τ dτ -i2πb R ε(τ )(τ -t) 2 g(τ -t)e -2iπϕ ′ (t)τ dτ .

Then the auto-correlation function of G reads:

E [G(t)G(t -x) * ] = e -iπbx
L(τ ) = σ 2 2π g ′′ (τ ) -i2πbτ 2 g(τ ) e iπbτ 2
since one has:

E [G(t)G(t -x) * ] = σ 2 ε (1 + b 2 σ 4 ) 1 2 R L(τ - x 2 )L(τ + x 2 ) * dτ.
So, computing the Fourier transform of this auto-correlation, one gets: (1-ibσ 2 )

1 2 = -σ 3 2π
(1-ibσ 2 )

3 2
η 2 e -π( σ 2 1-ibσ 2 )η 2

.

From this, we finally deduce that the power spectral density of G 1 is:

P G 1 (η) = σ 2 ε σ 6 4π 2 η 4 (1 + b 2 σ 4 ) 2 e -2πσ 2 η 2 1+b 2 σ 4 .

XIII

Abstract

Time-frequency analysis is commonly used to study real world signals. These can often be described as multicomponent signals made of the sum of frequency and amplitude modulated modes. This thesis describes time-frequency techniques to study such signals in heavy noise situations, by dealing with three major issues. The first one is to design robust ridge detection technique and linear chirp approximation to improve instantaneous frequency estimation and mode reconstruction. The results demonstrate that an accurate ridge detection is a necessary but not a sufficient condition to ensure an accurate mode reconstruction. The second issue is the identification and the separation of interfering modes. To address this issue, the approach we propose focuses on ridge detection to localize patterns, coined time frequency bubbles, associated with interference in the timefrequency plane. The third issue is on the adaption of the synchrosqueezing transform to the frequency modulation of the modes and to noise. Regarding the first aspect, an energy based criteria is defined to measure the concentration of a time-frequency representation, which we use to adapt the synchrosqueezing technique. On the noise issue, based on a theoretical study of the effect of noise on the chirp rate estimator used in synchrosqueezing transforms, a new chirp rate denoising technique is proposed improving the estimation. The fourth and last issue is the heart rate estimation on ECG signals using time-frequency analysis, for which we design a specific algorithm. We show that the choice of the representation has huge consequences and demonstrate that it should be taken into account.

Résumé

L'analyse temps-fréquence est fréquemment utilisée pour étudier les signaux du monde réel. Ceux-ci peuvent souvent être décrits comme des signaux multicomposantes, représentés sous la forme d'une somme de modes, modulés en amplitude et en fréquence. Cette thèse traite de techniques temps-fréquence pour l'étude de tels signaux dans des situations où le niveau de bruit est élevé pour trois principales problématiques. La première est l'élaboration de détection de ligne de crêtes et l'approximation basées sur un modèle de chirp linéaire dans l'objectif de gagner en précision sur l'estimation de la fréquence instantanée et la reconstruction des modes. La deuxième problématique est l'identification et la séparation de modes qui interfèrent. Cette problématique est adressée, nous proposons une approche basée sur la détection de ligne de crête pour identifier des structures appelées bulles temps-fréquence, associées à l'interférence sur le plan temps-fréquence. La troisième problématique porte sur l'adaptation du synchrosqueezing à la modulation en fréquence des modes ainsi qu'au bruit. Sur ce premier aspect, nous définissons un critère basé sur l'énergie pour mesurer la concentration des représentations temps-fréquence, que nous utilisons pour adapter le synchrosqueezing. Concernant le problème du bruit, en nous basant sur l'étude théorique de l'effet du bruit sur l'estimateur de la modulation en fréquence, utilisée pour le synchrosqueezing, nous proposons une nouvelle technique de débruitage de manière à obtenir une meilleure estimation. Nous traitons une quatrième problématique, qui concerne l'estimation de la fréquence cardiaque en utilisant une approche temps-fréquence, pour laquelle nous construisons un algorithme. Nous montrons que le choix de la représentation à de grandes conséquences et qu'elle doit être prise en compte.

4 . 1 . 7 2. 1 .

 4171 Towards Adaptive and Robust Synchrosqueezing 55 4.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. Definition of FSSTN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. Adaptive selection of the order . . . . . . . . . . . . . . . . . . . . . . . . 4.4. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. Study of the Chirp Rate Estimate in Case of Noise 65 5.1. Evaluating the CR Filtering Process . . . . . . . . . . . . . . . . . . . . . 5.2. Frequency Variation of a Voice Signals . . . . . . . . . . . . . . . . . . . . 5.3. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .6. SST Based Heart Rate Estimation 79 6.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2. Sparsity of Harmonic Signals and Application to ECG Signals . . . . . . . 6.3. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7. Conclusion 89 7.1. Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2. Future research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Proof of Proposition 3.2.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . I A.2. Proof of Proposition 3.2.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . III B. Towards Robust and Adaptive Synchrosqueezing IX B.1. Proof of Proposition 5.0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . IX B.2. Proof of Proposition 5.0.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . X B.3. Proof of Proposition 5.0.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . XI vi ContentsList of AcronymsTFR time frequency representation TF time frequency FT Fourier transform STFT short-time Fourier transform RM reassignment method SST synchrosqueezing transform FSST Fourier-based synchrosqueezing transform FSST2 second order Fourier-based synchrosqueezing transform FSST4 fourth-order Fourier-based synchrosqueezing transform FSSTN high order Fourier-based synchrosqueezing transform MCS multicomponent signal IF instantaneous frequency CR chirp rate ITF ideal time frequency representation MR mode retrieval Acronyms of Chapter 2 RD ridge detection LMMF local maximum along the frequency axis S-RD simple ridge detector MB-RD modulation based ridge detector RRP relevant ridge portion RP ridge portion vii Contents RRP-RD relevant ridge portion based ridge detector SNR signal to noise ratio BA basin of attraction LCR linear chirp based retrieval HT hard thresholding SSR-HT shifted-symmetrized-regularized hard-thresholding S-LCR-MR MR based on S-RD and LCR RRP-LCR-MR MR based on RRP-RD and LCR RRP-MR MR based on RRP-RD and HT FSST4-MR MR based on FSST4 Acronyms of Chapter 3 LMF local maximum of the spectrogram along the frequency axis TFB time frequency bubble LEF local extremum of the spectrogram along the frequency axis R-RD redundant ridge detector R-RDN redundant ridge detector in noise GT ground truth PH pure harmonic Acronyms of Chapter 4 PSD power spectral density Acronyms of Chapter 6 HR heart rate viii Contents FF fundamental frequency ECG electrocardiogram NMF non-negative matrix factorization EMD earth mover distance bpm beats per minute ix List of Figures 1.1. (a) STFT of f (t) = exp(2iπ(50(t 4 + t 3 + t 2 ) + 100t)); (b) instantaneous chirp rate CR(t) = 600t 2 + 300t + 100 . . . . . . . . . . . . . . . . . . . . 5 1.2. Spectrogram of a bat echolocation signal. . . . . . . . . . . . . . . . . . . 6 1.3. (a): STFT of a MCS, (b) ITF of the signal of (a). . . . . . . . . . . . . . . LMMFs corresponding to the three largest STFT modulus maxima for each time instant along with the true IF (dashed line) of the linear chirp (SNR= -10 dB ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. (a): from top to bottom: STFT of a linear chirp, of a signal with cosine phase, and of a signal with oscillatory phase; (b): P(β) for these three signals when the input SNR equals -10 dB as well as the proportion of time indices n, at which the relevant LMMF with ordinate in I[n] corresponds to the largest STFT modulus maximum (curves denoted by "Max, linear chirp", "Max, cosine", and "Max, modulated cosine"). The results are averaged over 100 noise realizations. . . . . . . . . . . . . . . . . . . . . . 2.3. (a): STFT modulus of two noisy linear chirps (SNR = -10 dB, σ = 0.0188); (b): STFT modulus of two noisy modes with cosine phase, with different modulation (SNR = -10 dB, σ = 0.0175); (c): STFT modulus of a signal made of a linear chirp and a mode with exponential phase (SNR = -10 dB, σ = 0.0241); (d): (s f in p ) p=1,2 computed for the signal displayed in (a); (e): same as (d) but for the signal whose STFT modulus is displayed in (b); (f): same as (d) but for the signal whose STFT modulus is displayed in (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi List of Figures 2.4. (a): Comparison between S-RD, MB-RD and RRP-RD, for the mode f 1 of the signal of Figure 2.3 (a), computation of output SNR between IF ϕ ′ 1 and estimated IF with respect to input SNR (the results are averaged over 40 realizations of noise); (b): same as (a) but for the mode f 1 of the signal of Figure 2.3 (b); (c): same as (a) but for the mode f 1 of the signal of Figure 2.3 (c); (d): same as (a) but the mode f 2 of the signal of Figure 2.3 (a); (e): same as (b) but the mode f 2 of the signal of Figure 2.3 (b); (f): same as (c) but the mode f 2 of the signal of Figure 2.3 (c). . . . . . . 2.5. (a): STFT of f (t) = e 2iπ(250t+568 t 2 2 ) ; (b): STFT of f (t) = e 2iπ(1200t+60 cos(2πt)) ; (c): STFT of f (t) = e 2iπ(2000t+238 cos(2πt)) . . . . . . . . . . . . . . . . . . . 2.6. (a): Signal reconstruction results for the signal of Figure 2.5 (a) using HT, SSR-HT and LCR on the ridge provided by S-RD. the mode is reconstructed using either (1.7, page 3) or (1.8, page 3) ; (b): Same as (a) but for the signal depicted in Figure 2.5 (b); (c): Same as (a) but for the signal depicted in Figure 2.5 (c). . . . . . . . . . . . . . . . . . . . . . . . 2.7. (a): For each mode p = 1, 2, output SNR between mode f p of signal of Figure 2.3 (a) and reconstructed mode using either S-LCR-MR or RRP-LCR-MR (the results are averaged over 40 noise realizations); (b): same but with signal of Figure 2.3 (b); (c): same but with signal of Figure 2.3 (c); . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 .

 2 3 (a); (b): computation of the SNR associated with the estimation of ϕ ′′ with (s f in ) ′ (computed for different values of λ) or with q f for the first mode of Figure 2.3 (a). . . . . . . . . . . . . . . . . . . . . 2.9. (a): STFT of signal made of two close chirps; (b): Reconstruction results using (1.7, page 3) for each of the two modes depicted in (a) using S-RD.

2. 10 .Figure 3 . 4 ,

 1034 Figure 3.4, denoted by zone 1; (b): Interference mask associated with the linear chirp and the mode with sinusoidal phase of Figure 3.4, denoted by zone 2; (c): top graph: percentage of detection of TFBs in zone 1, zone 2, and elsewhere (false detections) when the noise level varies, bottom graph: percentage of points associated with the detection of TFBs in the different zones when the noise level varies (the number of noise realization equals 30). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7. (a): STFT of a voice signal (recording 62 of the Saarbruken Voice Database [65]) along with the ridges associated with 95% of the total energy (no TFB is detected); (b): same as (a) but when 98% of the total energy is located in the vicinity of the ridges (a TFB corresponding to the white star symbol is detected); (c): same as (a) but in the presence of noise (a TFB is detected is detect in the same TF region as in (b). . . . . . . . . . 3.8. Illustration of the zero based Delaunay triangulation, where the red circles represent the zeros of the spectrogram. . . . . . . . . . . . . . . . . . . . .

3. 9 .

 9 (a): Magnitude of the STFT of two PHs affected by noise, with superimposed Voronoi cells; (b): Zoom on (a) on a Voronoi cell associated with noise; (c): Zoom on (a) on a Voronoi cell associated with interference; (d): Zoom on (a) on a Voronoi cell associated with signal. . . . . . . . . . . . . 3.10. Histogram associated with the selection of the maximum value on ≈ 80 000 segments of Voronoi cells in the case of pure noise. . . . . . . . . . . . 4.1. (a): magnitude of the STFT of a mode with an oscillating phase; (b): E N for 2 ≤ N ≤ 4 for the signal whose STFT is depicted in (a) . . . . . . . . 4.2. (a): STFT magnitude of a mode with dumped sinusoidal frequency and amplitude A = 1; (b),(c),(d): E 2 ,E 3 and E 4 for the mode corresponding to (a), respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. (a): output SNR corresponding to SNR(f, R N f -f ) for N = 2, • • • , 4 or to SNR(f, R N f -f ) for the signal ofFigure 4.1 (a); (b): same as (a) but for the signal of Figure 4.2 (a). . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4. (a): modulus of the STFT of the sum of two modes with varying amplitude; (b): mode reconstruction results as a function of input SNR. . . . .

4. 5 .

 5 (a): STFT modulus of a bat echolocation call; (b): SNR corresponding to signal reconstruction from the coefficients on the ridges of the 3 most energetic detected modes. . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv List of Figures 5.1. (a): STFT modulus of a noisy linear chirp with the STFT ridge superimposed; (b): same as (a) but for a noisy chirp with fourth order polynomial phase; (c): same as (a) but for a noisy mode with oscillating phase; (d): CR estimator q f +ε or qf+ε estimated on the STFT ridge of the noisy signal whose STFT modulus is represented in (a), along with the true CR, for t between 0.3 and 0.5; (e): same as (d) but for the signal whose STFT modulus is represented in (b); (f): same as (d) but for the signal whose STFT modulus is represented in (c). The input SNR is 10 dB. . . . . . . . 68 5.2. (a): boxplot corresponding to the output SNR associated with CR estimators q f +ε or F ( q f +ε ) evaluated on the STFT ridges of the noisy signal corresponding to Figure 5.1 (a) over 100 realizations of the noise when the input SNR varies; (b): same as (a) but for the noisy signal corresponding to Figure 5.1 (b); (c): same as (a) but for the noisy signal corresponding to Figure 5.1 (c); (d): q f +ε and F ( q f +ε ) computed for a noisy version of the signal of Figure 5.1 (a), along with the ground truth ϕ ′′ (t) (input SNR = 10 dB); (e): same as (d) but for the signal of Figure 5.1 (b); (f): same as (d) but for the signal of Figure 5.1 (c). . . . . . . . . . . . . . . . . . . 72 5.3. (a): magnitude of the DFT of ϕ ′′ for the signals of Figure 5.1 (a) and

Figure 5 .

 5 Figure 5.1 (b); (b): same as (a) but for the signal of Figure 5.1 (c). . . . . 73 5.4. (a): percentage of improvement for the LCR technique computed using (5.27) for the signal of Figure 5.1 (a) where the estimates for (ϕ ′′ , ϕ ′ ) are

  or (ϕ ′′ , ϕ ′ );(b): same as (a) but for the signal of Figure 5.1 (b) ;(c): same as (a) but for the signal of Figure 5.1 (c) . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.5. ridge corresponding to the fundamental frequency superimposed on the STFT of the recording 62 taken from [65]. . . . . . . . . . . . . . . . . . . 76 5.6. (a): for 100 realizations of the noise, the range corresponding to the 5% and 95% quantiles of q f +ε evaluated on the STFT ridge corresponding to the fundamental frequency of the noisy signal whose STFT modulus is displayed in Figure 5.5, q f on the ridge corresponding to the fundamental frequency of the noiseless signal being also superimposed (input SNR = 5 dB for the noise); (b): same as (a) but with input SNR= 10 dB; (c):

  (e): Same as (d) but with input SNR= 10 dB; (f): same as (d) but with input SNR= 15 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 xv List of Figures 6.1. (a): STFT of an ECG generated following [92]; (b): FSST2 of the same ECG; (c): Dictionary matrix W 1 used with STFT-based approach; (d):

6. 2 .

 2 In these figures, the HR varies around 65 bpm, which corresponds to the index 35 of the dictionary. (a): modulus of STFT, the white curve corresponds to the indices associated with the minimal EMD ( B = 39.61%); (b):same as (a) but for FSST ( B = 94.34%) (c):same as (a) but for FSST2 ( B = 75.4%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3. (a): Estimation of HR mean from î. (b): Estimation of the corresponding standard deviation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4. (a) FSST based EMD of a real ECG signal. (b) HR detection on the same signal using γ = 2.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi List of Tables 6.1. Computation of B on syntethic ECG signals with different HR means when either STFT, FSST or FSST2 are used to compute X . . . . . . . . 84 xvii Notations ∥.∥ p L p norm, where p ∈ {1, 2}

Figure 1

 1 Figure 1.1.: (a) STFT of f (t) = exp(2iπ(50(t 4 + t 3 + t 2 ) + 100t)); (b) instantaneous chirp rate CR(t) = 600t 2 + 300t + 100

Figure 1 . 2 .

 12 Figure 1.2.: Spectrogram of a bat echolocation signal.

Figure 1

 1 Figure 1.3.: (a): STFT of a MCS, (b) ITF of the signal of (a).

13 CHAPTER 2 .Figure 2

 1322 Figure 2.1.: LMMFs corresponding to the three largest STFT modulus maxima for each time instant along with the true IF (dashed line) of the linear chirp (SNR= -10 dB )

Figure 2 .

 2 2 (a), which, from top to bottom, are more and more modulated, and get the results of Figure2.2 (b). These latter correspond to experiments with an input SNR of -10 dB, meaning that the noise is very strong. The quantity P(β) is almost constant for β between 0 and 2, before decreasing rapidly. We are therefore going to choose β = 2 assuming that there will be enough RRPs to represent the ridge of each mode. From now on, β equals 2 unless mentioned otherwise.

Figure 2 .

 2 Figure 2.2.: (a): from top to bottom: STFT of a linear chirp, of a signal with cosine phase, and of a signal with oscillatory phase; (b): P(β) for these three signals when the input SNR equals -10 dB as well as the proportion of time indices n, at which the relevant LMMF with ordinate in I[n] corresponds to the largest STFT modulus maximum (curves denoted by "Max, linear chirp", "Max, cosine", and "Max, modulated cosine"). The results are averaged over 100 noise realizations.

Figure 2 .Figure 2

 22 Figure 2.3.: (a): STFT modulus of two noisy linear chirps (SNR = -10 dB, σ = 0.0188); (b): STFT modulus of two noisy modes with cosine phase, with different modulation (SNR = -10 dB, σ = 0.0175); (c): STFT modulus of a signal made of a linear chirp and a mode with exponential phase (SNR = -10 dB, σ = 0.0241); (d): (s f in p ) p=1,2 computed for the signal displayed in (a); (e): same as (d) but for the signal whose STFT modulus is displayed in (b); (f): same as (d) but for the signal whose STFT modulus is displayed in (c).

Figure 2 . 3 .

 23 Figure 2.3. In particular, we here compare how close the ridge provided by any of the methods is close to the true IF for all the modes. For that purpose, we define the SNR by,

  Figure 2.3 (a) are the best in terms of SNR, as displayed in Figure 2.4 (a) and (d). Then for the most modulated signals of Figure 2.3 (c) and (d), the smoothing effect of λ becomes less desirable as the noise level diminishes.

Figure 2 .

 2 Figure 2.4.: (a): Comparison between S-RD, MB-RD and RRP-RD, for the mode f 1 of the signal of Figure 2.3 (a), computation of output SNR between IF ϕ ′ 1 and estimated IF with respect to input SNR (the results are averaged over 40 realizations of noise); (b): same as (a) but for the mode f 1 of the signal of Figure 2.3 (b); (c): same as (a) but for the mode f 1 of the signal of Figure 2.3 (c); (d): same as (a) but the mode f 2 of the signal of Figure 2.3 (a); (e): same as (b) but the mode f 2 of the signal of Figure 2.3 (b); (f): same as (c) but the mode f 2 of the signal of Figure 2.3 (c).

  Figure 2.5 (a), (b) and (c); the results being depicted in Figure 2.6. As the input SNR varies between 0 and 10 dB, we use S-RD as it is working fine at such noise levels.

28 2Figure 2 Figure 2

 2822 Figure 2.5.: (a): STFT of f (t) = e 2iπ(250t+568 t 2 2 ) ; (b): STFT of f (t) = e 2iπ(1200t+60 cos(2πt)) ; (c): STFT of f (t) = e 2iπ(2000t+238 cos(2πt)) .

Figure 2 . 7 .Figure 2 . 8 .

 2728 Figure 2.7.: (a): For each mode p = 1, 2, output SNR between mode f p of signal of Figure 2.3 (a) and reconstructed mode using either S-LCR-MR or RRP-LCR-MR (the results are averaged over 40 noise realizations); (b): same but with signal of Figure 2.3 (b); (c): same but with signal of Figure 2.3 (c);

Figure 2 . 9 .

 29 Figure 2.9.: (a): STFT of signal made of two close chirps; (b): Reconstruction results using (1.7, page 3) for each of the two modes depicted in (a) using S-RD.

Figure 2 .

 2 Figure 2.10.: (a): STFT modulus (σ = 0.0105) of the Hanford signal along with the ridge obtained with RRP-RD and FSST4; (b): illustration of signal reconstruction based on RRP-LCR-MR and the numerical relativity; (c): SNR corresponding to the reconstruction of the signal using either RRP-LCR-MR, RRP-MR or FSST4-MR (the ground truth is assumed to be the one produced by numerical relativity). The results are averaged over 40 realizations of noise.

  RD and HT (RRP-MR), RRP-LCR-MR or MR based on FSST4 (FSST4-MR). An illustration of the reconstructed signal obtained with RRP-LCR-MR superimposed on the one given by the numerical relativity is displayed in Figure2.10 (b). Then, we 2.10. CONCLUSION estimate the quality of signal reconstruction with the three techniques when the input SNR varies. The results depicted in Figure2.10 (c) show that the three methods behave similarly, though RRP-MR and RRP-LCR-MR are always slightly better. What is very specific to the studied gravitational wave signal is that the part associated with the strong frequency modulation is very energetic, and a slight inaccuracy in IF estimation using the spline approximation at the locations corresponding to strong frequency modulation has a strong impact on mode reconstruction with RRP-LCR-MR. For this reason, it

. 2 )

 2 Denoting the local extrema of the spectrogram along the frequency axis by local extrema of the spectrogram along the frequency axis (LEFs), the study of |V g f (t k , .)| 2 leads to: Proposition 3.2.1. |V g f (t k , .)| 2 has three LEFs (two LMFs and a local minimum) if and only if α := π 2 σ(ξ 2 -ξ 1 ) > 1 and

  Figure 3.1 (c), A ̸ = 1 and the conditions of Proposition 3.2.2 are fulfilled. In each case, we also plot the zeros and the local maxima of the spectrogram along with the spectrogram ridges corresponding to the chaining along the time axis of LMFs. As expected, the local maxima are located on spectrogram ridges at t = t k .

r 2 σ 2 cosFigure 3

 23 Figure 3.1.: (a): spectrogram of two parallel pure harmonic modes with the same amplitude, when there exists three LEFs at each time instant (the ridges corresponding to the chaining of LMFs are also superimposed); (b): spectrogram of two parallel pure harmonic modes exhibiting TFBs; (c): same as (a) but when the pure harmonic modes have different amplitude (three LEFs are still present at each time instant but the symmetry no longer exists). In each case, we also plot the zeros and the local maxima of the spectrogram.

a 2 +CB 1 k 1 - 2 a 2 -a 1 rσ > 1 ,

 112211 CD , and thus one can deduce the following Proposition 3.2.3. |V g f (t, η)| 2 has three extrema (two maxima and one minimum) along 3.2. ON THE EXISTENCE OF TIME FREQUENCY BUBBLES the line t = Dη + B 1 k if and only if α C = π and then

2 . 1 . 2 a 2 -a 1 rσ > 1 ,

 212211 Then one generalizes Proposition 3.2.2 as follows: Proposition 3.2.4. |V g f (t, η)| 2 has three extrema (two maxima and one minimum) along the line t = Dη + B λ k if and only if α C = π and then

Figure 3

 3 Figure 3.2.: (a): Spectrogram of two parallel linear chirps with the same amplitude, when there always exists three extrema in the direction D; (b): Spectrogram of two parallel linear chirps when TFBs are present. Three extrema no longer exist in the direction D; (c): same as (a) but when the linear chirps have different amplitude (three extrema are still present along the direction D). In (a) to (c), we also plot the zeros and the local maxima of the spectrogram.

40 3. 3 .

 403 LOCALIZING TFBS BASED ON A NEW RD initialization at time index n 0 ,

Figure 3 .

 3 Figure 3.3 (a) and (b), corresponding to the detection of two ridges with S-RD and to the application of R-RD respectively. As explained above, S-RD enforces the detection of two separate ridges that are not made of LMFs, contrary to R-RD with which the extracted ridge portions can be on top of one another and are necessarily made of LMFs.

Figure 3

 3 Figure 3.3.: (a): ridges detected using S-RD; (b): ridges detected using R-RD.

CHAPTER 3 .Figure 3 . 4 .

 334 Figure 3.4.: (a): STFT of a three mode signal along with the ridges superimposed, in the definition of the Gaussian window σ = 0.02, the stars corresponds to the locations of the merging of ridges (R-RD is used for ridge detection); (b): same as (a) but with σ = 0.03 ; (c): same as (a) but with σ = 0.04 ;(d): study of the same signal as in (a) with some white Gaussian added noise (input SNR 20 dB), and σ = 0.04 in the definition of the Gaussian window; (e): same as (d) except the input SNR equals 10 dB;(f): same as (d) except the input SNR equals 5 dB.

Figure 3 . 2 - 1 .

 321 1 (c), where the modes are with different amplitude, GT has to be modified because TFBs are present either if α ≤ 1 or if α > 1 and | log(A)| ≥ -2 arccosh(α) + 2α √ α In Figure 3.5 (b), we display GT, as well as

Figure 3

 3 Figure 3.5.: (a): TFB detection for the signal of Figure 3.1 (a) either in the noiseless case or for different noise levels, using the function R, along with the ground truth; (b): TFB detection for the signal of Figure 3.1 (c) using either R or the ground truth. The results are averaged over 10 realizations of the noise.

. 10 )Figure 3 . 6 .

 1036 Figure 3.6.: (a): Interference mask associated with the constant and linear chirps of Figure 3.4, denoted by zone 1; (b): Interference mask associated with the linear chirp and the mode with sinusoidal phase of Figure 3.4, denoted by zone 2; (c): top graph: percentage of detection of TFBs in zone 1, zone 2, and elsewhere (false detections) when the noise level varies, bottom graph: percentage of points associated with the detection of TFBs in the different zones when the noise level varies (the number of noise realization equals 30).

Figure 3 .Figure 3

 33 Figure 3.7.: (a): STFT of a voice signal (recording 62 of the Saarbruken Voice Database [65]) along with the ridges associated with 95% of the total energy (no TFB is detected); (b): same as (a) but when 98% of the total energy is located in the vicinity of the ridges (a TFB corresponding to the white star symbol is detected); (c): same as (a) but in the presence of noise (a TFB is detected is detect in the same TF region as in (b).

Figure 3 . 8 .

 38 Figure 3.8.: Illustration of the zero based Delaunay triangulation, where the red circles represent the zeros of the spectrogram.

  Figures 3.9 (b), (c) and (d) respectively illustrate Voronoi cells associated with noise, interference and the signal. In these figures, one can observe that the energy in a segment can be significantly larger when it crosses signal information.

Figure 3

 3 Figure 3.9.: (a): Magnitude of the STFT of two PHs affected by noise, with superimposed Voronoi cells; (b): Zoom on (a) on a Voronoi cell associated with noise; (c): Zoom on (a) on a Voronoi cell associated with interference; (d): Zoom on (a) on a Voronoi cell associated with signal.

Figure 3 .

 3 Figure 3.10.: Histogram associated with the selection of the maximum value on ≈ 80 000 segments of Voronoi cells in the case of pure noise.

  From 4.1 (b),we remark that E 4 is stable around one and is the best estimate for this signal. In Fig

Figure 4 .E 2 ,

 42 Figure 4.1.: (a): magnitude of the STFT of a mode with an oscillating phase; (b): E N for 2 ≤ N ≤ 4 for the signal whose STFT is depicted in (a)

Figure 4 . 2 .

 42 Figure 4.2.: (a): STFT magnitude of a mode with dumped sinusoidal frequency and amplitude A = 1; (b),(c),(d): E 2 ,E 3 and E 4 for the mode corresponding to (a), respectively.

Figure 4 .

 4 4 (a), which contains two modes where, this time, the amplitude varies. The output SNR for the different techniques is depicted in Figure4.4 (b). As expected, using R N p fp leads to a higher output SNR, which

Figure 4 . 3 .

 43 Figure 4.3.: (a): output SNR corresponding to SNR(f, R N f -f ) for N = 2, • • • , 4 or to SNR(f, R N f -f ) for the signal ofFigure 4.1 (a); (b): same as (a) but for the signal of Figure 4.2 (a).

Figure 4 .

 4 Figure 4.4.: (a): modulus of the STFT of the sum of two modes with varying amplitude; (b): mode reconstruction results as a function of input SNR.

  fp , N a meaning we consider an adaptive N in the MCS case. Similarly to R Na f , we respectively denote signal reconstruction from the three detected modes for, 1 ≤ N ≤ 4, by R N f . We then compare R Na f and R N f , for different N , to f in terms of SNR, and we again notice that the first technique is much better than the others (see the Table of

  Figure 4.5 (b)).

Figure 4 .

 4 Figure 4.5.: (a): STFT modulus of a bat echolocation call; (b): SNR corresponding to signal reconstruction from the coefficients on the ridges of the 3 most energetic detected modes.

|Figure 5 . 1 .

 51 Figure 5.1.: (a): STFT modulus of a noisy linear chirp with the STFT ridge superimposed; (b): same as (a) but for a noisy chirp with fourth order polynomial phase; (c): same as (a) but for a noisy mode with oscillating phase; (d): CR estimator q f +ε or qf+ε estimated on the STFT ridge of the noisy signal whose STFT modulus is represented in (a), along with the true CR, for t between 0.3 and 0.5; (e): same as (d) but for the signal whose STFT modulus is represented in (b); (f): same as (d) but for the signal whose STFT modulus is represented in (c). The input SNR is 10 dB.

. 18 )

 18 Then, similarly to the PH case, we set η c,b at which H G (η c,b ) = cH G (η m ) and solve its corresponding equation,

Figure 5 . 2 .

 52 Figure 5.2.: (a): boxplot corresponding to the output SNR associated with CR estimators q f +ε or F ( q f +ε ) evaluated on the STFT ridges of the noisy signal corresponding to Figure 5.1 (a) over 100 realizations of the noise when the input SNR varies; (b): same as (a) but for the noisy signal corresponding to Figure 5.1 (b); (c): same as (a) but for the noisy signal corresponding to Figure 5.1 (c); (d): q f +ε and F ( q f +ε ) computed for a noisy version of the signal of Figure 5.1 (a), along with the ground truth ϕ ′′ (t) (input SNR = 10 dB); (e): same as (d) but for the signal of Figure 5.1 (b); (f): same as (d) but for the signal of Figure 5.1 (c).

Figure 5 . 3 .

 53 Figure 5.3.: (a): magnitude of the DFT of ϕ ′′ for the signals of Figure 5.1 (a) and Figure 5.1 (b); (b): same as (a) but for the signal of Figure 5.1 (c).

Figure 5 . 4 .

 54 Figure 5.4.: (a):percentage of improvement for the LCR technique computed using (5.27) for the signal of Figure5.1 (a) where the estimates for (ϕ ′′ , ϕ ′ ) are (F ( q f +ε (t, R f +ε (t))), ω

  or (ϕ ′′ , ϕ ′ );(b): same as (a) but for the signal of Figure 5.1 (b) ;(c): same as (a) but for the signal of Figure 5.1 (c)

[

Figure 5 . 6 .

 56 Figure 5.6.: (a): for 100 realizations of the noise, the range corresponding to the 5% and 95% quantiles of q f +ε evaluated on the STFT ridge corresponding to the fundamental frequency of the noisy signal whose STFT modulus is displayed in Figure5.5, q f on the ridge corresponding to the fundamental frequency of the noiseless signal being also superimposed (input SNR = 5 dB for the noise); (b): same as (a) but with input SNR= 10 dB; (c): same as (a) but with input SNR= 15 dB. (d): same as (a) but using the filtered estimator F ( q f +ε ) instead of q f +ε , F ( q f ) is also superimposed ; (e): Same as (d) but with input SNR= 10 dB; (f): same as (d) but with input SNR= 15 dB.

  (a) and (b), where most of the energy is located in the lower part of the TFR.

Figure 6 . 3 .

 63 Figure 6.3.: (a): Estimation of HR mean from î. (b): Estimation of the corresponding standard deviation.

Figure 6

 6 Figure 6.4.: (a) FSST based EMD of a real ECG signal. (b) HR detection on the same signal using γ = 2.5.

B. 1 .

 1 Proof of Proposition 5.0.1

4π 2 σ 2 εσ 2 ε σ 4 4π 2 Fσ 2 ε σ 4 4π 2 g 3 B. 3 .

 224233 g ′′ * g ′′ (-x) (η) = g ′′ * g ′′ (η) * = ′′ (η) 2 = σ 2 ε σ 4 4π 2 η 4 g(η) 2 = σ 2 ε σ 6 4π 2 η 4 e -2πσ2 η 2 . X B.3. PROOF OF PROPOSITION 5.0.Proof of Proposition 5.0.3

2 σ 2 ε( 1 + b 2 σ 4 ) 1 2 σ 4 4π 2 Rg 2 R τ 2 2 R2 b 2 Rτ 2 2 ε( 1 +

 2122222221 ′′ (τ )g ′′ (τ + x)e -2iπbxτ dτ ibσ g(τ )g ′′ (τ + x)e -2iπbxτ dτ -ibσ(τ + x) 2 g ′′ (τ )g(τ + x)e -2iπbxτ dτ 4π g(τ )(τ + x) 2 g(τ + x)e -2iπbxτ dτ .XI APPENDIX B. TOWARDS ROBUST AND ADAPTIVE SYNCHROSQUEEZINGSo the process is wide-sense stationary. Now, making a last change of variables, we obtain:E [G(t)G(t -x) * ] = σ b 2 σ 4 )

RE 2 ε( 1 + b 2 σ 4 ) 2 ε( 1 + b 2 σ 4 ) 1 2 RW 2 = 2 + 1 σ 2 +ib)τ 2

 21421222122 [G(t)G(t -x) * ] e -2iπxη dx = σ L (τ, η)dτ, in which W L is the Wigner-Ville distribution of L. From the properties of the Wigner-Ville distribution, we deduce that the modulus squared of the Fourier transform of L reads: | L(η)| 2 = R W L (τ, η)dτ. XII B.3. PROOF OF PROPOSITION 5.0.3Now as we use a Gaussian window we can have an analytic expression for L. Indeed, some simple computations lead to:L(τ ) = (( 2π σ 2 τ 2 -1) -i2πbτ 2 )g(τ )e iπbτ ib)τ 2 -1)e π(-So, if one considers the Fourier transform of this function, one gets:

  ′′ Based on FSST2 and RidgeDetectionThe objective of this section is to explain how the IF and CR estimators involved in the definition of FSST2, introduced in Section 1.2, can be used in the model introduced 2.6. LINEAR CHIRP BASED MODE RETRIEVAL TECHNIQUE above. In fact, to get an acceptable level of accuracy, those estimates are to be evaluated in TF regions associated with high energy called ridges. Therefore, denoting the ridge by Γ p , setting ψ ′

.26)

This last equation gives us a way to approximate each mode in the general case of a MCS. However, ϕ ′ p (t), and ϕ ′′ p (t) are not known, and we are now going to explain how the FSST2

[START_REF] Behera | Theoretical analysis of the second-order synchrosqueezing transform[END_REF]

provides reasonable values for ϕ ′ p (t) and ϕ ′′ p (t).

2.5. Approximation of ϕ ′ and ϕ p (t) = ω

  from values of the discrete STFT is needed. For that purpose, let us consider k 0 := ⌊ψ ′

	CHAPTER 2. MODE RETRIEVAL BASED ON ROBUST RIDGE DETECTION
	AND LINEAR CHIRP APPROXIMATION
	to
	p [n] N L ⌉, where ⌊x⌉ denotes the closest integer

  .38) CHAPTER 2. MODE RETRIEVAL BASED ON ROBUST RIDGE DETECTION AND LINEAR CHIRP APPROXIMATION Here again, similarly to LCR and HT, we respectively denote the two reconstruction techniques by SSR M 1 and SSR M 2 .

  1. Having studied the situation at t k , we can generalize this result to any time t through:

	Proposition 3.2.2. |V g f (t, .)| 2 has three LEFs (two LMFs and a local minimum) if and
	only if α

Table 6 .

 6 1.: Computation of B on syntethic ECG signals with different HR means when either STFT, FSST or FSST2 are used to compute X

  Now, one can check that the auto-correlation function can be rewritten in terms of the function:

												1
												2
						σ 4 4π 2	R	g ′′ (τ -	x 2	)g ′′ (τ +	x 2	)e -2iπbxτ dτ
			+ibσ 2	R	(τ -	x 2	) 2 g(τ -	x 2	)g ′′ (τ +	x 2	)e -2iπbxτ dτ
			-ibσ 2	R	(τ +	x 2	) 2 g ′′ (τ -	x 2	)g(τ +	x 2	)e -2iπbxτ dτ
	+4π 2 b 2	R	(τ -	x 2	) 2 g(τ -	x 2	)(τ +	x 2	) 2 g(τ +	x 2	)e -2iπbxτ dτ .

This chapter is based on the work published in[START_REF] Laurent | A novel time-frequency technique for mode retrieval based on linear chirp approximation[END_REF] and[START_REF] Laurent | A novel ridge detector for nonstationary multicomponent signals: Development and application to robust mode retrieval[END_REF].

This chapter is based on the work published in[START_REF] Laurent | On local chirp rate estimation in noisy multicomponent signals: With an application to mode reconstruction[END_REF].

√ πσ . To solve such an equation, we make use of the Lambert W function, denoted by W,

This chapter is based on the work published in[START_REF] Laurent | A novel algorithm for heart rate estimation based on synchrosqueezing transform[END_REF].

6.3. CONCLUSIONrelevance of the proposed approach.

Acknowledgment

I would first like to express my warmest thanks to my supervisor, Sylvain Meignen, who was inspiring, dedicated and benevolent both as a teacher and as a supervisor. I would also like to thank my co-supervisors, Bertrand Rivet and Julie Fontecave-Jallon, for their availability, kindness and for the helpful discussions.

I thank Maria Sandsten and Roland Badeau, for their time reporting my thesis. I further thank them with the other members of the jury; Patrick Flandrin, Jérôme Mars and Pierre Chainais for rich interaction and their precious advices.

I thank my colleagues in the Laboratoire Jean Kuntzmann and my collaborators, for the very rich and friendly discussions. I also thank them for the advices they gave me for teachings. I also thank those who encouraged me in my studies, including some of my previous teachers and dear friends. Finally, I thank my family for all the trust and freedom they gave me.

iii

to get the following: Proposition 5.0.2. The power spectral density of G(t) is

The proof is available in Appendix B.2.

Case of a Linear Chirp

Now let us consider the more general case of a linear chirp f (t) = e 2iπ(at+ b 2 t 2 ) = e 2iπϕ(t) , for which we may write on the STFT ridge, using (5.8),

and thus:

leading to the following:

Proposition 5.0.3. The power spectral density of G(t) is

The proof is available in Appendix B.3.

Filtering the Oscillatory Error

Using the PSD computed above, we justify the use of a low-pass filter on q f . We first look for the maximum value of P G , by focusing on H G := √ P G , as it has its maximum at the same frequency and has a simpler expression. In the case of a PH, when ϕ ′′ = 0, we remark that

(5.16)

For the purpose of defining a low-pass filter, we first define a cutoff value that we denote η c,0 . We propose to choose c ∈ ]0, 1], and set the respective fraction of the maximum value such that H G (η c,0 ) = cH G

FREQUENCY VARIATION OF A VOICE SIGNALS the proposed technique to:

(A, B, C) = (A 0 , B 0 , F ( q f +ε (t, R f +ε (t)))), (5.25) while the other two references to (A, B, C) = (A 0 , ϕ ′ (t), F ( q f +ε (t, R f +ε (t))))

(A, B, C) = (A 0 , ϕ ′ (t), ϕ ′′ (t)).

(5.26)

For a monocomponent signal f , we define f (A,B,C) the reconstructed signal using LCR technique with the set of parameters (A, B, C). Then, using the definition of SN R out introduced in (2.21, page 22), we compute the percentage of improvement brought by using a set (A, B, C) different from (A 0 , B 0 , C 0 ) as follows:

(5.27)

The results are depicted in Figure 5.4 for the three signals of Figure 5.1 and when (A, B, C) are either defined by (5.25) or (5.26). First, these clearly show the improvement brought by using

Second, these results also tell us that to know the actual ϕ ′ (t) does not result in an extra improvement when compared with the proposed method. Finally, we notice that to know the actual ϕ ′′ enables a better mode reconstruction, but only at high SNRs and only for the first two signals. On the contrary, the gain of knowing the true IF and CR does not result in a significant improvement in terms of reconstruction for the last signal, except at very low SNR, but then ridge detection may fail, making the interpretation of the results more difficult.

So, we reach the conclusion that the proposed estimator

f +ε (t, R f +ε (t)) for (ϕ ′′ , ϕ ′ ) leads to nearly optimal reconstruction. If one would like to further improve mode reconstruction with LCR technique one should investigate how to find a better estimate of V g f (t, ϕ ′ (t)) than A 0 (t), but this is beyond the scope of this thesis.

Frequency Variation of a Voice Signals

In this section, the objective is to study the stability of the estimate F ( q f ) which we apply to the same voice signal as in Section 3.4.3, namely the recording 62 of the Saarbruken

Voice Database [START_REF] Woldert-Jokisz | Saarbruecken voice database[END_REF]. The signal depicted in Figure 5.5 shows its TFR with superimposed ridge. We recall that it contains sustained /a/ vowels of the type "low-high-low", and we want our CR estimate to preserve the change in pitch. On figure 5.6, we compare the stability of F ( q f ) and q f over 100 realizations of noise. We test this for three noise levels, in 5.6 (a), (d) for 5dB, 5.6 (b), (e) for 10dB and 5.6 (c), (f) for 15dB. In all cases, the 5% and 95% quantiles associated with the estimation of F ( q f ) are much closer to the estimation without noise and still preserve the variation of the pitch that we remark on Figure 5.5 around 0.75s and 1.65s. We may therefore conclude that F ( q f ) seems to be improved in terms of stability with respect to q f .

Conclusion

In this chapter, we introduced the filtering of the CR estimator associated with FSST2 to remove most of the oscillatory errors due to noise. We carried out the theoretical study of the estimator in both pure harmonic and linear chirp cases. We also evaluated the filtered estimate on numerical experiments, first on synthetic signals, where the filtered estimator gives significant improvement compared to the original one, in terms of output SNR. In a second part, we verified on a real voice signal that it is also more robust, meaning that the noise has a smaller effect on the estimation. From these results, we conclude that the filtered estimate is significantly more robust to noise. SIGNALS does not increase the accuracy of the estimation.

In order to make a decomposition, one needs to set the cost function D. In the case where D is the Euclidian distance, the decomposition is a least-square minimization problem,

and thus H satisfies W T WH = W T X, and is not necessarily unique because W T W is not invertible. Since the FF is associated with only one component of W, we impose stronger conditions to obtain the unicity of the solution. At each time instant, we compute the EMD [START_REF] Peleg | Fast and robust earth mover's distances[END_REF] between the i-th column of W and the n-th column of X to obtain a matrix H defined by,

where d is the EMD.

In this setting, finding a minimum in H is equivalent to finding the most relevant component to associate with the signal at a time index, and to estimate the associated FF. We denote this component by i 0 which satisfies,

In the following section, we propose an algorithm to estimate the HR using i 0 for which we take into account the accuracy induced by the choice of the representation for X.

Algorithm for HR Computation

We consider the HR estimate as the minimum on i of argmin

To assess the quality of this estimate, we set

and consider it as an estimate of the mean HR. We now seek to evaluate the quality of the HR estimate given by the minimum of the EMD. For that purpose, we define the following range of frequencies, Î := [ 3 îM 4 , 3 îM 2 ] that we consider as valid, and then define B as follows:

The choice for Î is supported by the fact that, assuming that the HR varies slowly, the upper harmonic admits 3 îM 2 as a lower bound and that no minimum should be detected under the lower bound of the FF, namely 3 îM 4 . Based on this model, we may define the probability of false detection as

However, in order to trust the median, we have to make sure that P f < 1 2 and that the signal is sufficiently long. For that purpose, imposing that the median of the detection computed from time indices 1, • • • , n med is exact with probability p 0 or above, we compute SIGNALS

Study of Real ECGs

We here study ECGs extracted from SiSEC database [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF], which correspond to thoracic recordings, which can be contaminated by different types of noises.Our HR estimation computed on such a signal, using γ = 2.5, fixed as explained in Section 6.2.4, is depicted in Figure 6.4 (a).

As FSST somehow reduces the noise while enhances the significant harmonics, the proposed HR estimation algorithm never gets trapped by minima that arise at subharmonic locations. We display HR estimations given by î in Figure 6. The IF estimator associated with FSSTN is sensitive to noise and its improvement in terms of accuracy would make FSSTN more stable. In Chapter 4, we made progress towards such an estimate in the case of FSST2, by remarking that the nature of the error associated with the CR estimator is oscillatory, enabling the design of a denoising technique. Further research on high order IF estimators may lead to the design of an adaptive technique in noisy case, which should match the order of FSSTN with the phase for each mode of the signal.

Concerning the particular case of ECG signals, more knowledge on interference and TF localization can also help to identify the HR. Particularly, the case of measuring the fetal HR from a single abdominal ECG channel is a challenging task, since mother and fetal heart beats can interfere. Moreover, it also has the added difficulty that the fetal signal have less power and its HR is therefore more difficult to read.

Appendices

A. Time Frequency Analysis of Interference

A.1. Proof of Proposition 3.2.1

Let us consider the function (Ae -πσ 2 (η-ξ 1 ) 2 + e -πσ 2 (η-ξ 2 ) 2 ) 2 having the same extrema as

and making the change of variables η = ν + ξ 1 +ξ 2 2 , one gets:

which we denote by q(ν). Differentiating q, one obtains:

So q ′ (ν) has the sign of h(ν) := -ν ξ + tanh(2πσ 2 ξν -log(A) 2 ) whose derivative is:

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, h ′ has two zeros, and one has the following table of variations:

If h(ν 1 ) ≥ 0, then h(ν) and hence q ′ (ν) annihilates once and changes signs at some ν ′

and strictly decreasing on [ν ′ 2 , +∞[, and thus has a maximum at ν = ν ′ 2 which is the unique extremum.

If h(ν 2 ) ≤ 0, h(ν) and hence q ′ (ν) annihilates and changes signs once for a certain ν ′ 1 in ] -∞, ν 1 [ and then q has a maximum at ν = ν ′ 1 which is its unique extremum.

If h(ν 1 ) < 0 and h(ν 2 ) > 0, h(ν) 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.

To determine when q has three extrema we need to compute h(ν 1 ) and h(ν 2 ). Remember that ν 1 and ν 2 are the roots of h ′ thus also of

and therefore

and thus as α = π 2 σ(ξ 2 -ξ 1 ) = √ 2πσξ, we may write:

From this we deduce that h(ν 1 ) < 0 and h(ν 2 ) > 0, and thus q has three extrema, if and only if:

A.

In any other case, |V g f (t, .)| 2 has a unique extremum.

Proof. Let us consider the function: 2 , one may define:

=e -2πσ 2 (ν 2 +ξ 2 ) (A 2 e -4πσ 2 ξν + e 4πσ 2 ξν + 2Aγ).

Putting ν = log(A)+µ 4πσ 2 ξ , enables us to define:

2A(cosh(µ) + γ).

The derivative of f 2 reads:

which has the sign of (assuming γ ̸ = -1):

Differentiating h we get

which has the same sign as

).

The term inside the parentheses is a second order polynomial in cosh(µ) whose discrim-

The roots of this polynomial are denoted by

To know the location of cosh(µ) with respect to X 1 and X 2 , one first compute Y (1) = (γ + 1)(1 -γ+1 2α 2 ) which has the same sign as

1 belongs to ]X 2 , +∞[. Finally, as cosh(µ) ≥ 1, Y (cosh(µ)) and thus h ′ (µ) are negative.

From this study one deduces the following table of variations:

and l(η) has a unique extremum which is a maximum when α ≤ γ+1 2 .

If one now assumes that α > γ+1 2 , then Y (1) > 0 and thus

] arccosh(X 2 ), +∞[, and Y (cosh(µ)) is positive otherwise, leading to the following table of variations:

vanishes with a change of sign only once at some µ = µ 1 if and only if h(-arccosh(X 2 )) ≥

V APPENDIX A. TIME FREQUENCY ANALYSIS OF INTERFERENCE

0 or h(arccosh(X 2 )) ≤ 0. In this case, we deduce that

meaning l(η) has a unique extremum (which is a maximum).

If h(-arccosh(X 2 )) < 0 and h(arccosh(X 2 )) > 0, h(µ) vanishes and changes signs three times at some µ = µ 1 , µ 2 and µ 3 leading to the following table of variations for

In such case, l(η) has three extrema: 2 maxima and a mimimum.

Finally, to specify in which situations l has three extrema, we remark that

So l(η) has three extrema if and only if α > 1+γ 2 and then if

To prove Proposition 3.2.2, assume that -1 < γ ≤ 1. f ′ 2 (µ) has the sign of:

which one denote by h(µ, γ), rewritten as:

We know from the lemma that f 2 has three extrema if and only if α > ∃ y 0 < 0 h(y 0 , γ) < 0 and y ′ 0 > 0 h(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 h(y 0 , γ) < 0 and h(y ′ 0 , γ) > 0. As lim 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 characterization using the function h of the existence of three extrema, assume there exist y 0 < 0 and y ′ 0 > 0 such that h(y 0 , 1) and h(y ′ 0 , 1) > 0. Then as γ ∈] -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, one deduces that ∀µ ≥ 0, h(µ, γ) ≥ h(µ, 1) and ∀µ ≤ 0, h(µ, γ) ≤ h(µ, 1), which means in particular that h(y ′ 0 , γ) ≥ h(y ′ 0 , 1), and that h(y 0 , γ) ≤ h(y 0 , 1), implying that h(y ′ 0 , γ) > 0 and h(y 0 , γ) < 0. Thus whatever γ ∈] -1, 1], the function h 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.

VII