Benjamin Baron

Marcelo Dias De Amorim

Vehicles as big data carriers: Road map space reduction and efficient data assignment

We advocate the use of a data shuttle service model to offload bulk transfers of delay-tolerant data from the Internet onto standard vehicles equipped with data storage capabilities. We first propose an embedding algorithm that computes an offloading overlay on top of the road infrastructure. The goal is to simplify the representation of the road infrastructure as raw maps are too complex to handle. In this overlay, each logical link maps multiple stretches of road from the underlying road infrastructure. We formulate then the data transfer assignment problem as a novel linear programming model that determines the most appropriate logical paths in the offloading overlay for a data transfer request. We evaluate our proposal using actual road traffic counts in France. Numerical results show that we can satisfy weekly aggregate requests in the petabyte range while achieving cumulative bandwidth above 10 Gbps with a market share of 20% and only one terabyte of storage per vehicle.

I. INTRODUCTION

We propose a disruptive offloading service model to handle the growth of demand for bandwidth-intensive background services that have strained the Internet infrastructure to its limits (such as maintenance activities, data migration, and online backups) [START_REF]Cisco Visual Networking Index[END_REF]. Our service exploits the delay-tolerant nature of the background Internet traffic to deliver massive amount of data over the conventional road infrastructure, relieving the Internet from this traffic while providing a reliable and robust transfer. The idea is to take advantage of the increasing number of journeys involving vehicles, which remain by far the most commonly used mode of transport. We argue in favor of equipping standard electric vehicles with data storage capability as part of the recent development in the automotive sector. Vehicles so equipped act as data shuttles delivering data in large quantities from the producer repositories to the destination points. We consider the model in which electric vehicles charge or replace their depleted batteries at charging or swapping stations, respectively. In our model, we propose that these stations serve as offloading spots where vehicles can load/deliver data (making all the process totally transparent to the driver). Thus, vehicles take turns carrying the data to its destination by loading and unloading it at offloading spots.

Related proposals suggest offloading Internet data onto other transmission media such as planes [START_REF] Keränen | DTN over aerial carriers[END_REF], postal services [START_REF] Cho | Budget-constrained bulk data transfer via internet and shipping networks[END_REF], or even birds [START_REF] Waitzman | A standard for the transmission of ip datagrams on avian carriers[END_REF]. These approaches allow saving some bandwidth from the Internet and even reducing the costs compared to other solutions such as renting a dedicated line. Nevertheless, they fail to achieve massive offloading and can be used only in very specific situations. As we will show later in this paper, our solution is more ambitious and can achieve rates of several gigabytes per second for a global amount of data reaching a petabyte per week.

A major challenge toward our vision is how to assign a data transfer request onto the road infrastructure. We provide a mathematical optimization formulation of this problem. We first propose a mapping algorithm that creates an offloading overlay network on top of the road infrastructure. Nodes in the overlay network are connected through virtual links, which correspond to logical paths connected by multiple stretches of road in the underlying road infrastructure. Each virtual link is characterized by a set of attributes such as delay, capacity, and data losses. The offloading overlay network provides an abstract view of the resources that may be allocated to a particular data transfer. Second, we formulate the data transfer assignment problem as a novel linear programming model that determines the optimal logical paths consisting of a set of offloading spots connected together by logical links in the offloading overlay.

We provide numerical results using actual road traffic counts in France. The offloading spots are determined to ensure full coverage of all data transfer requests. The search space is further reduced by restricting the physical paths that satisfy the delay requirements.

In summary, our main contributions in this paper are:

• Space reduction. We design a mapping algorithm that creates the offloading overlay that mitigates the complexity of the road infrastructure and gives a coherent representation of the underlying resources. • Data transfers. We formulate the data transfer assignment problem as a novel linear programing model. The model provides the optimal logical paths on the offloading overlay for a given data transfer request. • System capacity. We evaluate the offloading service model for the data transfer assignment problem using actual road traffic counts of the main roads of France. The results show impressive numbers, with potential transfers in the order of a petabyte during a one-week time window.

II. PRINCIPLES

We consider a bimodal data transportation service that combines the legacy Internet with flows of vehicles traveling roads. Our offloading service aims at overcoming the limitations of the Internet in terms of capacity when it comes to transfer large amounts of delay-tolerant data. We assume that vehicles are equipped with some storage capacity and a network of offloading spots provides the locations where vehicles can load/deliver data. Offloading spots could be battery charging and swapping stations as the ones Tesla is currently rolling out both in Europe and North America [START_REF] Motors | Supercharger[END_REF]. At offloading spots, data can be either transferred onto the storage devices of the vehicles or the latter can be exchanged for preloaded storage devices that match the destination of the vehicle. 1This offloading service allows to relieve the Internet of the background traffic by transferring it on the road infrastructure using vehicles that will take turns to carry the data to its destination. Hence, the background traffic will not hinder both the elastic and interactive traffic, which are under stricter constraints (in terms of delay and throughput) [START_REF] Hong | Achieving high utilization with software-driven wan[END_REF]. The service further provides a robust transfer of high volumes of data, as the transfer is spread over multiple independent vehicles going in different directions. Finally, we assume the on-board storage cannot be accessed by the drivers and the data is encrypted, preventing anyone from reading the content.

III. OFFLOADING SERVICE MODEL AND NOTATIONS

The problem we introduce is twofold: [START_REF]Cisco Visual Networking Index[END_REF] to propose a mapping algorithm to create the offloading overlay network on top of the road infrastructure and (2) to model the data transfer assignment problem with a linear programing model that aims at maximizing the revenue on the offloading overlay.

A. Offloading overlay

The road infrastructure is represented as a directed graph G R = (N R , L R ) where N R and L R refer to the set of nodes (junctions and stations where vehicles can charge/replace their depleted batteries) and links (stretches of road), respectively.

The road infrastructure has many characteristics to be taken into account, such as the flow of vehicles from an origin to a destination, the capacity, the free-flow speed, and the travel time of the stretches of road. We introduce an offloading overlay to aggregate these properties in a simplified representation of the road infrastructure. The offloading overlay is a bijection from the stations of the road infrastructure to offloading spots and a mapping from physical paths between pairs of stations to logical links connecting pairs of offloading spots.

We denote the offloading overlay by a directed graph

G O = (N O , L O ),
where N O and L O refer to the set of nodes (offloading spots mapping the stations) and links (logical links between the offloading spots), respectively. The logical links are characterized a logical by the weighted travel time t(i, j), aggregated capacity c(i, j), and data leakage L(i, j), that is the proportion of data lost on logical link (i, j) ∈ L O . In Fig. 1, we show an example of a possible realization of an overlay on top of the French road infrastructure.

B. Data transfer assignment problem

We consider data transfer requests, denoted R st , where data has to be delivered from a source offloading spot s ∈ N O to a target offloading spot t ∈ N O , within a maximal delay T st , and with a leakage tolerance L st . We denote by R the matrix of all requests indexed by all sources and destinations. The leakage tolerance 0 < L st 1 refers to the proportion of data that can be lost during the transfer. Finally, we denote by B st the total amount of data (including redundant data) transferred between s and t.

We assume that all offloading spots replicate data on multiple outgoing vehicles. For a given request R st , each offloading spot i replicates the data ρ st (i, j) times when allocated on logical link (i, j) ∈ L O (i.e., the same data is loaded onto ρ st (i, j) vehicles). Such replication is needed to mitigate the logical link leakage L(i, j) and satisfy the leakage tolerance L st of the request.

Given a request R st , our offloading service is subject to a gain factor γ st and an operational cost α i . γ st is the gain factor in terms of cost to send data from s to t using our offloading service over the legacy Internet (e.g., a dedicated connection between s and t). Without loss of generality, we assume that γ st is linear with the volume of data transferred (i.e., for one bit we transfer using our offloading service, we save γ st compared to the legacy Internet). α i is the operational cost at offloading spot i to handle the storage of data that is unloaded from incoming vehicles.

In the following, we state our data delivery problem as a linear optimization problem. We define the revenue as the difference between the gain we generate using the offloading service to transfer B st data within duration T st and the cost of handling the data at the offloading spots. This problem aims at maximizing the revenue, which is constrained by the delay tolerance T st and the capacity c(i, j) of the logical links.

IV. MAPPING OF THE OFFLOADING OVERLAY

A logical link (i, j) ∈ L O (see Section III-A) is characterized by a travel time t(i, j), an aggregated capacity c(i, j), and a data leakage L(i, j). In the following, we consider O-D flows between two stations, i.e., flows that go from station i ∈ N S to station j ∈ N S . Let r be the number of "reasonable" paths between i and j, i.e., paths that do not backtrack. 2 The value r depends on the traffic assignment strategy that we consider (e.g., All-or-Nothing, Wardrop's equilibrium, or Dial's stochastic assignments) [START_REF] De Dios Ortúzar | Modelling transport[END_REF], [START_REF] Sheffi | Urban transportation networks: equilibrium analysis with mathematical programming methods[END_REF].

A. Travel time of logical links

For stretch of road (a, b) ∈ L R , let v ab be its nominal volume of vehicles (vehicles per unit of time), c ab be its capacity (vehicles per unit of time), and t ab (0) be its corresponding travel time at free-flow speed (i.e., when v ab = 0). The travel time of (a, b) is given by the BPR (Bureau of Public Roads) function defined as [START_REF]Traffic assignment manual for application with a large, high speed computer[END_REF]:

t ab (v ab ) = t ab (0) 1 + α v ab c ab β , (1)
where α and β are BPR parameters that depend on the road profile (α = 0.15 minutes and β = 4.0 are typically used) [START_REF]Transportation Research Board, Highway Capacity Manual[END_REF].

We can deduce from Eq. ( 1) the travel time of physical path p, denoted t p , which is the sum of all travel times of the stretches of road the path consists of (we do not consider any turning delays at junctions):

t p = (a, b)∈p t ab (v ab ).
(

) 2
From Equation (2), we deduce an expression of the average travel time t(i, j) experienced on the r physical paths between nodes i and j, weighted by the road traffic flow v p on each path p:

t(i, j) = p∈P ij t p v p r p∈P ij v p , (3)
where P ij is the set of all simple physical paths between i and j (i.e., with no cycles in the path).

B. Capacity at logical links

The capacity c(i, j) of the overlay link (i, j) ∈ L O depends on the sum of the traffic flows v p of the simple paths between offloading spots i and j (i.e., the number of vehicles per unit of time going from i to j on path p). It also depends on the market share M of the vehicles participating in the offloading service and the storage size S on each vehicle:

c(i, j) = M × S p∈P ij v p . (4)

C. Leakage of logical links

The leakage L(i, j) (comprised between 0 and 1) of logical link (i, j) ∈ L O represents the proportion of data that is lost on logical link (i, j). The leakage increases as more vehicles carrying data prematurely exit the road (e.g., the vehicles may exit the highway before reaching the offloading spot or an accident may have occurred).

V. OFFLOADING SERVICE REVENUE MAXIMIZATION

We denote by p a logical simple path (i.e., without any cycles) defined as a set of offloading spots connected by logical links. We denote by P st the set of all logical simple paths between s and t.

We define the travel time t(p) experienced on logical path p ∈ P as the sum of the travel times t(i, j) on the logical links that compose path p and the waiting time t i at each intermediate offloading spot on the path:

t(p) = (i, j)∈p t(i, j) + t i . (5)
We denote by f (p) the resulting flow (bandwidth) on logical path p ∈ P st for a given request R st . For the same request, the transfer of an amount of data B st is constrained by the delay tolerance T st and does not depend on replication of data. Hence, the delay constraint is the same for both models. The delay to transfer B st between s and t is expressed as the sum of the duration to transfer this quantity over all logical paths between s and t and the average travel time experienced on these paths, weighted by the flow on each path. The delay constraint is expressed as:

p∈P st f (p)t(p) p∈P st f (p) + B st p∈P st f (p) T st . (6)
Eq. ( 6) is not a linear expression, which can lead to some difficulties to solve the linear optimization problem. But it can be rewritten as the following linear expression:

p∈P st f (p) T st -t(p) B st . (7)
In the model we consider, offloading spots store and replicate data. Consider a given request R st with leakage tolerance L st and p ∈ P st a path between s and t. If ρ st (i, j)f (p) data is transmitted on overlay link (i, j) ∈ p, f (p) data is received at destination offloading spot j. Here, the replication factor ρ st (i, j) is calculated by each offloading spot i that sends data to remote offloading spot j on logical link (i, j) ∈ L O as a function of L st and L(i, j) as follows:

L(i, j) ρ st (i, j) L st =⇒ ρ st (i, j) log(L st ) log(L(i, j)) • (8)
The flow on logical link (i, j) ∈ L O is expressed as the sum of all flows sent by offloading spot i to offloading spot j on logical link (i, j) for data transfer request R st . The flow is constrained by the capacity of the logical link c(i, j): s, t∈R p∈P st p∋(i, j)

ρ st (i, j) f (p) c(i, j)• ( 9
)
The revenue is expressed as the gain of our offloading service compared to an Internet technique (using the gain factor γ st ) minus the operational cost (costs α i at each offloading spot i). Our objective is to maximize the total revenue, expressed as:

s, t∈R B st γ st - i∈N O α i p∋i f (p). ( 10
)
We note that if Eq. ( 10) is optimal (i.e., maximized), then Eq. ( 7) becomes an equality since the quantities B st are constrained only in this equation and γ st > 0. Now, if we consider an equality instead of an inequality for the delay constraint Eq. ( 7), when replacing B st with its expression, the total revenue objective can be rewritten as:

s, t∈R p∈P st f (p) γ st T st -t(p) -α(p) , ( 11
)
where α(p) = i∈p α i . In Eq. ( 11), if γ st T stt(p)α(p) is negative, then f (p) is null (objective function decreases otherwise). Hence, from Eq. ( 11), we have:

t(p) + α(p) γ st T st =⇒ f (p) = 0. ( 12
)
Since Eq. ( 12) is a weight on the logical paths (resulting from the weights of the arcs) that depends only on the overlay graph, it allows narrowing the search space of the paths that do not satisfy Eq. ( 12). Finally, the formulation can be reduced to the total revenue maximization objective that is subjected to the capacity constraint on each logical link

(i, j) ∈ L O : Maximize s, t∈R p∈P st f (p)ψ st (p)

Subject to

s, t∈R p∈P st p∋(i, j)

ρ st (i, j)f (p) c(i, j) ∀(i, j) ∈ L O
where ψ st (p) = γ st T stt(p)α(p) can be considered as a weight on logical path p ∈ P st . The amount of data transferred B st within the delay tolerance imposed by request R st is deduced from Eq. ( 7):

B st = p∈P st f (p) T st -t(p) • (13)
We can finally obtain from Eq. ( 13) the average bandwidth of request R st by dividing the amount of data B st by the duration of the transfer T st .

VI. EVALUATION

We evaluate our offloading service model with real traffic datasets measured in France's road infrastructure.

A. Dataset and planning of the network of stations

We use road traffic data collected in 2011 on the major stretches of road in France [START_REF]Recensement de la circulation sur le réseau routier national en 2011[END_REF]. The dataset shows AADT (Annual Average Daily Traffic), i.e., the total volume of traffic measured on both directions of the stretches of road for a year divided by the number of days in the year. The graph is composed of 3,310 edges covering over 20,000 km of roads. We consider a network of stations equivalent to the one Tesla is currently rolling out in North America [START_REF] Motors | Supercharger[END_REF]. Since there is no such a network in France for the time being, we plan a simple yet realistic network of stations over the road infrastructure of France. We consider a facility-allocation problem that minimizes the number of facilities to allocate, a problem we adapted from the maximal covering location problem [START_REF] Church | The maximal covering location problem[END_REF]. The problem takes demand points and candidate locations as inputs: (i) the demand points are the 9,555 cities of France with a population greater than 1,000, and (ii) the candidate locations are the 1,024 gas stations of Total, a major oil company in France.

The facility-allocation algorithm chooses the stations such that maximum demand points are allocated to the stations within a range of 150 km, while minimizing the number of chosen stations. The resulting allocation gives 38 stations scattered on the French road infrastructure as shown in Fig. 2a.

B. Mapping of the offloading overlay

To compute the overlay depicted in Fig. 2b, we consider the All-or-Nothing traffic assignment strategy that assigns all the traffic between a source and a destination to the shortest path, here defined as the path with the shortest travel time [START_REF] De Dios Ortúzar | Modelling transport[END_REF], [START_REF] Sheffi | Urban transportation networks: equilibrium analysis with mathematical programming methods[END_REF]. The traffic flow on this path is set to the minimum AADT of the stretches of road composing it. Then, we use the results of the most recent French travel household survey made in 2008 (ENTD), which distinguishes two kinds of travel [START_REF] Medde | Mobilité -déplacement[END_REF]: local travel is less than 80 km from home, while long-distance travel is greater than 80 km from home. The latter one accounts for 2.8% of the total amount of recorded travel. We consider a market share M of 20% (which is the case for Renault in France) and a per-node storage unit S of one TB (assumed to be the same for each electric vehicle).

Since the AADT is measured on both directions of a stretch of road, we divide the measured values by two. To compute the bandwidth of the logical link, we multiply the resulting flow by S. We use the BPR (Eq. ( 1)) to calculate the travel time on each stretch of road with: α = 0.15 minutes, β = 4.0, the practical capacity given by the dataset, and the AADT for the assigned volume. We use then Eq. (3) to deduce the travel time of the logical links using the shortest physical road paths. Fig. 3 shows the benefits of using the offloading overlay compared to the road infrastructure. We note that the number of simple paths grows exponentially with the travel time cutoff and is much larger on the road infrastructure. Also, the difference in the number of paths grows exponentially.

C. Data transfer assignment problem

We interface our offloading overlay with IBM ILOG CPLEX 12.5.1 using the linear optimization models we presented in Section V. We consider a scenario with the following three different requests (distances are Euclidean):

• Request R1: from Paris to Lyon (384 km).

• Request R2: from Paris to Bordeaux (492 km).

• Request R3: from Paris to Marseille (646 km). It is important to note that requests R1 and R3 will compete for the flows since they share some common subpaths. We use a breadth-first search algorithm to generate all the simple paths for each request. The cutoffs of the simple paths are given by Eq. [START_REF] Church | The maximal covering location problem[END_REF]. However, the generation of the simple overlay paths between a source and a destination is exponential, as seen in Fig. 3. To solve this issue, we reduce our search space by applying a default cutoff of 12 hours on the travel time of the simple logical paths we generate for our experiments.

We express the gain factor γ st as an exponential function of the Euclidean distance (in km) between s and t: γ st = dist(s, t) β , where β ∈ R. In our analyses, we investigate the impact of the delay tolerance T st . By default, for all requests, we set leakage tolerance L st to 10 -2 and link leakage L(i, j) to 0.05 for each logical link (i, j) ∈ L O . The operational cost α i is weighted by the demands allocated to offloading spot i, resulting from the facility-allocation problem. Since we want a fair flow allocation in our offloading service, we choose β such that the standard deviation of the bandwidth of all flows is minimized. Finally, for each offloading spot i, we consider a waiting time t i = 20 minutes, which corresponds to the duration of a charge that provides a 300-km range to the vehicles [START_REF] Motors | Supercharger[END_REF]. Fig. 4 shows the average bandwidth of each request as a function of the delay tolerance T st . We notice that the bandwidth stabilizes as the duration of the transfer increases. Also, requests with longer distances are favored (request R2 and R3) over requests with shorter distances (request R1) since the total revenue is increased. VII. CONCLUSION We proposed a bulk-data transfer service to offload delaytolerant Internet traffic onto the road infrastructure network connecting geographically distant offloading spots. We evaluated our system using real traffic numbers in France. Our offloading service allows data transfers in the petabyte range per week while achieving cumulative bandwidth exceeding 10 Gbps with a market share of 20% and only one terabyte of storage per vehicle. These results confirm the offloading potential of our service which can help operators handle big data, especially during busy periods.

As future work, we plan to replicate the data partially at some offloading spots or at the source offloading spot. We also intend to address security and privacy concerns, as well as the adaptation of our models so as to capture diurnal and seasonal variations of the road traffic.

Fig. 1 :

 1 Fig. 1: France main roads infrastructure (the darker and bolder the stretches of road, the higher the AADT) and the offloading overlay above. The AADT (Annual Average Daily Traffic) is the volume of traffic over a year divided by the number of days in the year.

  Offloading overlay with the average capacity of overlay links.

Fig. 2 :

 2 Fig. 2: Facility-allocation result and offloading overlay. The big dots are the chosen stations.

Fig. 3 :

 3 Fig. 3: Total number of simple paths in the road infrastructure and logical simple paths in the offloading overlay as a function of the travel time cutoff (y-axis is in a logarithmic scale).

Fig. 4 :

 4 Fig. 4: Bandwidth as a function of delay tolerance T st

The authors are conscious that there are some privacy issues involved in this process, but they are out of the scope of this paper.

A path that does not backtrack is a path that takes a traveler away from origin and closer to destination.

ACKNOWLEDGMENT

Benjamin Baron, Promthe Spathis, and Marcelo Dias de Amorim carried out part of the work at LINCS (http://www.lincs.fr).


Thomas Oberlin
email: thomas.oberlin@inria.fr

Sylvain Meignen
email: sylvain.meignen@imag.fr

Valérie Perrier
email: valerie.perrier@imag.fr

Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations

Keywords: time-frequency, reassignment, synchrosqueezing, AM/FM, multicomponent signals. EDICS Category: DSP-TFSR, DSP-TRSF

This paper considers the analysis of multicomponent signals, defined as superpositions of real or complex modulated waves. It introduces two new post-transformations for the short-time Fourier transform, that achieve a compact time-frequency representation while allowing for the separation and the reconstruction of the modes. These two new transformations thus provide the benefits of both the synchrosqueezing transform (which allows for reconstruction) and the reassignment method (which achieves a compact time-frequency representation). Numerical experiments on real and synthetic signals demonstrate the efficiency of these new transformations, and illustrate their differences.

Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations I. INTRODUCTION Most real signals can be accurately modeled as superpositions of locally bandlimited, amplitude and frequency-modulated (AM-FM) waves. Containing several modes they are often called multicomponent signals (MCS). Many approaches have been built to decompose these signals into their constituent modes either in time-frequency (TF) domain, like the synchrosqueezing transform [START_REF] Daubechies | A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models[END_REF], or in time domain, as for instance the empirical mode decomposition (EMD) [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], [START_REF] Rilling | On empirical mode decomposition and its algorithms[END_REF]. Yet, very few results on the accuracy of the mode retrieval are available. In this regard, a new approximation result on the retrieval of weakly modulated modes of an MCS has been shown using the synchrosqueezing transform (SST) [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF]. To obtain such a result, very strong assumptions were made on the modes making up the signal among which the most limitating one is the weak modulation hypothesis on the modes. Indeed, most MCS contain strongly frequency modulated modes, as for instance chirps involved in radar [START_REF] Skolnik | Radar handbook[END_REF], speech processing [START_REF] Pitton | Applications of positive time-frequency distributions to speech processing[END_REF], or graviational waves [START_REF] Candes | Detecting highly oscillatory signals by chirplet path pursuit[END_REF]. Although many TF transforms such as quadratic distributions [START_REF] Cohen | Time-frequency analysis: theory and applications[END_REF], [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF] or the reassignment method (RM) [START_REF] Kodera | A new method for the numerical analysis of non-stationary signals[END_REF], [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF] manage to accurately represent this kind of signals, they do not allow for mode separation and reconstruction. Recently, some work has been done on how to deal with high FM with SST. For instance, in [START_REF] Li | A generalized synchrosqueezing transform for enhancing signal time-frequency representation[END_REF], authors use a two-steps algorithm, which first computes SST, estimates the FM, and finally recomputes SST on the demodulated mode. In the same spirit, authors in [START_REF] Wang | Matching demodulation transform and synchrosqueezing in time-frequency analysis[END_REF] use an iterative procedure, where at each iteration a SST is computed with a finer resolution. A different and earlier attempt [START_REF] Gardner | Sparse time-frequency representations[END_REF] consisted in building an invertible extension of RM. Unfortunately, the proposed reconstruction was only an approximation, valid for Gaussian windows with a small temporal width.

The aim of this paper is to extend SST to MCS containing strongly frequency modulated modes by generalizing the SST in two different ways. The first attempt is inspired by RM, which is known to optimally sharpen the TF representation of a linear chirp. However, in the original RM the sharpened representation obtained from the spectrogram did not allow for mode retrieval. We will see that by fully exploiting the information contained in the so-called reassignment vector, one can derive a new May 22, 2014 DRAFT generation of SST that provides a sharpened TF representation as good as the one given by RM, while enabling exact retrieval of linear chirps. Our second approach relies on a second-order local estimate of the instantaneous frequency, which is used to improve modes localization and reconstruction using SST.

The outline of the paper is as follows: in section II, we recall important notations that are subsequently used in the paper and we introduce the concept of weakly modulated modes. We then define SST and RM in section III, putting the emphasis on the differences in terms of reconstruction and representation.

Then, section IV is devoted to the presentation of the new extensions of SST. Section V finally delivers numerical results on both synthetic and real data, comparing the proposed methods with standard SST and RM.

II. TIME-FREQUENCY ANALYSIS OF MULTICOMPONENT SIGNALS

A. Short-Time Fourier Transform

Definition II.1. The Fourier transform of a function f ∈ L 1 (R) is defined by

f (η) = R f (t)e -2iπηt dt. (1)
For a given η, f (η) represents the part of f that oscillates at frequency η on the whole time domain.

The need for time-localized frequency descriptors leads to the following definition:

Definition II.2. Let g ∈ L 2 (R) a real-valued and even function, with unit norm. The short-time Fourier transform (STFT) of f ∈ L 2 (R), with respect to the window g, is defined by

V g f (η, t) = f, g η,t = R f (τ )g(τ -t)e -2iπη(τ -t) dτ. ( 2
)
The representation of the positive quantity |V g f (η, t)| 2 in the TF plane is called the spectrogram of f . Without ambiguity, V g f will be denoted by V f .

Remark:

The usual definition of STFT differs from equation ( 2) by a factor e 2iπηt and considers the conjugate of g. However, to be consistent with the continuous wavelet transform formulation that follows we stick to equation [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. Also, as we assume g to be real-valued, we do not need the conjugate g * in (2).

It is well known from the properties of the Fourier transform that a well-localized TF representation requires a window g well localized both in time and frequency. We shall then mention that STFT is invertible on L 2 (R) using either of the following formulae:

f (τ ) = R 2 V f (η, t)g(τ -t)e 2iπη(τ -t) dtdη (3)
f (τ ) = 1 g(0) R V f (η, t)dη. (4)
Note that equation ( 4) is valid only if g is non-zero and continuous at 0. If the signal is analytic (i.e.

η ≤ 0 ⇒ f (η) = 0) then the integral domain for η is R + .

B. Continuous Wavelet Transform

STFT can be interpreted as the scalar product of f onto the translated and dilated versions of window g, each basis functions having the same TF resolution. On the contrary, the continuous wavelet transform (CWT) is associated with the decomposition of a signal on a family of functions with varying TF resolution [START_REF] Grossmann | Decomposition of Hardy functions into square integrable wavelets of constant shape[END_REF], [START_REF] Daubechies | The wavelet transform, time-frequency localization and signal analysis[END_REF]. Let us recall briefly its definition:

Definition II.3. An admissible wavelet is a function ψ ∈ L 2 (R) satisfying the condition 0 < C ψ = ∞ 0 | ψ(η)| 2 dη η < ∞. (5)
Definition II.4. Let f ∈ L 2 (R) and ψ be an admissible wavelet, the continuous wavelet transform (CWT) of f is defined for each (a, t) ∈ R * + × R by:

W f (a, t) = f, ψ a,t = 1 a R f (τ )ψ τ -t a * dτ. ( 6
)
The representation |W f (a, t)| 2 in the time-scale plane is called the scalogram of f .

It is well known [START_REF] Grossmann | Decomposition of Hardy functions into square integrable wavelets of constant shape[END_REF] that when ψ is real, the transform is invertible:

f (τ ) = 1 C ψ R * + ×R W f (a, t) 1 a ψ τ -t a dadt a . (7)
As in the STFT case, another reconstruction formula exists, which consists in the following integral:

f (τ ) = 1 C ′ ψ ∞ 0 W f (a, τ ) da a , (8)
where the constant

C ′ ψ = R ψ(η)
η dη is supposed to be finite. Similarly to the STFT case, a good TF localization imposes some constraints on the wavelet ψ.

To properly deal with real signals, one needs to associate with these signals only positive frequencies, that is to consider the analytic signal f a = (I + iH)f , where Hf is the Hilbert transform of f . Then, May 22, 2014 DRAFT when one uses analytic wavelets (i.e. satisfying η ≤ 0 ⇒ ψ(η) = 0), f can be reconstructed from the CWT of its analytic signal f a since:

Proposition II.1. Let f ∈ L 2 (R) be a real-valued signal, and ψ an admissible analytic wavelet (satisfying ν ≤ 0 ⇒ ψ(ν) = 0). Denote by f a = (I + iH)f the analytic signal associated to f , then W f = 1 2 W fa . In particular, the following synthesis stands in L 2 (R):

f (τ ) = Re 2 C ψ ∞ 0 R W f (a, t) 1 a ψ τ -t a dtda a , (9)
This result shows that the use of an analytic wavelet enables to deal with real signals the same as complex ones.

C. MCS and Ridges

Definition II.5. An AM-FM mode is an oscillating function h(t) = A(t)e 2iπφ(t) , with A and φ ′ positive and slow-varying.

A(t) is called the instantaneous amplitude of h at time t, φ ′ (t) being its instantaneous frequency. Note that the slow-varying condition is not quantified, since as we will see later, it is model-dependent. Thanks to that condition, a first order expansion of the phase combined with a zeroth order expansion of the amplitude lead to the approximation, for τ close to a fixed time t:

h(τ ) ≈ ht (τ ) = A(t)e 2iπ[φ(t)+φ ′ (t)(τ -t)] . (10)
STFT and CWT of mode h are then approximated by:

V h (η, t) ≈ V ht (η, t) = h(t)ĝ(η -φ ′ (t)) W h (a, t) ≈ W ht (a, t) = h(t) ψ(aφ ′ (t)) * . (11)
This shows that the STFT (resp. CWT) of an AM-FM mode is non-zero on a TF (resp. time-scale) strip, centered on the ridge corresponding to its instantaneous frequency and defined by η = φ ′ (t) for STFT and a = ηψ φ ′ (t) for CWT, η ψ = arg max η | ψ(η)| being the center frequency of the wavelet. Without loss of generality, we shall assume in the sequel that η ψ = 1. The width of this strip depends on that of the support of ĝ (resp. ψ). Let us now define MCS, that are superpositions of AM-FM modes.

May 22, 2014 DRAFT Definition II.6. A MCS reads

f (t) = K k=1 f k (t) = K k=1 A k (t)e 2iπφk(t) ,
where the modes

f k satisfy φ ′ K > • • • > φ ′ k > • • • > φ ′ 1 > 0.
In order to retrieve each mode from the STFT (resp. CWT) of f , we assume the modes occupy nonoverlapping regions in TF (resp. TS) plane. This requires the following frequency separation condition:

Definition II.7. If supp ĝ ⊂ [-∆, ∆], a MCS f is separated for STFT if the instantaneous frequencies satisfy for each k ∈ {1, • • • , K -1}: φ ′ k+1 -φ ′ k > 2∆. If supp ψ ⊂ [1 -∆, 1 + ∆], a MCS f is separated for the CWT if for each k ∈ {1, • • • , K -1}: φ ′ k+1 -φ ′ k φ ′ k+1 + φ ′ k > ∆.
Note that this separation condition is of a linear (resp. logarithmic) type for STFT (resp. CWT). Note also that for non-compactly supported windows (e.g. the Gaussian window) or wavelets, these conditions can be adapted by defining ∆ as the "essential support".

III. SYNCHROSQUEEZING AND REASSIGNMENT IN A NUTSHELL

The so-called reassignment method (RM) [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF] is a general way to sharpen a TF representation towards its ideal TF representation which corresponds for MCS given by definition II.6 to:

TI f (η, t) = K k=1 A k (t) s δ(η -φ ′ k (t)), (12)
s being equal to 1 or 2 depending on whether a linear or a quadratic representation is used (e.g. s = 2

for the spectrogram and s = 1 for STFT).

May 22, 2014 DRAFT

A. Reassignment of the Spectrogram

RM presented here in the spectrogram context consists of the computation of so-called reassignment operators:

ωf (η, t) = 1 2π ∂ t arg V f (η, t) = η + Im 1 2π V g ′ f (η, t) V g f (η, t) τf (η, t) = t - 1 2π ∂ η arg V f (η, t) = t + Re V tg f (η, t) V g f (η, t) , (13)
where ωf can be viewed as an approximation of the instantaneous frequency and τf as the group delay.

Then, the reassignment step aims to retrieve the ideal TF representation by moving the coefficients of the spectrogram according to the map (η, t) → (ω f (η, t), τf (η, t)), which reads:

S f (ω, τ ) = R ∞ 0 |V f (η, t)| 2 δ(ω -ωf (η, t))δ(τ -τf (η, t)) dηdt. (14)

B. STFT-Based SST

The aim of SST based on STFT (FSST), introduced in [START_REF] Daubechies | A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models[END_REF] in the CWT context, is to retrieve the ideal TF representation, i.e. T I f with s = 1, from its STFT. Following [START_REF] Daubechies | A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models[END_REF], the authors in [START_REF] Wu | Adaptive analysis of complex data sets[END_REF], [START_REF] Oberlin | The Fourier-based synchrosqueezing transform[END_REF] proposed to neglect operator τf , and to reassign only the complex coefficients V f (η, t) according to the map (η, t) → (ω f (η, t), t), which read:

T f (ω, t) = 1 g(0) ∞ 0 V f (η, t)δ (ω -ωf (η, t)) dη. ( 15
)
Knowing φ ′ k , the k th mode can then be retrieved by considering:

f k (t) ≈ |ω-φ ′ k (t)|<d T f (ω, t)dω. (16)
In a nutshell, SST reassigns the information in the TF plane and then makes use of this sharpened representation to recover the modes. The rationale behind equations ( 15) and ( 16) is the reconstruction formula

f (t) = 1 g(0) ∞ 0 V f (η, t) dη, (17)
but only considering a mode-related integration domain instead of R + . This enables us to retrieve the modes which was not possible with RM based on the spectrogram, but the representation is unfortunately less sharp as soon as the frequency modulation is not negligible. This difference between FSST and RM is further illustrated in section III-D.

May 22, 2014 DRAFT

C. Perfect Localization and Reconstruction

We now recall the main properties of FSST and RM in terms of TF localization and mode reconstruction. We first highlight for what type of signals they actually lead to the ideal TF representation and then consider a more general case.

1) Synchrosqueezing for Perturbations of Pure Waves: One can show that ωf actually matches the instantaneous frequency only for a pure tone Ae 2iπωt . In this case, FSST provides the ideal TF representation given by equation [START_REF] Li | A generalized synchrosqueezing transform for enhancing signal time-frequency representation[END_REF]. More interesting is the study of the behavior of FSST on perturbations of pure waves as was done in [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF] for SST based on CWT. The corresponding result is available in [START_REF] Oberlin | The Fourier-based synchrosqueezing transform[END_REF] and can be summarized as follows. Let us start by defining precisely the class of MCS we study.

Definition III.1. Let ε > 0 and ∆ ∈ (0, 1). We define the set B ∆,ε of MCS where

• for all k, f k satisfies: A k ∈ C 1 (R) L ∞ (R), φ k ∈ C 2 (R), sup t φ ′ k (t) < ∞ and for all t, A k (t) > 0, φ ′ k (t) > 0, |A ′ k (t)| ≤ ε and |φ ′′ k (t)| ≤ ε. • the f k are separated with resolution ∆ i.e. for all k ∈ {1, • • • , K -1} and for all t, φ ′ k+1 (t) -φ ′ k (t) > 2∆.
We can then define the synchrosqueezing operator, using a function ρ ∈ D(R), the space of compactly supported smooth functions, such that ρ = 1, a threshold γ > 0 and an accuracy parameter δ > 0:

T δ,γ f (ω, t) = 1 g(0) |Vf (η,t)|>γ V f (η, t) 1 δ ρ ω -ωf (η, t) δ dη, (18)
and finally state the main result.

Theorem III.1. Let f ∈ B ∆,ε and ν ∈ (0, 1 2 ). Consider a window g ∈ S(R), the space of smooth functions with fast decaying derivatives, such that

supp ĝ ∈ [-∆, ∆]. Then, if ε is small enough, • |V f (η, t)| > ε ν only if there exists k ∈ {1 • • • , K} such that (η, t) ∈ Z k := {(η, t) / |η-φ ′ k (t)| < ∆}. • For all k ∈ {1, • • • , K} and all pair (η, t) ∈ Z k such that |V f (η, t)| > ε ν , one has |ω f (η, t) -φ ′ k (t)| ≤ ε ν . ( 19
)
• For all k ∈ {1, • • • , K} there exists a constant C such that for all t ∈ R,

lim δ→0 |ω-φ ′ k (t)|<ε ν T δ,γ f (ω, t) dω -f k (t) ≤ Cε ν . (20)
Equation [START_REF] Carmona | Multiridge detection and time-frequency reconstruction[END_REF] shows that the synchrosqueezing operator has to be summed up around the ridge, the May 22, 2014 DRAFT width of the interval being 2ε ν . As in practice ε is not known, one chooses a parameter d instead. Mode f k is hence estimated at time t from FSST T f according to:

f k (t) ≈ φ ′ k (t)+d φ ′ k (t)-d T f (ω, t) dω. (21)
Note that a similar approximation result was stated before in [START_REF] Wu | Adaptive analysis of complex data sets[END_REF], but with stronger (global) assumptions on the mode separation f k , instead of the local one in Def III.1.

2) Reassignment of Perturbations of Linear Chirps: It was shown in [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF] that the reassignment operators perfectly localize linear chirps, i.e. modes with a linear instantaneous frequency and constant amplitude. More precisely, one can easily prove the following result.

Theorem III.2. Consider a pure linear chirp h(t) = Ae 2iπφ(t) , where φ is a polynomial of degree 2 and A > 0. Then, whenever the reassignment operators are defined one has ωh (η, t) = φ ′ (τ h (η, t)).

As the operators are local, this suggests that whenever the instantaneous frequency is quasi-linear and the amplitude slow-varying, the reassignment provides an accurate approximation of the ideal TF representation.

D. Illustrations

To illustrate the differences between FSST and RM, we use two synthetic test signals throughout this paper. One (signal 1) is made of low-order polynomial chirps, that behave locally as linear chirps, and the other one (signal 2) contains strongly nonlinear sinusoidal frequency modulations and singularities (points where φ ′′ = 0). Figure 1 displays FSST and RM for these signals. Note that wherever the slope of the ridge φ ′′ (t) is strong, RM gives a more concentrated representation than FSST.

To illustrate the specificity of FSST regarding invertibility, we display on Figure 2 the reconstruction process of each test signals from FSST. To compute an estimation of the ridges (t, φ ′ k (t)), knowing the number K of modes, we use the same algorithm as in [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF] or [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], which computes a local minimum of the functional

E f (ϕ) = K k=1 - R |T f (t, ϕ k (t))| 2 dt + R λϕ ′ k (t) 2 + βϕ ′′ k (t) 2 dt, (22)
where (t, ϕ k (t)) is an estimator of the ridge (t, φ ′ k (t)), and λ and β are two positive parameters tuning the level of regularization. Other ridge detectors such as the one developed in [START_REF] Carmona | Multiridge detection and time-frequency reconstruction[END_REF] may of course be used instead but this is beyond the scope of the present paper. Looking at the results of Figure 2, it is worth noting that the quality of the reconstruction of a mode is highly dependent on its modulation. For instance mode 2 and 3 of signal 1 are properly retrieved but so is not mode 3 since it contains stronger modulations. Our goal is thus in what follows to improve mode reconstruction based on FSST when the modes are strongly modulated.

May 22, 2014 DRAFT

Test signal 1

IV. TOWARDS A SECOND-ORDER FSST

Recalling that FSST was designed to process superpositions of perturbed pure waves, by second order FSST we mean that we want to adapt SST to superpositions of perturbed linear chirps. This should enable us to obtain an invertible sharpened TF representation of the same quality of the one provided by RM. To build our extension of SST, we are going to first consider the case of linear chirps, i.e. modes with constant amplitude and quadratic phase.

May 22, 2014 DRAFT

A. Handling Strong Modulations: Invertible Reassignment or Second-Order Synchrosqueezing?

Let us investigate here the case of linear chirps h(t) = Ae 2iπφ(t) with φ(t) = a + bt + 1 2 ct 2 . For such signals, one has for each t and τ :

h(τ ) = h(t)e iπc(τ -t) 2 e 2iπφ ′ (t)(τ -t) , (23)
which leads to:

V h (η, t) = h(t) R g(τ )e iπcτ 2 e -2iπ(η-φ ′ (t))τ dτ. (24)
In such a case, the first reassignment operator writes

ωh (η, t) = φ ′ (t) - c 2π u(η -φ ′ (t)), ( 25
)
where u is the function defined by:

u(ξ) = d dξ arg R g(x)e iπcx 2 e -2iπξx dx . (26)
In the case of a Gaussian window g, one can show that u is linear [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], [START_REF] Oberlin | Analysis of strongly modulated multicomponent signals with the short-time Fourier transform[END_REF]. However, for a general smooth window g the exact nature of u is unknown, but one can remark that it is a smooth function.

Formula (25) is directly derived from the first equation of [START_REF] Wang | Matching demodulation transform and synchrosqueezing in time-frequency analysis[END_REF]. In the case of a linear chirp, one may write using (23, 24):

ωh (η, t) = 1 2π ∂ t argV h (η, t) = 1 2π ∂ t arg h(t) R g(τ )e iπcτ 2 e -2iπ(η-φ ′ (t))τ dτ = 1 2π ∂ t 2πφ(t) + arg R g(τ )e iπcτ 2 e -2iπ(η-φ ′ (t))τ dτ = φ ′ (t) + 1 2π ∂ t (η -φ ′ (t)) u(η -φ ′ (t)) = φ ′ (t) - c 2π u(η -φ ′ (t)).
Equation (25) clearly shows that ω h (η, t) = φ ′ (t) does not hold as soon as c = 0. We propose two different ways to deal with the non-negligible modulation c = φ ′′ (t). The first one consists in modifying the first reassignment operator by taking into account second-order terms, and then use this operator to define a new SST. Since reassignment will only take place in the direction of the frequencies, we will call it VSST (V for vertical) in the sequel.

May 22, 2014 DRAFT

The second one uses the local time delay τh defined in RM, which writes for linear chirps:

τh (η, t) = t - 1 2π u(η -φ ′ (t)), (27)
which is derived from the second equation of ( 13):

τf (η, t) = t - 1 2π ∂ η arg V f (η, t) = t - 1 2π ∂ η arg R g(τ )e iπcτ 2 e -2iπ(η-φ ′ (t))τ dτ = t - 1 2π u(η -φ ′ (t))
Remembering that φ ′ (t) = b + ct, we get from ( 25) and ( 27) that:

ωh (η, t) = φ ′ (τ h (η, t)), (28)
which means ωh (η, t) is the exact value of φ ′ but at a shifted time location given by τ (η, t). This is not surprising since RM is well known to perfectly localize linear chirps. The challenge here lies in designing a complex reassignment formula which allows for mode reconstruction. Since the proposed technique will use the two components of the reassignment vector we will call it OSST (O for oblique).

B. Local Instantaneous Modulation

Both VSST and OSST will use an estimation of the frequency modulation, qf , whose definition follows:

Definition IV.1. Let f be in L 2 (R). qf (η, t) is defined wherever V f (η, t) = 0 and ∂ t τf (η, t) = 0 by:

qf (η, t) = ∂ t ωf (η, t) ∂ t τf (η, t) . (29)
The following result shows that the local operator qf (η, t) perfectly estimates the frequency modulation c = φ ′′ (t) for a linear chirp.

Proposition IV.1. If h(t) = Ae 2iπφ(t) is a linear chirp, ie φ ′′ (t) = φ ′′ = 0, then for each (η, t) where qh is defined, one has

qh (η, t) = φ ′′ . ( 30
)
Otherwise qh (η, t) can be computed through the general formula: Regarding the computation of qh , recall that for any signal f ,

qh (η, t) = Im 1 2π (V g ′ h ) 2 -V g ′′ h V g h V xg h V g ′ h -V xg ′ h V g h , (31)
∂ t V g f (η, t) = -V g ′ f (η, t) + 2iπηV g f (η, t).
To lighten the notation, we will note V g h for V g h (η, t). Equation ( 13) gives

∂ t ωh = ∂ t η + Im 1 2π V g ′ h V g h = Im 1 2π ∂ t V g ′ h V g h -V g ′ h ∂ t V g h (V g h ) 2 = Im 1 2π (-V g ′′ h + 2iπηV g ′ h )V g h -V g ′ h (-V g ′ h + 2iπηV g h ) (V g h ) 2 = Im 1 2π (V g ′ h ) 2 -V g ′′ h V g h (V g h ) 2
.

Similarly one has

∂ t τh = ∂ t t + Re V xg h V g h = 1 + Re ∂ t V xg h V g h -V xg h ∂ t V g h (V g h ) 2 = 1 + Re (-V g h -V xg ′ h + 2iπηV xg h )V g h -V xg h (-V g ′ h + 2iπηV g h ) (V g h ) 2 = Re V xg h V g ′ h -V xg ′ h V g h (V g h ) 2
.

May 22, 2014 DRAFT

C. Vertical Second-Order Synchrosqueezing

Let h(t) = Ae 2iπφ(t) be a linear chirp with constant FM φ ′′ . Equations ( 25) and (27) show that for any time t,

φ ′ (t) = φ ′ (τ h (η, t)) + φ ′′ (τ h (η, t))(t -τh (η, t)) = ωh (η, t) + qh (η, t)(t -τh (η, t)), (32)
which means that the first order Taylor expansion of φ ′ can be written in terms of the first reassignment and frequency modulation operators. Using locally the chirp approximation for a given signal f , one can now improve the estimation of the instantaneous frequency as follows:

Definition IV.2. Let f ∈ L 2 (R)
, we define the second-order instantaneous frequency estimator of f by:

ωf (η, t) =      ωf (η, t) + qf (η, t)(t -τf (η, t)) if ∂ t τf (η, t) = 0 ωf (η, t) otherwise.
The vertical second-order synchrosqueezing (VSST) then consists in replacing ω by ω in standard SST:

T V f (η, t) = 1 g(0) R V f (η, t)δ (ω -ωf (η, t)) dη (33)
Since the coefficients are only moved vertically, the reconstruction is still achievable by summing up the SST coefficients in the vicinity of the ridge. Moreover, as all the operators are local, VSST should remain efficient for small perturbations of linear chirps. This will be further studied numerically in the next section.

D. Oblique Synchrosqueezing

An alternative way to modify FSST to take into account FM consists in moving the coefficients according to the reassignment vector field while keeping the phase information so as to allow signal reconstruction. A first attempt in this direction defined the following complex reassigned transformation [START_REF] Gardner | Sparse time-frequency representations[END_REF]:

T C f (ω, τ ) = R 2 V f (η, t)e iπ(ω+η)(τ -t) δ (ω -ωf (η, t)) δ (τ -τf (η, t)) dηdt. ( 34
)
Unfortunately, this definition allowed for an approximative reconstruction only for Gaussian windows, and with non-negligible errors as demonstrated by numerical results.

Let us now propose an alternative reassignment technique computing the phase shift in a different way

May 22, 2014 DRAFT from [START_REF] Gardner | Sparse time-frequency representations[END_REF] and that we will call oblique synchrosqueezing (OSST). The starting point is to remark that the STFT of a linear chirp h is invariant by translations along the direction of the ridge (ct, t) , t ∈ R modulo a phase shift. Indeed, one can write, using (24) and φ

′ (t) = b + ct, V h (η, t) = h(t) R g(s)e iπcs 2 e -2iπs(η-φ ′ (t)) ds = e 2iπφ(t) R g(s)e iπcs 2 e -2iπs(η+c(τ -t)-φ ′ (τ )) ds = V h (η + (τ -t)c, τ )e 2iπ(φ(t)-φ(τ )) = V h (η + (τ -t)c, τ )e 2iπ(t-τ )(φ ′ (τ )+ c 2 (t-τ )) . ( 35
)
Instead of reassigning the real value |V h (η, t)| 2 as in RM, we propose to move the complex quantity

V h (η, t)e 2iπ(τ -t)(φ ′ (τ )-c 2 (τ -t)) = V h (η + (τ -t)c, τ )
Contrary to VSST where one builds a first order Taylor expansion of the derivative of the phase, OSST directly involve the phase shift. For a linear chirp and when τ = τh (η, t), one has φ ′ (τ ) = ωh (η, t) so that one can then write:

V h (η, t)e 2iπ(τ -t)(ωh(η,t)-qh (η,t) 2 (τ -t)) = V h (η + (τ -t)c, τ ).
This phase-modified coefficients belong to the vertical line I ′ τ = {(η, τ ) / η ∈ R}. Finally, summing up all the phase shifted coefficients along the frequency axis will provide an approximation of h at time τ .

The principle of OSST is illustrated on Figure 3: the coefficients making up the dashed curve are phase-shifted while moved to point (ω, τ ). To obtain an approximation of the signal h(τ ), the phase of the coefficients are modified in order to be equal to their projection on the dashed dotted line. Let us now give a more formal definition of OSST.

Definition IV.3. Let f ∈ L 2 (R), its oblique synchrosqueezing (OSST) is the function

T O f (ω, τ ) = 1 g(0) R 2 V f (η, t)e iπ(2ω-qf (η,t)(τ -t))(τ -t) δ (ω -ωf (η, t)) δ (τ -τf (η, t)) dηdt. (36)
An approximation of mode f k is then given by its OSST on the ridge, i.e.

f k (t) ≈ T O f (φ ′ k (t), t). (37)
As for FSST and VSST, we would like OSST to remain efficient for perturbations of linear chirps, when the amplitude (resp. phase) is not exactly constant (resp. quadratic). Recall that for SST, the strategy to handle such perturbations was to integrate the coefficients around the ridge, so as to compensate for the error made by estimating the instantaneous frequency by ωf . However, in the case of OSST the coefficients are moved both in time and frequency, making this regularization impossible. We will see in the next section how this impacts the performance of the proposed technique.

V. NUMERICAL RESULTS

To achieve a quantitative comparison between the different methods, we will reuse the test signals defined in section III-D. Recall that test signal 1 is made of quasi-linear chirps that should be well dealt with by VSST or OSST, while test signal 2 contains stronger nonlinear frequency modulations as well as singularities where φ ′′ (t) vanishes, which makes the study of such signals more difficult. For the sake of simplicity, we only consider the case of FSST, but similar results could be obtained within the wavelet framework (see for instance [START_REF] Oberlin | The Fourier-based synchrosqueezing transform[END_REF] for a comparison). Note that in our simulations, we use a Gaussian window with parameter σ = 0.05, defined by g(t) = σ -1 2 e -πt 2 σ 2 .

On Figure 4 SST, VSST and OSST of both signals 1 and 2 are depicted. Comparing with Figure 1, one notices that both VSST and OSST provide a nice TF representation, as sharp as RM. In the next subsections, we will compare more quantitatively the three methods in terms of the sparseness of the representation (i.e. the number of nonzero coefficients) and accuracy of mode reconstruction. We also numerically investigate the robustness to noise, and finally compare the methods on some real data.

A. Quality of Representation

The main purpose of RM was to make TF representations more readable, more easily interpretable.

Since readability is somewhat subjective, we propose here to measure the readability as the amount May 22, 2014 DRAFT Indeed, the representation of test-signal 1 with each of these three methods is almost perfect, since one needs only 3 coefficients per time instant to recover the signal energy, which is consistent with the three modes making up the test signal. Comparatively, the energy of FSST grows rapidly to reach a value around 80% and then stagnates which means the signal energy cannot be retrieved with a reasonnable number of coefficients. The results are similar but less sharp for the second test-signal, since it contains stronger nonlinear frequency modulations.

In order to visualize the influence of noise, we carry out the same experiments when the test signals are contaminated by white Gaussian noise. We use a strong noise level (0 dB) and display the results on Figure 6. We observe that the increase of the energy is slower, since the coefficients corresponding to noise are spread out in the whole TF plane. However, OSST, RM and VSST still behave better than FSST, even though VSST is slightly less competitive.

May 22, 2014 DRAFT

B. Reconstruction of the Modes

The main advantage of FSST over RM is its invertibility. Let us recall that there are two possible reconstructions of the mode from FSST. The easiest one assumes T f is almost ideal, so that the whole information is concentrated on the ridge. The mode f k at time t is then estimated following

f k (t) ≈ T f (φ ′ k (t), t). (38)
FSST allows to achieve a more accurate reconstruction by using an integration around the ridge as a regularization step, according to the formula

f k (t) ≈ d -d T f (φ ′ k (t) + ξ, t) dξ. ( 39
)
Provided d is of the same order of magnitude as the perturbation defined in Theorem III.1, the latter ensures an asymptotically perfect reconstruction. We investigate here the accuracy of the reconstruction for each method starting with the reconstruction technique defined in (38 of direct reconstruction from FSST is very poor. Looking at the results for each mode, one observes that all the methods achieve a good reconstruction for mode 3 since it is not modulated, whereas FSST is unadapted for modes 1 and 2 which contain stronger modulations. The results are less convincing for test-signal 2: even if OSST and VSST improve the reconstruction compared to FSST, the accuracy is far from satisfactory. Indeed, as this signal contains strongly nonlinear frequency modulations, both OSST and VSST fail to represent the modes as a perfect sharp curve in the TF plane, some information being spread out around the ridges. This result is in line with the simulation of Figure 6 right, where on average 4 coefficients per time instant are needed to recover the most part of the energy.

In order to compensate for that default, one has to integrate the transformations when possible around the ridge, following equation (39). The results for such a reconstruction are displayed on Figure 7. For each test signal, the SNR of the reconstruction is displayed with respect to parameter d. This time we display the SNR associated with the reconstruction of the whole signal but considering the modes individually would lead to the same results. Figure 7 shows that VSST allows for a quasi-perfect reconstruction from very few coefficients (SNR = 50 dB for d ≥ 3). Conversely, with FSST, even if one uses a large d, i.e. many coefficients, one cannot retrieve the modes with a high accuracy.

We already saw that in the noise-free case, the reconstruction is all the better that parameter d is large. This is obviously not true when the signal is contaminated by noise: d must be large enough to take into account the signal coefficients, while remaining small enough so as not to include too much noise in the reconstruction. We illustrate this quantitatively by making the same experiments as for Figure 7 but in the noisy case. Note that in this case, the estimation of the ridges using formula (22) is worsened due to the presence of noise. The results of the ridge estimation are depicted on Figure 8, together with the reconstruction accuracy, for both FSST and VSST and with a moderate noise level (SNR = 12 dB). While the ridge estimation can be successfully carried out on FSST and VSST for test-signal 1, only VSST gives a correct estimation for test-signal 2 since the latter contains stronger FMs. As a result, the reconstruction of test-signal 2 from SST cannot be satisfactory. This study tells us that to take into account second order terms in the definition of VSST not only improves the quality of mode reconstruction but also of ridge detection.

C. Illustration on Real Signals

We now illustrate the qualitative improvement of the two new transformations on real data. We start by considering a bat echolocation call, available at http://dsp.rice.edu/sites/dsp.rice.edu/files/software/ batsignal.zip. spectrogram and FSST (parameter σ being set to 0.04). By zooming around a portion of one ridge, the improvement brought about these new methods in terms of concentration and compactness of the representation is undeniable. Note that the result provided by FSST remains satisfactory since the frequency modulation is quite slow compared to the size of the window.

To better illustrate that improvement, we then considered a portion of a speech signal, with parameter σ set to 0.025. Figure 10 shows the corresponding TF representations associated with the four different methods, as well as a zoom on a modulated portion. For both VSST and OSST, the energy of the representation is almost perfectly concentrated along the ridges whereas FSST does not provide an as sharp representation. Note that the representation given by OSST seems slightly better but we know by now that VSST provides with a more reliable reconstruction technique.

VI. CONCLUSION

This paper introduced the vertical and oblique second-order synchrosqueezing transformations to circumvent the limitations of standard SST and reassignment. We emphasized that OSST can be viewed as a complex version of the reassignment method, whereas VSST keeps the fixed-time structure of the standard synchrosqueezing transform, allowing for a better reconstruction through a regularization step.

We showed that both methods provide the ideal invertible representation in the case of a linear chirp, while remaining efficient on general multicomponent signals containing strong frequency modulations. All the experiments suggest that OSST is slightly more competitive for the representation purpose, whereas VSST allows for a better reconstruction.

Note that both the OSST and VSST can be extended readily to the wavelet setting. Further work should be devoted to the theoretical analysis of both methods for perturbations of linear chirps, as well as a precise understanding of the influence of noise on the reassignment operators and on the ridge estimation.

Figure 2 .

 2 Figure 2. Illustration of mode retrieval based on SST on test signals 1 (left) and 2 (right). From top to bottom: the three modes making up the test signals, their SST with the estimated ridges superimposed, and the reconstruction error. SNR after reconstruction for modes 1, 2 and 3 equals 14, 51 and 70 dB for test signal 1, and 5, 7 and 70 dB for test signal 2, respectively. We use d = 5 in the reconstruction (see equation (21)).

Figure 3 .

 3 Figure 3. Illustration of oblique synchrosqueezing (OSST).

Figure 4 .

 4 Figure 4. Test signals number 1 (top) and 2 (bottom). From left to right, magnitude of FSST, OSST and VSST.

Figure 5 .

 5 Figure 5. Normalized energy as a function of the number of coefficients for test-signal 1 (left) and 2 (right). Abscissa gives the number of coefficients over the size N of the signal, i.e. the mean number of coefficients for each column of the TF plane.

Figure 6 .

 6 Figure 6. Normalized energy as a function of the number of coefficients for the noisy (0 dB) test-signal 1 (left) and 2 (right).

Figure 7 .

 7 Figure 7. Reconstruction of the modes as a function of d, which is proportional to the size of the interval of integration. The results are in SNR, for test signal 1 (left) and 2 (right).

Figure 8 .

 8 Figure 8. From left to right: the ridges estimated from FSST, those estimated from VSST, and the quality of reconstruction measured as a function of d. We use again test signal 1 (top) and 2 (bottom).

Figure 9 .

 9 Figure 9. TF representations of a bat echolocation call. The representations are (from left to right) STFT, FSST, OSST and VSST. A zoom on strong modulations is displayed in the second row.

Figure 10 .

 10 Figure 10. TF representations of a speech signal. The representations are (from left to right) STFT, FSST, OSST and VSST. A zoom on a strongly modulated part is displayed on the second row.

  Let us first recall that for a linear chirp equation (28) holds, i.e ωh (η, t) = φ ′ (τ h (η, t)). Then differentiating this equality with respect to t, one obtains: ∂ t ωh (η, t) = φ ′′ (τ h (η, t))∂ t τh (η, t) or equivalently qh (η, t) = φ ′′ (τ h (η, t)). Now since φ ′′ is constant, this leads to the expected result.

	where V g ′′ h and V xg ′
	May 22, 2014	DRAFT

h are the STFT of h with windows g ′′ and x → xg ′ (x).

Proof.

  ). The results measuring the accuracy of mode reconstruction displayed in terms of SNR are displayed in TableI. We check that both OSST and VSST enable a good reconstruction of polynomial chirps (signal-test 1), whereas the quality Direct reconstruction for each method, and both test signals. The accuracy of mode retrieval is expressed in SNR.

	SST OSST VSST Gardner
	Test signal 1 3.5	17	18	11
	Mode 1 1.5	15	16	8.8
	Mode 2 2.2	17	19	12
	Mode 3 19	19	19	19
	Test signal 2 2.2	8.9	7.6	5.6
	Mode 1 0.92	7.4	5.9	3.8
	Mode 2 1.4	8.4	7.2	5.3
	Mode 3 57	23	57	23
		Table I
	May 22, 2014			DRAFT

May 22, 2014 DRAFT

The authors acknowledge the support of the French Agence Nationale de


K Abed-Meraim
email: karim.abed-meraim@univ-orleans.fr

M A Al-Radhawi

Ph Ravier

A New Windowing Method for the Parameter Estimation of Damped Sinusoids

Keywords: Exponentially Damped Sinusoids, Colored Noise, Cyclo-stationarity, Complex Exponential Windowing, Cramer-Rao Bound

This paper presents a preprocessing technique based on exponential windowing (EW) for parameter estimation of superimposed exponentially damped sinusoids. It is shown that the EW technique significantly improves the robustness to noise over two other commonly used preprocessing techniques: subspace decomposition and higher order statistics. An ad-hoc but efficient approach for the EW parameter selection is provided and shown to provide close to CRB performance.

I. INTRODUCTION

Parameter estimation of superimposed exponentially damped sinusoids from a set of noisy observation data is a problem one frequently encounters in many applications such as power systems [START_REF] Lovisolo | Modeling of electric disturbance signals using damped sinusoids via atomic decompositions and its applications[END_REF], audio modeling [START_REF] Boyer | Audio modeling based on delayed sinusoids[END_REF], mechanical systems [START_REF] Sarma | An adaptive nonlinear filter for online parameter estimation of undamped and damped sinusoids[END_REF], physics [START_REF] Berti | Mining information from binary black hole mergers: a comparison of estimation methods for complex exponentials in noise[END_REF], chemistry [START_REF] Papy | Exponential data fitting using multilinear algebra: The decimative case[END_REF], nuclear magnetic resonance imaging [START_REF] Rubtsov | Time-domain Bayesian detection and estimation of noisy damped sinusoidal signals applied to NMR spectroscopy[END_REF] and others. This problem has received considerable attention in the signal processing community in the last decades and is now considered by many to be "solved". Indeed, there are plenty of "solutions" available to this problem. The best known methods include the Prony's method [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise[END_REF], [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF], the iterative quadratic maximum likelihood (IQML) method [START_REF] Bresler | Exact maximum likelihood parameter estimation of superimposed exponential signals in noise[END_REF], [START_REF] Hua | The most efficient implementation of the IQML algorithm[END_REF], the matrix pencil (MP) method [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF], and state space-based method [START_REF] Lu | Internal model principle-based control of exponentially damped sinusoids[END_REF]. In addition to the above methods, there are many tools available to preprocess (i.e., clean up) the raw data before any of the above methods is applied. The most popular ones of these tools are subspace decomposition (SD) [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise[END_REF], [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF] and higher order statistics (HOS) [START_REF] Papadopoulos | Parameter estimation of exponentially damped sinusoids using higherorder statistics[END_REF], [START_REF] Ruiz | Parameter estimation of exponentially damped sinusoids using a higher order correlation-based approach[END_REF].

In the last few years, renewed interest in the sinusoidal estimation problem has been observed either for dedicated applications [START_REF] Boyer | Audio modeling based on delayed sinusoids[END_REF], [START_REF] Rubtsov | Time-domain Bayesian detection and estimation of noisy damped sinusoidal signals applied to NMR spectroscopy[END_REF] or for adverse situations such as irregular sampling [START_REF] Gudmundson | Spectral Estimation of Damped Sinusoids in the Case of Irregularly Sampled Data[END_REF], impulsive noise [START_REF] Altınkaya | Subspace-based frequency estimation of sinusoidal signals in alphastable noise[END_REF] or colored noise [START_REF] Stoica | Subspace-based frequency estimation in the presence of moving-average noise using decimation[END_REF].

In this work, we consider a colored noise context and we discuss the idea of preprocessing the data with complex exponential windowing. Although there are various ways to modify and improve the Prony's method, the IQML method, or the MP method with a preprocessing step, using SD or HOS is known to be the most effective in reducing the noise effect. However, the SD method, which separates "signal subspace" from "noise subspace", requires the assumption that the noise is either white, or colored with known covariance up to a scalar. The HOS method works well only with Gaussian noise and long available data length.

We show that the concept of cyclostationarity (CS) [START_REF] Gardner | Exploitation of spectral redundancy in cyclostationary signals[END_REF], [START_REF] Abed-Meraim | Parameter estimation of exponentially damped sinusoids using second order statistics[END_REF] or more generally exponential windowing (EW) can be easily applied to better preprocess the superimposed exponentially damped sinusoids. The CS or EW methods preserve the signal structure but mitigate the noise effectively as long as it is stationary.

Compared to the SD and HOS methods, the CS or EW methods require the least assumption on the noise and hence perform the best. This advantage will be clearly illustrated by simulation results. We propose simple approaches to optimize the EW parameters in such a way the signal power or the post-processing signal to noise ratio (SNR) is maximized. As a benchmark, we derive the Cramér-Rao bound (CRB) expression and use it to illustrate the closeness of the proposed method's performance to the CRB.

The organization of this paper is as follows. In section II, the idea behind using cyclo-stationarity for pre-processing and a generalized exponential windowing method are introduced. Section III is devoted to performance investigations. First, we review the CRB expression for the colored noise case. Then, we evaluate the asymptotic covariance matrix of the pre-processed noise term. The latter is evaluated for white noise and used to express the averaged signal-to-noise ratio (after data pre-processing) that is used for efficient EW parameter selection. In section IV, simulation results are presented demonstrating the superiority of the proposed method. Finally, we state the conclusion in section V.

II. PROPOSED SIGNAL PREPROCESSING METHODOLOGY

A. Data Model

The data y(n) under consideration is modeled as follows:

y(n) = x(n) + w(n) = L ∑ m=1 h m e bmn + w(n), n = 0, 1, • • • , N -1 (1)
where h m = a m e jθm , b m = α m + jω m with α m < 0, and w(n) denotes the stationary noise 1 . a m and θ m are respectively the amplitude and the initial phase of the m th sinusoid; its damping and frequency factors are respectively α m and ω m .

B. Preprocessing Using Cyclo-Stationarity

Considering the fact that exponentially damped signals have relatively short (effective) length, we apply the CS concept only in the context of second order statistics (SOS) as opposed to HOS. Define the k th lag (k ≥ 0) cyclo-correlation r β (k) at the cyclo-frequency β as:

r β (k) := N -k-1 ∑ n=0 y(n + k)y * (n)e jβn ( 2
)
where * denotes the complex conjugate transpose operator. Assuming large sample size2 , and using (1) in (2) yields

r β (k) = L ∑ m,l=1 h m h * l e kbm N -1-k ∑ n=0 e (bm+b * l +jβ)n + v(k) ≈ L ∑ m=1 A β (m)e kbm + v(k) (3)
where

v(k) = N -1-k ∑ n=0 w(n + k)w * (n)e jβn , A β (m) = L ∑ l=1 h m h * l 1 -e bm+b * l +jβ (4)
It is seen from ( 4) that the cyclo-correlation as function of k also consists of superimposed exponentially damped sinusoids and a noise term. The main advantage, in dealing with cyclo-correlation instead of correlation function as considered in standard sinusoidal estimation techniques [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF], [START_REF] Bresler | Exact maximum likelihood parameter estimation of superimposed exponential signals in noise[END_REF], is that noise contribution is considerably reduced in the former case. In fact, due to the stationarity assumption of the noise process, we have:

lim N →∞ 1 N N -1-k ∑ n=0 w(n + k)w * (n)e jβn = 0 when lim N →∞ 1 N N -1-k ∑ n=0 w(n + k)w * (n) = ρ(k)
ρ(k) being the k th correlation factor of the noise process. Generally, for additive colored noise, the signal to noise ratio (SNR) gain, can be considerable since, for ρ(k) ̸ = 0, we have

lim N →∞ ∑ N -1-k n=0 w(n + k)w * (n) ∑ N -1-k n=0 w(n + k)w * (n)e jβn = ∞.
Remark 1: In the case where the data length is relatively small, one can enhance the signal component by applying the integral of the cyclo-correlation over an interval of the cyclo-frequency:

∫ β1 β0 r β (k)dβ = L ∑ m=1 B β (m)e bmk (5)
where

B β (m) = L ∑ l=1 h m h * l [j(log(1 -e bm+b * l +jβ1 ) -log(1 -e bm+b * l +jβ0 )) + β 1 -β 0 ]. (6)
Note that this integration is equivalent to replacing the unit-norm exponential weight w(n) = e jβn by the weight coefficient

w ′ (n) = ∫ β1 β0 e jβn dβ = e j β 1 +β 0 2 n sinc ( β 1 -β 0 2 n )
representing the modulation of a unit-norm exponential with a sinc function (where sinc(x) = sin(x)/x).

Such integration has also the advantage to provide certain robustness against a 'bad' choice of the cyclo- Remark 2: In the case where the noise is complex circular which implies E(w(n)w(n + k)) = 0, we also can use the following cyclo-correlation: [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise[END_REF] in which the noise term is mitigated due to both the stationarity and the circularity.

r ′ β (k) def = N -1-k ∑ n=0 y(n + k)y(n)e jβn

C. Generalized Preprocessing Using Exponential Windows

A generalization of the previous preprocessing method is the use of a complex cyclo-frequency or a complex exponential window (EW). Therefore, we define the k th lag exponentially windowed signal r W (k) as:

r W (k) := N -1-k ∑ n=0 y(n + k)y * (n)W n ( 8
)
where W is a complex windowing parameter satisfying |W | ≤ 1. We retrieve the previous case when

|W | = 1.
Again by assuming large sample size, the expression of the exponentially windowed signal is:

r W (k) = L ∑ m=1 A W (m)e kbm + v(k) (9)
where

v(k) = N -1-k ∑ n=0 w(n + k)w * (n)W n , A W (m) = L ∑ l=1 h m h * l 1 -W e bm+b * l (10)
The rational behind using a complex windowing parameter is that the exponentially damped signal has a short effective length determined by the damping factor, and the rest of the data samples are almost noise. Therefore, using an exponential window will reduce accumulated noise in the resulting signal and provide higher weight to data samples with higher signal amplitude. This can lead, as will be observed later in section V, to non-negligible performance gain. Another reason is that if the signal frequencies were pure sinusoids, then the formula (3) will no longer be valid without a damped window W . The damped window |W | < 1 will guarantee the geometric series to converge always.

III. PERFORMANCE INVESTIGATION

In this section, we review the CRB for the estimation based on the exponentially windowed signal (EW-CRB). This will serve two objectives, the first is to compare the EW-CRB with the genuine CRB computed from the raw data directly, and to observe the amount of information lost by the EW preprocessing. The second objective is tp usethe EW-CRB to optimize the EW parameters. More precisely, the optimal value will be chosen as the one minimizing the EW-CRB.

A. CRB for the damped sinusoid estimation problem with colored noise

Consider the data model [START_REF] Lovisolo | Modeling of electric disturbance signals using damped sinusoids via atomic decompositions and its applications[END_REF]. We assume that w(n) is zero-mean stationary Gaussian circular noise with unknown covariance. Define:

Θ d = [α 1 α 2 ... α L ω 1 ω 2 ... ω L ] T , Θ n = [a 1 a 2 ... a L θ 1 θ 2 ... θ L ] T ( 11
)
where Θ d are the desired parameters and Θ n are the nuisance parameters 3 . Define Θ = [Θ T d Θ T n ] T to refer to the desired and nuisance parameters. Assume that the noise follows some parametric model and Θ w represents the parameters of the noise model. Let Ω = [Θ T Θ T w ] T denotes all model parameters. The CRB, which is given by the diagonal terms of the Fisher Information Matrix (FIM) inverse [START_REF] Stoica | Spectral analysis of signals[END_REF], is a lower bound 4 on the variance of the model parameters, i.e., MSE(Ω) ≥ CRB(Ω) = F -1 Ω , where F Ω denotes the FIM for parameter Ω and MSE stands for mean-squared error.

We will utilize the following lemma to reduce CRB computation to F Θ only: Lemma 3.1 ( [START_REF] Francos | Bounds for estimation of complex exponentials in unknown colored noise[END_REF]): The elements of the FIM corresponding to the cross terms of Θ and Θ w are zeros.

The FIM is given elementwise by:

[F Θ ] ij = E [ ∂L(Y ) ∂Θ i ∂L(Y ) H ∂Θ i ] ( 12
)
where Y is the observation vector, i.e Y = [y(0) y(1) . . . y(N -1)] T , and L(Y ) is the log-likelihood function with respect to Y .

For a circular Gaussian process, the FIM can be written as [START_REF] Stoica | Spectral analysis of signals[END_REF]:

[F Θ ] ij = 2 Re { ∂X * ∂Θ i Γ -1 ∂X ∂Θ j } + tr ( Γ -1 ∂Γ ∂Θ i Γ -1 ∂Γ ∂Θ j ) = 2 Re { ∂X * ∂Θ i Γ -1 ∂X ∂Θ j } ( 13
)
where X = [x(0) x(1) ... x(N -1)] T and Γ is the noise covariance. Note that the second term in [START_REF] Ruiz | Parameter estimation of exponentially damped sinusoids using a higher order correlation-based approach[END_REF] is zero since the noise covariance is assumed to be independent from the damped sinusoids parameters.

[F Θ ] ij = 2 Re { ∂X * ∂Θ i Γ -1 ∂X ∂Θ j } . ( 14
)

B. Parameter Selection

The EW method has an important parameter W = e d+jβ to be selected. We use the expressions [START_REF] Hua | The most efficient implementation of the IQML algorithm[END_REF] for that purpose. An intuitive idea consists in enhancing the signal-to-noise ratio after pre-processing.

In the appendix, we derive the asymptotic covariance of the noise, which turns to have a complicated expression depending on the noise auto-covariance function. We compute instead an explicit expression in the white noise case: Lemma 3.2: If w(k) was zero-mean circular white Gaussian noise. Then, the EW noise v(k) is asymptotically circular white Gaussian noise with auto-covariance function:

ρ(k) = σ 4 1 -|W | 2 δ(k)
Thus, we can select the EW parameter by maximizing the averaged post-processing SNR given by:

max d,β (1 -|W 2 |)∥A W ∥ 2 (15)
where

A W = [A W (1)...A W (L)] T .
Note that the noise covariance is independent of β. Therefore, for a fixed d, the maximization problem is interpreted as:

max β ∥A W ∥ 2 . ( 16
)
However, if we have unknown noise, we cannot compute the SNR. Nevertheless, it turns out that selecting β according to [START_REF] Stoica | Subspace-based frequency estimation in the presence of moving-average noise using decimation[END_REF] and choosing d as the slowest damping can serve as a good choice for a wide range of noise statistics. Indeed, the latter choice for d is motivated by the fact that the 'effective life duration' of the exponentially damped signal e αn is approximately N = 1/(1 -e α ). Hence, it is natural to choose a window size corresponding to the desired signal life duration (i.e. choosing d as the slowest damping factor).

As a result, we provide the following iterative algorithm:

1) Initialize d = 0 and β randomly 5 .

2) Compute the cyclo-correlation as in [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise[END_REF] or [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF].

3) Apply a standard method (e.g MP or IQML) to r W to compute the complex frequencies. Compute the complex amplitude via solving a least square problem.

4) Set d as the slowest nonpositive damping and choose β by solving [START_REF] Stoica | Subspace-based frequency estimation in the presence of moving-average noise using decimation[END_REF].

5) If an end criterion6 is reached, then stop. Otherwise, go to step 2.

C. The CRB for the exponentially windowed signal (EW-CRB)

In order to compute the EW-CRB, we need to computed to EW noise covariance, this is done in Appendix I.

IV. OPTIMIZATION OF THE EW PARAMETERS V. SIMULATION

A. Parameter Selection

The data we consider here is as in [START_REF] Papadopoulos | Parameter estimation of exponentially damped sinusoids using higherorder statistics[END_REF], [START_REF] Abed-Meraim | Parameter estimation of exponentially damped sinusoids using second order statistics[END_REF]:

y(n) = e b1n + e b2n + w(n)
where b 1 = -0.1 + 0.52(2π)i, b 2 = -0.2 + 0.42(2π)i. We assume, at first, zero-mean circular white Gaussian noise. We apply the cyclo-correlation [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise[END_REF]. For each experiment, the sample size is set to N = 100. N r = 500 independent Monte-Carlo runs are performed to evaluate the statistics.

Figure 1, shows the plot of the cost function. The optimal value is (d, β) = (-0.55, 0.37). Note that β = 0.37 is optimal for every fixed d.

B. Performance Comparisons

To illustrate the performance of the EW preprocessing method, we consider the MP algorithm (which uses implicitly the SD technique) [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF]. By combining the MP method with HOS [START_REF] Papadopoulos | Parameter estimation of exponentially damped sinusoids using higherorder statistics[END_REF], we have a method to be called HOS-MP method. By applying the MP method to the exponentially window signal, we have the EW-MP method. We use CS-MP to refer to the EW-MP with d = 0. In the sequel, we compare the MP, HOS-MP, EW-MP and CS-MP methods.

The data we consider is same as previous subsection. The performance is measured by the meansquare error (MSE) defined as MSE= The same β is used for the CS-MP method. We set the processing window size for all MP methods to 15 (see [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF] for more details).

The following cases are considered:

Case 1: White Gaussian circular noise: Figure 3 shows the MSE versus SNR. The EW-MP method performs the best for low SNRs even in this case where the noise assumptions required by all methods are satisfied. For medium and high SNR, all methods perform closely.

where the white Gaussian noise has the power σ 2 (SNR is defined as -10 log 10 σ 2 ). Figure V-B shows the MSE versus SNR for a = 0.9. In this case, the EW-MP method performs significantly better for low and medium SNRs. The HOS-MP seems to be highly affected by the small sample size. Note that EW-MP has better performance than CS-MP which illustrates the utility of the damped exponential windowing.

VI. CONCLUSION

In this paper, we discussed a sinusoids estimation method that works acceptably with wide range of stationary noise. The idea is based on computing the cyclo-correlation (or a generalized exponential window) of the signal. We provided a simple algorithm for EW parameter selection and performance optimization. We showed through simulations that it has a superior performance to other methods especially in the case of non-Gaussian and/or colored noise. Comparisons with the CRB showed that it can give close to CRB performance at moderate and high SNRs. Proof: From (8), we have

v(k) = ∞ ∑ n=0 w(n + k)w * (n)W n (19)
First, we compute the mean of v(k). Using [START_REF] Stoica | Spectral analysis of signals[END_REF] we have:

µ k = E[v(k)] = ∞ ∑ n=0 E[w(n + k)w * (n)]W n = ρ(k) 1 -W = ρ(k)L, L = 1 1 -W .
Therefore,

C v (k, ℓ) = E[v(k)v * (ℓ)] -µ k µ * ℓ = E[v(k)v * (ℓ)] -|L| 2 ρ(k)ρ * (ℓ)
Using noise gaussianity and circularity, we compute the first term as:

E[v(k)v * (ℓ)] = E [ ∑ m,
C v (s) = σ 4 ∞ ∑ m=0 ∞ ∑ n=0 W n W * m δ(n -m + s)δ * (n -m) = σ 4 ∞ ∑ n=0 |W | 2n = σ 4 1 -|W | 2 .

  frequency parameter β since we consider all cyclo-frequencies in the range [β 0 , β 1 ]. Indeed, one can observe in Fig. V-A (section IV) that a 'bad' selection of parameter β might lead to poor estimation performance at low or moderate SNRs.

Figure V-A 2 ) A 2 1Fig. 1 .

 221 Figure V-A shows the sensitivity of the method with respect to β. It is observed that the maximum of the cost function corresponds to the best performance obtained. Note also that optimization is highly needed since bad a choice of β causes bad performance.

  Fig. 2. (a) The MSE curves with different values of β with fixed d = -0.1 in the case of white Gaussian noise. (b) MSE curves comparing the performance between the optimization (15) and the proposed algorithm in the case of white Gaussian noise.

Figure

  Figure V-A, shows the result when varying d between the best value suggested by (15) versus slowest damping factor suggested by the proposed algorithm in the case of white Gaussian noise. It can be seen that although d = -0.55 results in slightly improved results for mid and high SNRs, our proposed algorithm gives very reasonable performance. The comparison with the case d = 0 is given next.

√ 1

 1 Nr ∑ Nr r=1 ∥ br -b∥ 2 , where b denotes the parameter vector, i.e., b = [b 1 , b 2 ] T and br is the estimate of b at the r-th run. For the EW-MP method, we chose β = 0.37 and d = -0.1.

Fig. 3 .

 3 Fig. 3. White Gaussian circular noise case.

Case 2 :

 2 Colored Gaussian noise: The colored noise is generated by the AR(a) filter w(n) = aw(n -1) + white Gaussian noise

Figure V-BCase 3 :Case 4 :Fig. 4 .

 344 Figure V-B shows the MSE performance between the three methods versus the value of a with SNR=20 dB. Clearly, the EW preprocessing increases the robustness of MP against AR noise. Case 3: White non-Gaussian noise: The noise is complex with real and imaginary parts uniformly distributed over [-σ/ √ 2, σ/ √ 2]. Figure V-B, shows also the superiority of EW-MP method. Case 4: Colored non-Gaussian noise: In this case, the colored noise is generated filtering the white uniformly distributed by a AR(a) filter, where a = 0.9. As shown in Figure V-B, both the MP and HOS-MP methods fail to compete with the EW-MP method because the noise assumptions required by the former two are not satisfied.

Fig. 5 .W

 5 Fig. 5. (a) The case of white uniformly distributed noise. (b) The case of AR(0.9) uniformly distributed noise.

W 2 :

 2 n W n W * m w(n + k)w * (n)w * (m + ℓ)w(m) n W * m E [w(n + k)w * (n)w * (m + ℓ)w(m)] = ∑ m,n W n W * m (E [w(n + k)w * (n)] E [w * (m + ℓ)w(m)] + E [w(n + k)w * (m + ℓ)] E [w * (n)w(m)] + E [w(n + k)w(m)] E [w * (n)w * (m + ℓ)]) = ∑ m,n W n W * m (ρ(k)ρ * (ℓ) + ρ(n -m + k -ℓ)ρ * (n -m)) = |L| 2 ρ(k)ρ * (ℓ) + ∑ m,n W n W * m ρ(n -m + k -ℓ)ρ * (n -m)and the result follows.We can proof Lemma 3.2 now: Proof of Lemma 3.Set ρ(k) = δ(k). If s ̸ = 0, we have C v (s) = 0. If s = 0 then:

Note that, in contrary to most existing methods, the noise can be colored and non-Gaussian. The only assumption is its wide sense stationarity.

More precisely, we assume that e αmN << 1, ∀m = 1, ..., L.

This paper focuses only on the frequency and damping parameters estimation.

This lower bound applies for all unbiased estimators.

If any prior knowledge is available it can be used for initialization, e.g. if it's known that all sinusoids are damped, then d can be initialized to be less than zero.

We observed that 4-6 iterations yield acceptable performance.