Vinh-Thong Ta

Abderrahim Elmoataz

Olivier Lézoray

Nonlocal Graph Morphology

Keywords: Nonlocal, graphs, morphology

Mathematical morphology (MM) operators can be defined in terms of algebraic sets or as partial differential equations. In this paper, we introduce a novel formulation of MM over weighted graphs of arbitrary topology. The proposed framework recovers local algebraic and PDEs-based formulations of MM and introduces nonlocal configurations. This enables to PDEs-based methods to process any discrete data that can be described by a graph such as high dimensional data defined on irregular domains.

Introduction

Two formulations of morphological operators can be found in literature: the algebraic formulation and the Partial Differential Equations (PDEs) formulation. Morphological algebraic flat dilation δ and erosion ε of a function f 0 :R n →R are usually formulated by: δ B (f 0 )(x) = sup{f 0 (x+y) : y∈B} and ε B (f 0 )(x) = inf{f 0 (x+y) : y∈B} with B a compact convex symmetric set (called structuring element). By using structuring sets B={x : x p ≤1}, the general PDEs generating these flat dilations and erosions [START_REF] Roger | Evolution equations for continuous-scale morphology[END_REF] are as follows:

∂δ(f ) ∂t =∂tf =+ ∇f p and ∂ε(f ) ∂t =∂tf =-∇f p ( 1
)
where f is a modified version of f 0 , ∇ is the gradient operator, . p corresponds to the L p -norm, and one has the initial condition ∂ t=0 f =f 0 . With different values of p, one obtains different structuring elements: a rhombus for p=∞, a disc with p=2 and a square with p=1 [START_REF] Roger | Evolution equations for continuous-scale morphology[END_REF]. Whatever the chosen formulation (algebraic or PDEs), there is actually no satisfying formulation of MM for the processing of multivariate unorganized data. The algebraic formulation has drawbacks for multivariate data (no natural ordering for vectors to form a complete lattice [START_REF] Aptoula | A comparative study on multivariate mathematical morphology. Pattern (a) (b) (c) Figure 2: Morphological Processing of a Mesh Graph (see text for details). Original data sets Dilation Erosion Opening Figure 3: Morphological Procesing of an Image Database Graph (see text for details)[END_REF]) and the PDEs formulation has drawbacks for unorganized data (spatial discretization is difficult). In addition, whatever the formulation, both only consider local interactions on the data while nonlocal schemes have recently received a lot of attention for arbitrary data processing [START_REF] Buades | Nonlocal image and movie denoising[END_REF][START_REF] Elmoataz | Nonlocal discrete regularization an weighted graphs: a framework for image and manifolds processing[END_REF]. In the spirit of [START_REF] Elmoataz | Nonlocal discrete regularization an weighted graphs: a framework for image and manifolds processing[END_REF], we propose to consider a discrete version of continuous MM over weighted graphs that naturally enables nonlocal processing of any multivariate data living on any domain.

Operators on graphs

We consider the general situation where any discrete domain can be viewed as a weighted graph. Let G=(V, E, w) be a weighted graph composed of a finite set of vertices V and a finite set of weighted edges E⊆V ×V . An edge (u, v)∈E connects two adjacent or neighbor vertices u and v of V . The set of neighbors of a vertex u is denoted by N (u)={v∈V \{u} : (u, v)∈E}. The graph can be associated with a weight function w:E→R + . Graphs are assumed to be simple, connected and undirected, implying that the weight function w is symmetric. For the sake of simplicity, w(u, v) will be denoted by w uv and notation v∼u will mean that vertices u and v are adjacent. Let f :V →R be a discrete real-valued function that assigns a real value f (u) to each vertex u∈V . Moreover, each vertex u∈V can be assigned with a feature vector F (f 0 , u)∈R m that is used to compute weights:

g1(u, v)= exp(-(f 0 , u) -F (f 0 , v) 2 /σ 2 ) with σ>0 , g2(u, v)=( F (f 0 , u) -F (f 0 , v) + ǫ) -1 with ǫ>0, ǫ→0 .
A typical feature vector is F (f 0 , u)=f 0 (u) (local processing) or F (f 0 , u) = F τ (f 0 , u) = {f 0 (v) : v∈N τ (u) ∪ {u}} with N τ (u) = {v∈V \{u} : uv ≤τ } (nonlocal processing). We denote A as the set of connected vertices with A⊂V such that for all u∈A, there exists a vertex v∈A with (u, v)∈E. We denote by ∂ + A and ∂ -A, the external and internal boundary sets of A, respectively. For a given vertex u∈V :

∂ + A = {u∈A c : ∃v∈A with (u, v)∈E} and ∂ -A = {u∈A : ∃v∈A c with (u, v)∈E} , (2)
where A c =V \ A is the complement of A.

Weighted morphological differences

All the basic operators considered in this paper are defined from the difference operator or the discrete derivative. In this work, we consider a definition that allows the expressions of discrete weighted gradient on graphs. The weighted difference operator [START_REF] Elmoataz | Nonlocal discrete regularization an weighted graphs: a framework for image and manifolds processing[END_REF] d w :H(V )→H(E) of a function f ∈H(V ) is the vector of all weighted discrete derivatives:

dwf = `(dwf )(u, v)
´(u,v)∈E where for all (u, v)∈E ,

(dwf )(u, v)=w 1/2 uv `f (v)-f (u)
´,

and

∂ v f (u) = (d w f )(u, v)
is the discrete partial derivative of f . This definition has the following properties for a function defined in a Euclidean space:

∂ v f (u)= -∂ u f (v), ∂ u f (u)=0 and if f (u) = f (v) then ∂ v f (u)=0.
Based on the difference operator, we define two new weighted morphological directional difference operators. The weighted morphological external and internal difference operator are respectively:

(d + w f )(u, v)=w 1/2 uv `max `f (u), f (v) ´-f (u) ´and (d - w f )(u, v)=w 1/2 uv `f (u)-min `f (u), f (v) ´´, (3)
with the following properties:

(d - w f )(u, v)=(d + w f )(v, u) (d + w f )(u, v)= max `0, (dwf )(u, v)
´,

(d - w f )(u, v)= -min `0, (dwf )(u, v)
´.

since one has max(a, b)

-a = max(0, a -b), a -min(a, b) = -min(0, b-a), and max(0, a-b) = -min(0, b-a).
The two proposed weighted morphological directional difference operators (3) recover classical directional difference operators on unweighted graphs (w=1) and extend them to weighted graphs that enables more adaptation in the difference computation. Finally, the corresponding external and internal partial derivatives are:

∂ + v f (u)=(d + w f )(u, v) and ∂ - v f (u) = (d - w f )(u, v).

Weighted morphological gradients

The weighted gradient of a function f ∈H(V ) at vertex u is the vector of all edge directional derivatives (∇wf

)(u) = `∂vf (u) ´(u,v)∈E .
Following this definition, we introduce two new weighted morphological (internal and external) gradients based on the internal and external partial derivatives such as

(∇ + w f )(u)= `∂+ v f (u) ´(u,v)∈E and (∇ - w f )(u)= `∂- v f (u) ´(u,v)∈E . (4)
The external gradient of a function is a directional difference operator defined as the difference between an extensive operator and the function (a typical one being the max). Similarly, the internal gradient uses an anti-extensive (the min) operator [START_REF] Rivest | Morphological gradients[END_REF]. In the sequel, we use the L p -norm of gradients (4) for a given vertex u∈V , we have

(∇ + w f )(u) p = h X v∼u w p/2 uv ˛max `0, f (v)-f (u) ´˛p i 1/p and (∇ - w f )(u) p = h X v∼u w p/2 uv ˛min `0, f (v)-f (u) ´˛p i 1/p ; (5)
and for the L ∞ -norm, we have

(∇ + w f )(u) ∞ = max v∼u " w 1/2 uv ˛max `0, f (v)-f (u) ´˛" and (∇ - w f )(u) ∞ = max v∼u " w 1/2 uv ˛min `0, f (v)-f (u) ´˛" . (6)
Same expressions can be obtained for the general weighted gradient. One can note that general definitions presented in this Section are defined on graphs of arbitrary topology. Hence, they can be used to process any discrete regular or irregular data sets that can be represented by a weighted graph. Moreover, local and nonlocal settings are directly handled in these definitions and both are expressed by the graph topology in terms of neighborhood connectivity [START_REF] Elmoataz | Nonlocal discrete regularization an weighted graphs: a framework for image and manifolds processing[END_REF].

Relations with algebraic morphology

The previously defined external and internal gradients operate on any graph structure. In the sequel, we show that in the particular case of an unweighted (w=1) graph and with p=∞, our gradient formulations recover algebraic morphological operators where the structuring element is provided by the graph neighborhood i.e. B=N (u) for all u∈V . The L ∞ -norms (6) of the proposed external and internal gradients ∇ + w and ∇ - w recover the classical definition of algebraic morphological external and internal gradients. Indeed, for the external gradient

(∇ + w f )(u) ∞ = max v∼u `max `0, f (v)-f (u) ´= max v∼u `f (u), f (v) ´-f (u) =δ(f )(u)-f (u) ,
and similarly for the internal one. With these latter relations, we immediately recover the algebraic classical morphological gradient and Laplace operators:

(∇ + w f )(u) ∞ + (∇ - w f )(u) ∞ =δ(f )(u)-ε(f )(u) and (∇ + w f )(u) ∞ -(∇ - w f )(u) ∞ =δ(f )(u)+ε(f )(u)-2f (u) .
Finally, our formulation recovers classical morphological gradient and Laplace operators and extend them to weighted graphs that define new families of local and nonlocal weighted morphological gradients and Laplace operators.

Relations with graph boundary

Intuitively from definitions (2), dilation over A can be interpreted as a growth process that adds vertices from ∂ + A to A. By duality, erosion over A can be interpreted as a contraction process that removes vertices from ∂ -A. The decomposition of f into its level sets is denoted f l =χ(f -l) where χ is the Heaviside function (a step function). Then, one can prove [START_REF] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF] that for any level set f l , at vertex u∈V , the L p -norm of the gradient (∇ w f l )(u) can be decomposed as

(∇wf l )(u) p p = (∇ + w f l )(u) p p + (∇ - w f l )(u) p p .
Moreover, we can also deduce that

(∇wf l )(u) p= ( (∇ + w f l )(u) p if u ∈ ∂ + A l , (∇ - w f l )(u) p if u ∈ ∂ -A l .

Non local graph morphology

Starting from PDEs-based dilation and erosion formulations (1), we define the discrete analogue of such definitions and obtain the following expressions over graphs. For a given initial function

f 0 ∈H(V ): ∂δ(f )(u) ∂t = ∂tf (u) = + (∇ + w f )(u) p and ∂ε(f )(u) ∂t = ∂tf (u) = -(∇ - w f )(u) p ∀u∈V
, with the initial condition ∂ t=0 f =f 0 (f is a modified version of f 0 ) and ∇ + w and ∇ - w are the weighted discrete morphological gradients defined in (4).

Dilation and erosion processes

As for the continuous case, a simple variational definition of dilation applied to f l can be interpreted as maximizing a surface gain proportional to the gradient norm + ∇ w f l p . Similarly, erosion can be viewed as minimizing a surface gain proportional to -∇ w f l p . Dilation of f l on A l corresponds to only consider the external boundary set ∂ + A l and can be expressed by ∂tf l (u)=+ (∇ + w f l )(u) p where the gradient (∇wf l )(u) p is reduced to the external gradient

(∇ + w f l )(u) p.
Similarly, erosion of f l on A l can be defined by ∂tf l (u)=-(∇ - w f l )(u) p. Finally, by extending these two processes for all the levels of f , we can naturally consider the following two families of dilation and erosion. These two processes are parameterized by p and w over any weighted graph G = (V, E, w):

δ(f )(u) = ∂tf (u) = + (∇ + w f )(u) p and ε(f )(u) = ∂tf (u) = -(∇ - w f )(u) p.
To solve these latter dilation and erosion processes, on the contrary to the PDEs case, no spatial discretization is needed thanks to derivatives that are directly expressed in a discrete form. Then, by using discretization in time, and with the usual notation f (u, n) ≈ f (u, n∆t), the general iterative schemes for dilation and erosion, can be defined at time n+1, for all u ∈ V , as

f n+1 (u)=f n (u)+∆t (∇ + w f ) n (u) p and (7)
f n+1 (u)=f n (u)-∆t (∇ - w f ) n (u) p . ( 8
)
The initial condition is f (0) =f 0 where f 0 ∈H(V ) is the initial function defined on the graph vertices. If dilation is considered and with the corresponding gradient norms, [START_REF] Soille | Morphological Image Analysis -Principles and Applications[END_REF] becomes for 0<p<+∞

f n+1 (u) (5) = f n (u)+∆t X v∼u w p/2 uv max 0, β n uv p 1 p (9)
and for p=∞

f n+1 (u) (6) = f n (u)+∆t max v∼u w 1/2 uv max 0, β n uv , (10)
where β n uv =f n (v)-f n (u). At each step of the algorithms, the new value at vertex u only depends on its value at step n and the existing values in its neighborhood. The proposed dilation and erosion formulations operate on any graph structures. It can be shown that with an adapted graph topology and an appropriated weight function, the proposed methodology for dilation and erosion is linked to well-known methods defined in the context of image processing. With an unweighted (w = 1) 4-adjacency grid graph associated with a grayscale image defined as function f 0 :V ⊂Z 2 →R, our formulation of dilation (p = 2) recovers the exact Osher-Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] upwind first order discretization scheme of dilation. The proposed dilation expression also recovers the classical algebraic flat morphological dilation formulation over graphs. In the case where p=∞ with a constant discretization time ∆t=1 and a constant weight function (w = 1), (10) becomes for u∈V :

f n+1 (u)=f n (u)+ max v∼u " max `0, β n uv ´" = max v∼u `f n (v), f n (u) ´. (11)
In that case, the structuring element is provided by both the graph topology and the vertices' neighborhoods. For instance, if we consider a 8-adjacency image grid graph, it is equivalent to a dilation by a square structuring element of size 3×3. It is important to note that for the case of weighted graphs, the proposed morphological operators are operators that adapt their behavior according to the importance of weights with their neighbors.

Results

Figure 1 presents dilation and closing of an original scalar grayscale image considered as a function f 0 : V ⊂ Z 2 →R that defines a mapping from vertices to grayscale pixel values. This example shows the adaptative aspect of our morphological framework. First row of Fig. 1 shows results for the case where p=∞, ∆t=1. This corresponds to the algebraic processing. In that case, the graph is a 25-adjacency grid graph (equivalent to a circle structuring element of radius 2, denoted as G 2 ). The last three next rows show results for the case where p=2. Second row of Fig. 1 presents results for an unweighted (w=1) 4-adjacency grid graph (denoted as G 0 . This case corresponds to PDEs-based morphological processing. Difference between second and third row is that in the latter case, we use a weighted (w=g 1 ) graph instead of an unweighted one. These results show the benefits of weights that enable to better preserve edge information as compared to the unweighted cases. Fourth row presents nonlocal patchbased results. In this example the graph is a 25-adjacency grid graph weighted by function g 1 associated with patches of size 5×5 as pixel features. These results clearly show that this morphological processing better preserves image components such as edges, fine structures and repetitive elements. Figure 2 shows application of morphological alternated sequential filters on 3D data (a mesh, see Fig. 2(a)).

The graph corresponds here to the natural graph representation of meshes where each vertex is associated to the 3D spatial coordinates of each mesh point. This latter graph is weighted by function g This database consists in grayscale images scanned from digits 0 to 9. Each image is of size 16×16. In these experiments, we use three randomly subsampled test sets of 100 samples for each database. Graphs associated with the original data sets are k-NN graphs weighted with function g 2 . It is important to note that each vertex of the graphs corresponds to an image sample and is described by a 256-dimensions (R 16×16 ) feature vector. Dilation, erosion and opening operations are respectively shown in Figure 3. One can note the filtering and simplification effects of opening operation on the data. Indeed, this operation tends to reduce the data to new artificial and uniform image samples. Finally, such morphological operations on databases can be used as pre-processing steps for data mining purposes.

Conclusion

In this paper, a novel formalism for nonlocal morphology on graphs has been proposed. This provides a framework that extends PDEs-based methods to discrete local and nonlocal schemes. Moreover, this enables to process by morphological means any high dimensional unorganized multivariate data represented by a weighted graph. The potentiality and the flexibility of our approach has been illustrated for the morphological processing of organized data (regular domain) and unorganized data (irregular domain).

1 .

 1 Figure1presents dilation and closing of an original scalar grayscale image considered as a function f 0 : V ⊂ Z 2 →R that defines a mapping from vertices to grayscale pixel values. This example shows the adaptative aspect of our morphological framework. First row of Fig.1shows results for the case where p=∞, ∆t=1. This corresponds to the algebraic processing. In that case, the graph is a 25-adjacency grid graph (equivalent to a circle structuring element of radius 2, denoted as G 2 ). The last three next rows show results for the case where p=2. Second row of Fig.1presents results for an unweighted (w=1) 4-adjacency grid graph (denoted as G 0 . This case corresponds to PDEs-based morphological processing. Difference between second and third row is that in the latter case, we use a weighted (w=g 1 ) graph instead of an unweighted one. These results show the benefits of weights that enable to better preserve edge information as compared to the unweighted cases. Fourth row presents nonlocal patchbased results. In this example the graph is a 25-adjacency grid graph weighted by function g 1 associated with patches of size 5×5 as pixel features. These results clearly show that this morphological processing better preserves image components such as edges, fine structures and repetitive elements. Figure2shows application of morphological alternated sequential filters on 3D data (a mesh, see Fig.2(a)). The graph corresponds here to the natural graph representation of meshes where each vertex is associated to the 3D spatial coordinates of each mesh point. This latter graph is weighted by function g 1 . Figures2(b) and 2(c) show the evolution of the simplification for iterations n=10 and n=50, respectively. It is important to note that the number of vertices does not change during the filtering process and vertices have moved to similar spatial coordinates. Figure3shows examples of real world databases processing within our morphological framework. The first row of this figure shows the original data sets. Images come from the United States Postal Service (USPS) handwritten database. This database consists in grayscale images scanned from digits 0 to 9. Each image is of size 16×16. In these experiments, we use three randomly subsampled test sets of 100 samples for each database. Graphs associated with the original data sets are k-NN graphs weighted with function g 2 . It is important to note that each vertex of the graphs corresponds to an image sample and is described by a 256-dimensions (R 16×16 ) feature vector. Dilation, erosion and opening operations are respectively shown in Figure3. One can note the filtering and simplification effects of opening operation on the data. Indeed, this operation tends to reduce the data to new artificial and uniform image samples. Finally, such morphological operations on

Figure 1 :

 1 Figure 1: Morphological Image Processing with graphs (see text for details).


Wenting Li

François Auger

Zhuosheng Zhang
email: zszhang@mail.xjtu.edu.cn

Xiangxiang Zhu

Newton time-extracting wavelet transform: an effective tool for characterizing frequency-varying signals with weakly-separated components and theoretical analysis

Keywords: Time-frequency analysis, Time reassignment operator, Group delay, Ridge, Newton time-extracting wavelet transform, Signal reconstruction

In this paper, we propose a high resolution time-frequency (TF) analysis method called Newton time-extracting wavelet transform (NTEWT), which is designed to analyze frequency-varying signals with fast varying group delay (GD). Firstly, we discuss the relationship among the fixed points of time reassignment operator, the ridge points of wavelet transform and GD of the signal. Combining the above relations and Newton algorithm, we propose a Newton GD estimator. By only retaining the TF information most related to frequency-varying features of the signal and removing the weakly-related TF coefficients, we further introduce the NTEWT, which can not only achieve a more concentrated TF representation, but also enable signal reconstruction. Meanwhile, we develop a theoretical analysis of NTEWT under the mathematical framework. Firstly, we introduce a precise mathematical definition of a class of weakly-separated frequency-varying chirp-like components, and we prove that Newton GD estimator can accurately estimate GD of arbitrary function in this class, and NTEWT does indeed succeed in decomposing these functions. Finally, we use numerical experiments to evaluate the performance of the proposed NTEWT in terms of TF concentration, GD estimation and signal reconstruction.

Introduction

Non-stationary signals widely appear in various fields, such as radar and sonar [START_REF] Cohen | Time-frequency analysis: theory and applications[END_REF][START_REF] Amin | Radar signal processing for elderly fall detection: the future for in-home monitoring[END_REF], seismic [START_REF] Fomel | Seismic data decomposition into spectral components using regularized nonstationary autoregression[END_REF], biomedicine [START_REF] Park | Time-frequency analysis of EEG asymmetry using bivariate empirical mode decomposition[END_REF][START_REF] Yang | Synchrosqueezed wave packet transforms and diffeomorphism based spectral analysis for 1D general mode decompositions[END_REF], and mechanical engineering [START_REF] Zhu | Two robust approaches to multicomponent signal reconstruction from STFT ridges[END_REF][START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF], etc. Time-frequency analysis (TFA) can convert onedimensional signals into two-dimensional time-frequency (TF) representations (TFR) efficiently, which helps to detect the TF feature of non-stationary signals in the TF plane visually [START_REF] Mallat | A wavelet tour of signal processing: the sparse way (3rd edn)[END_REF][START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF]. Many effective TFA methods mainly include linear TFA methods, nonlinear TFA methods and bilinear TFA methods.

Linear TFA methods, such as the short-time Fourier transform (STFT) [START_REF] Mallat | A wavelet tour of signal processing: the sparse way (3rd edn)[END_REF] and the wavelet transform (WT) [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF], calculate the inner product between the signal and the basis function to characterize the TF structure of signal. Generally, the better the basis function matches the signal, the more suitable it is to characterize TF features with concentrated energy. Nonlinear TFA methods aim to construct new basis functions to characterize the non-linearity of signals, such as linear chirplet transform [START_REF] Mann | The chirplet transform: physical considerations[END_REF], local polynomial transform [START_REF] Li | Local polynomial Fourier transform: a review on recent developments and applications[END_REF], and matched demodulation transform [START_REF] Wang | Matching demodulation transform and synchrosqueezing in time-frequency analysis[END_REF], etc. However, restricted by the Heisenberg uncertainty principle, linear and nonlinear TFA methods often generate blurry TFRs, and they had difficulty providing accurate TF descriptions for signals. Besides, bilinear TFA methods, such as the Wigner-Ville distribution [START_REF] Cohen | Time-frequency analysis: theory and applications[END_REF] and its variants, which can effectively improve the TF resolution of mono-component signals, but introduce interference terms for multi-component signals and thus reduce the readability of the TFRs. In order to overcome these shortcomings, some advanced methods have been proposed successively.

The reassignment method (RM), proposed by Kodera et al. to the spectrogram [START_REF] Kodera | A new method for the numerical analysis of nonstationary signals[END_REF][START_REF] Kodera | Analysis of time-varying signals with small BT values[END_REF], and generalized by Auger and Flandrin to any bilinear time-frequency or time-scale distribution [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF], is an old but under-explored method [START_REF] Flandrin | Time-frequency reassignment: from principles to algorithms[END_REF][START_REF] Chassande-Mottin | Time frequency/time-scale reassignment[END_REF][START_REF] Chassande-Mottin | Differential reassignment[END_REF]. In parallel, Daubechies and Maes proposed the synchrosqueezing wavelet transform (SWT) [START_REF] Daubechies | A nonlinear squeezing of the continuous wavelet transform[END_REF] based on a phase technique to analyse auditory signals, which is a special case of the RM. And Huang et al. introduced an empirical mode decomposition (EMD) [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] to decompose the signal into a series of intrinsic mode functions in a data-driven manner. Though attractive because of its effectiveness and simplicity [START_REF] Dragomiretskiy | Variational mode decomposition[END_REF][START_REF] Gilles | Empirical wavelet transform[END_REF], EMD lacks solid mathematical foundations. In this context, SWT has resurfaced as an alternative theoretical way to understand the principle of EMD with a convenient mathematical framework [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF]. RM has, in turn, received new attention [START_REF] Auger | Timefrequency reassignment and synchrosqueezing: an overview[END_REF][START_REF] Auger | On phase-magnitude relationships in the short-time Fourier transform[END_REF][START_REF] Brynolfsson | Parameter estimation of oscillating Gaussian functions using the scaled reassigned spectrogram[END_REF]. Although the RM relocates TF energy aiming to sharpen the TFR, and employs instantaneous frequency (IF) and group delay (GD) estimators, it is not invertible unfortunately. While SWT can reconstruct the interesting signals but suffers from blurred TFRs for signals with fast varying IF. Furthermore, new extensions of SWT have been developed [START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF][START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF][START_REF] Fourer | Chirp rate and instantaneous frequency estimation: application to recursive vertical synchrosqueezing[END_REF][START_REF] Yu | Local maximum synchrosqueezing transform: an energyconcentrated time-frequency analysis tool[END_REF][START_REF] Yu | Multi-synchrosqueezing transform[END_REF][START_REF] Li | Adaptive synchrosqueezing transform with a time-varying parameter for nonstationary signal separation[END_REF][START_REF] Yu | Synchroextracting transform[END_REF], among these improved methods, recently synchroextracting transform (SET) [START_REF] Yu | Synchroextracting transform[END_REF] has attracted great attention. Differing from the squeezing manner of SWT, SET only utilizes a small number of TF coefficients, but achieves a more concentrated TFR. Similar to SWT, SET also develops the corresponding higher-order extensions [START_REF] Bao | Second-order synchroextracting transform with application to fault diagnosis[END_REF][START_REF] Li | Self-matched extracting wavelet transform and signal reconstruction[END_REF][START_REF] Zhu | Three-dimension extracting transform[END_REF]. On the other hand, according to different direction of post-processing, TFA post-processing methods are divided into three categories. The first class considers both time and frequency directions, such as reassigned scalogram (RS), although RS can improve both the time resolution and the frequency/scale resolution of the TFR, unfortunately, it cannot reconstruct the signal. The second class considers only frequency direction, such as SWT [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF], matching synchrosqueezing wavelet transform (MSWT) [START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF], SET [START_REF] Yu | Synchroextracting transform[END_REF] and self-matched extracting wavelet transform (SMEWT) [START_REF] Li | Self-matched extracting wavelet transform and signal reconstruction[END_REF]. These TFA methods can reconstruct signals, and they are more suitable for characterizing slowly or strongly time-varying signals, but they may be no longer applicable for characterizing frequency-varying signals. In order to accurately analyze frequency varying signals, the third class of TFA post-processing along the time direction have emerged, including the time-reassigned synchrosqueezing transform (TSST) [START_REF] He | Time-reassigned synchrosqueezing transform: the algorithm and its applications in mechanical signal processing[END_REF] and transient-extracting transform (TET) [START_REF] Yu | A concentrated time-frequency analysis tool for bearing fault diagnosis[END_REF], and our recently proposed time-reassigned synchrosqueezing wavelet transform (TSWT) [START_REF] Li | Theoretical analysis of time-reassigned synchrosqueezing wavelet transform[END_REF] and time-extracting wavelet transform (TEWT) [START_REF] Li | Time-extracting wavelet transform for characterizing impulsive-like signals and theoretical analysis[END_REF], where we provide rigorous theoretical analysis of TSWT and TEWT for frequency-varying signals under the mathematical framework. Moreover, secondorder TSST (TSST2) [START_REF] Fourer | Second-order time-reassigned synchrosqueezing transform: application to draupner wave analysis[END_REF][START_REF] He | Gaussian-modulated linear group delay model-application to second-order time-reassigned synchrosqueezing transform[END_REF] and second-order TET [START_REF] He | Second-order transient-extracting transform with application to time-frequency filtering[END_REF][START_REF] Yu | Second-order transient-extracting transform for the analysis of impulsive-like signals[END_REF] have been proposed. Indeed, these TFA methods reassigning or extracting along the time direction can reconstruct signals and are suitable for characterizing slowly or strongly frequency-varying signals.

It can be concluded that one of the advantages of TET/TEWT is to extract the TF coefficients of the fixed points of time reassignment operator, because the time reassignment operator can accurately estimate the GD of impulsive signal, and the fixed points of time reassignment operator are equivalent to ridge points of WT, where ridge is the most relevant to the TF features of signals. Therefore, for more general signals, the relationship among the fixed points of the time reassignment operator, ridge points of the WT and GD of signals is particularly important. In this paper, we discuss the relationship among these three, and develop a new GD estimator for strongly frequency-varying signal, which can precisely describe the ridge of WT and accurately estimate the GD of signal. Further we introduce a new frequency-varying TFA method and develop the rigorous theoretical analysis. The main innovations of this paper are: i) We convert the GD estimation problem into solving the fixed point problem and propose a new GD estimator. We firstly develop a strict equivalence condition between the fixed points of the time reassignment operator and ridge points of WT, and further discuss the relationship among the fixed points of the time reassignment operator, ridge points of the WT and GD of a frequency-varying linear chirp signal. Then, combining these relations with Newton algorithm, we propose a Newton GD estimator, which not only accurately estimates the GD of the signal, but also accurately describes the ridge of WT.

ii) Based on the proposed Newton GD estimator, we propose a Newton time-extracting wavelet transform (NTEWT), which can generate a highly concentrated TFR and reconstruct signal. The main idea of NTEWT is to retain the TF information that are most related to the frequency-varying features of the signal and to remove the weakly-related TF coefficients.

iii) We define a class of weakly-separated frequency-varying chirp-like components (WFCC), which relaxes the strict well-separated condition between components to some extent. The function class we define includes possible interference components and is applicable to a wider range of practical signals.

iv) For such function class, we provide a theoretical analysis for NTEWT under a strict mathematical framework, and prove that the Newton GD estimator indeed successfully estimates the GD and the NTEWT reconstructs the signal with high accuracy.

The remainder of this paper is organized as follows. In section 2, we discuss the relationship among the fixed points of the time reassignment operator, ridge points of the WT and GD of the signal. Then in section 3, we combine these relations and Newton algorithm to propose a Newton GD estimator. And further, we introduce the NTEWT. In section 4, we analyze the performance of the NTEWT theoretically, including a theorem and its proof about the NTEWT for a class of WFCC. Section 5 is devoted to the comparative study of the NTEWT and other TFA methods on simulated and real sigals with three quantitative indicators. Finally, conclusions are drawn in section 6.

Time reassignment operator

In this section, we first recall several basic definitions, then study the relationship among the fixed points of the time reassignment operator, ridge points of the WT and GD of the signal.

Basic Notations and Definitions

Definition 1. [START_REF] Mallat | A wavelet tour of signal processing: the sparse way (3rd edn)[END_REF] The Fourier transform (FT) of a given signal x(t) ∈ L 1 (R), i.e. its correlation with a sinusoidal wave e iωt , is defined as:

X(ω) = R x(t)e -iωt dt, (1)
and its inverse FT is defined by [START_REF] Mallat | A wavelet tour of signal processing: the sparse way (3rd edn)[END_REF]:

x(t) = 1 2π R X(ω)e iωt dω. (2)
Definition 2.

[41] The frequency-varying signal is defined in the frequency domain

X(ω) = A(ω)e -iφ(ω) , (3)
where A(ω) is the signal amplitude, φ(ω) and its derivative φ (ω) are the phase and the GD, respectively.

Definition 3. [8] Consider an admissible wavelet ψ(t) ∈ L 2 (R), its FT satisfies C ψ = R + | ψ(ω)| 2 ω dω <
+∞, and then the WT of a signal x(t) ∈ L 2 (R) is defined as a complex-valued function of time b and scale a > 0 by

W ψ x (b, a) = 1 a R x(t)ψ * ( t -b a )dt, ( 4
)
where z * is the complex conjugate of z. And W ψ x (b, a) can be expressed in the frequency domain by Plancherel's theorem

W ψ X (b, a) = 1 2π R X(ω) ψ * (aω)e iωb dω. (5)
where X(ω) and ψ(ω) are FTs of x(t) and ψ(t), respectively.

Denote W ψ X (b, a) = M ψ X (b, a)e iΦ ψ X (b,a) , M ψ X (b, a)
and Φ ψ X (b, a) stand for magnitude and phase of WT, respectively.

Definition 4. [44]

A ridge point of WT W ψ X (b, a) is a time-scale pair (b, a) satisfying the two conditions:

∂ b ln M ψ X (b, a) = 0, ∂ 2 ∂b 2 ln M ψ X (b, a) < 0. ( 6
)
where ∂ b , ∂ 2 ∂b 2 are the first and second order partial derivatives with respect to b, respectively. These conditions state that for each fixed scale a, a time ridge point b corresponds to the time at which a local maximum of magnitude M ψ X (b, a) occurs.

Definition 5. [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF] The complex time reassignment operator tx (b, a) and the time reassignment operator tx (b, a) can be defined respectively by

     tx (b, a) = b + a W tψ x (b,a) W ψ x (b,a) tx (b, a) = b + a { W tψ x (b,a) W ψ x (b,a) }, (7)
where W ψ x (b, a) = 0 and {z} denotes the real part of z. Definition 6. [START_REF] Smart | Fixed point theorems[END_REF] Metric space (Λ, ρ) is a set Λ together with a metric ρ on the set Λ, we say that a mapping F on (Λ, ρ) is a contraction mapping from Λ to Λ if there exists a costant k such that 0 < k < 1 and ρ(F (x), F (y)) ≤ kρ(x, y), ∀x, y ∈ Λ. (8)

Time reassignment operator and ridge

In this subsection, we discuss the relationship between the fixed points of the time reassignment operator and ridge points of WT for frequency-varying signals, and we further provide the equivalent conditions between the fixed points and ridge points.

Theorem 1. For any signal X(ω) = A(ω)e -iφ(ω) , consider a Morlet wavelet expressed by

ψ(ω) = (4πσ 2 ) 1 4 e -σ 2 (ω-ω ψ ) 2 2 , (9)
then for any (b, a) ∈ R × R + , the following statements hold: --→ i ψ and its derivative ψ (ω) = -σ 2 (ω -ω ψ ) ψ(ω), it is easy to derive that by Eq. ( 7)

tX (b, a) = b + a { W ψ X (b, a) W ψ X (b, a) } = b -aσ 2 { W ω ψ X (b, a) W ψ X (b, a) }, ( 10
)
where {z} is the imaginary part of z. On the other hand, noting that

∂ b W ψ X (b, a) = i a W ω ψ X (b, a), it follows that {∂ b ln W ψ X (b, a)} = -{ W ω ψ X (b, a) aW ψ X (b, a) }. (11)
According to Eqs. [START_REF] Mann | The chirplet transform: physical considerations[END_REF] and [START_REF] Li | Local polynomial Fourier transform: a review on recent developments and applications[END_REF], and noticing that ∂

b ln M ψ X (b, a) = {∂ b ln W ψ X (b, a)}, we get tX (b, a) = b + a 2 σ 2 ∂ b ln M ψ X (b, a). (12)
Then, Theorem 1 (a) holds.

|F (b n+1 ) -F (b n )| ≤ k • |b n+1 -b n |, n = 1, 2, • • • (13)
Thanks to the regularity of the WT, the function F (b) is continuously differentiable, then this follows that there exists a ξ n being between b n and b n+1 , such that

|F (ξ n )| ≤ k, n = 1, 2, • • • (14)
Due to lim Theorem 1 well reveals the existence of the fixed points of the time reassignment operator, explains the relationship between the fixed points of the time reassignment operator and ridge points of WT, and provides equivalent condition of the two, i.e. tx (b, a) is a contraction operator. Since the fixed points of the contraction operator exist and are unique, the ridge points exist and are unique. In the next subsection, we will use a frequency-varying linear chirp signal to verify Theorem 1.

n→+∞ b n = b * , it deduces that |F (b * )| ≤ k, that is |∂ b tX (b * , a)| ≤ k < 1. Assume (b,

Time reassignment operator and GD

For a frequency-varying linear chirp signal X(ω), the following lemma shows that tX (b, a) is a contraction operator and further gives the relations among the fixed points of tX (b, a), ridge points of WT and GD of X(ω).

Lemma 1. For a frequency-varying linear chirp signal,

X(ω) = Ae -iφ(ω) , (15)
where φ(ω) = α 2 ω 2 + βω + γ, consider the Morlet wavelet ψ(ω) = (4πσ 2 )

1 4 e -σ 2 (ω-ω ψ ) 2 2
. Then, for any given scale a > 0, the following hold: Proof. (a) Since ω ψ is the frequency center of the Morlet wavelet ψ(ω), the frequency ω and scale a are related by

ω = ω ψ a . (16)
Thus, the GD of signal X(ω) can be expressed as

φ ( ω ψ a ) = α ω ψ a + β. (17)
A straightforward computation shows that the WT W ψ X (b, a) and W ψ X (b, a) of X(ω) can be expressed as

W ψ X (b, a) = Ae -i( α 2 ( ω ψ a ) 2 +β ω ψ a +γ e i ω ψ a b σ(a 2 σ 2 -iα) π 1 2 (a 4 σ 4 + α 2 ) e - a 2 (a 2 σ 2 -iα) 2(a 4 σ 4 +α 2 ) ( α ω ψ a +β-b a ) 2 = 1 a X( ω ψ a )C x (b, a), W ψ X (b, a) = iσ 2 (a 2 σ 2 -iα) a 4 σ 4 + α 2 (α ω ψ a + β -b)e i ω ψ a b X( ω ψ a )C x (b, a), (18)
where

C x (b, a) = a 2 σ(a 2 σ 2 -iα) π 1 2 (a 4 σ 4 +α 2 ) e - a 2 (a 2 σ 2 -iα) 2(a 4 σ 4 +α 2 ) ( α ω ψ a +β-b a ) 2 .
Therefore, by the definition of complex time reassignment operator tX (b, a) in Eq. ( 7) and Eq. ( 18), we have

tX (b, a) = b -ia W ψ X (b, a) W ψ X (b, a) = b - a 2 σ 2 (a 2 σ 2 -iα) a 4 σ 4 + α 2 (b -α ω ψ a -β). (19)
It follows that the time reassignment operator tX (b, a) = { tX (b, a)} can be expressed as

tX (b, a) = b - a 4 σ 4 a 4 σ 4 + α 2 (b -α ω ψ a -β) = α 2 b + a 4 σ 4 (α ω ψ a + β) a 4 σ 4 + α 2 . ( 20
)
Thus, for any given a, we have

0 ≤ ∂ b tX (b, a) = α 2 a 4 σ 4 + α 2 < 1, (21)
that is, tX (b, a) is a contraction operator.

(b) For any given scale a, Eqs. ( 17) and [START_REF] Huang | The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] shows that tX (b, a) is a biased GD estimator.

(c) According to Eq. ( 20), the fixed point curve of tX (b, a) can be expressed as

b = α ω ψ a + β, a > 0. ( 22
)
Hence, by Eq. ( 22) and the GD expression of signal X(ω) in Eq. ( 17), one can get immediately that Lemma 1 (c) holds.

(d) The derivative of ln W ψ X (b, a) with respect to time b is calculated as

∂ b ln W ψ X (b, a) = i ω ψ a - a 2 σ 2 -iα a 4 σ 4 + α 2 (b -α ω ψ a -β). (23)
This, together with the fact

∂ b ln M ψ X (b, a) = {∂ b ln W ψ X (b, a)}, implies ∂ b ln M ψ X (b, a) = - a 2 σ 2 a 4 σ 4 + α 2 (b -α ω ψ a -β), (24)
and

∂ 2 ∂b 2 ln M ψ X (b, a) < 0. (25)
Thus, the ridge trajectory of WT W ψ X (b, a) can be written as

b = α ω ψ a + β, a > 0. ( 26
)
This, together with Eq. ( 22), lead to the conclusion that Lemma 1 (d) holds.

It can be concluded from Lemma 1 that for the frequency-varying linear chirp signal X(ω):

Newton time-extracting wavelet transform

Newton group delay estimator

As can be seen from the previous section, for the frequency-varying linear chirp signal X(ω), the time reassignment operator tX (b, a) can't estimate the GD of X(ω) accurately, but the fixed points of the time reassignment operator can. Indeed, the fixed points of the time reassignment operator can not only accurately estimate the GD of the signal, but also precisely describe the ridge of WT. Therefore, the GD estimation problem is converted into the problem for solving the fixed points of the time reassignment operator, then we can try to employ some root finding algorithms to derive a new GD estimator. Inspired by the relations between Eqs. ( 17), ( 19) and ( 22), we use Newton algorithm to estimate GD in the following.

Lemma 2. For a frequency-varying linear chirp signal,

X(ω) = Ae -iφ(ω) , (27)
where

φ(ω) = α 2 ω 2 + βω + γ, consider the Morlet wavelet ψ(ω) = (4πσ 2 ) 1 4 e -σ 2 (ω-ω ψ ) 2 2 . Let tX (b, a) = b - b -tX (b, a) 1 -∂ b tX (b, a) . ( 28
)
Then the following hold: 19), we immediately get

∂ b tX (b, a) = 1 - a 2 σ 2 (a 2 σ 2 -iα) a 4 σ 4 + α 2 , ( 29
)
it means that we have 1 -∂ b tX (b, a) = 0. Thus, by Eq. ( 28), tX (b, a) = b if and only if tX (b, a) = b.

(b) Combining Eq. ( 19) with Eq. ( 29), we have

tX (b, a) = b -(1 -∂ b tX (b, a))(b -α ω ψ a -β).
Thus, for any (b, a

) ∈ R × R + , we have tX (b, a) = b - b -tX (b, a) 1 -∂ b tX (b, a) = α ω ψ a + β = φ ( ω ψ a ). ( 30
)
that is, tX (b, a) is an accurate GD estimation of the signal X(ω).

(c) According to Eq. ( 30), one has tX (b, a) = b if and only if the point (b, a) is at the GD trajectory

of signal X(ω), i.e. b = α ω ψ a + β. (31)
(d) Combining Eq. ( 26) with Eq. ( 30), we have tX (b, a) = b if and only if the point (b, a) is a ridge point of WT.

Eq. ( 28) can be considered as the first iteration of the Newton algorithm for solving the fixed point equation tX (s, a) = s, in which a sequence defined in the scalar case as

s n+1 = s n -f (sn) f (sn)
, is hoped to converge to a root of the equation f (s) = 0, where f (s) = s -tX (s, a), with s 0 = b. Thus, we have the following definition: Definition 7. A Newton GD (NGD) estimator tX (b, a) can be defined by

tX (b, a) = b - b -tX (b, a) 1 -∂ b tX (b, a) . (32)
It is noteworthy that we convert GD estimation problem into solving the fixed point problem, and propose the NGD estimator tX (b, a) combined with Newton algorithm. It is a new attempt to estimate GD, which also allows to design other GD estimators from other root finding algorithms.

It can be found that various IF estimators, presented in [START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF], [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF], [START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF], [START_REF] Fourer | Chirp rate and instantaneous frequency estimation: application to recursive vertical synchrosqueezing[END_REF], [START_REF] Yu | Synchroextracting transform[END_REF], [START_REF] Li | Self-matched extracting wavelet transform and signal reconstruction[END_REF]

NTEWT

In this subsection, we extend TEWT to the second-order case based on the Newton GD estimator, and propose a new TFA post-processing method extracting along the time direction. Definition 8. The Newton time-extracting wavelet transform (NTEWT) is defined as

N T e(b, a) = W ψ X (b, a) • δ(b -tX (b, a)). (33)
It is well known that SWT [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF] is designed for local harmonic signal in the time domain, MSWT [START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF] is designed for local constant amplitude-modulated linear chirp signal in the time domain, and SMEWT [START_REF] Li | Self-matched extracting wavelet transform and signal reconstruction[END_REF] is designed for local Gaussian amplitude-modulated linear chirp in the time domain. Thus, these TFA methods reassigning or extracting along the frequency direction are suitable for dealing with slowly or strongly varying signals in the time domain. However, the TSST [START_REF] He | Time-reassigned synchrosqueezing transform: the algorithm and its applications in mechanical signal processing[END_REF], TSWT [START_REF] Li | Theoretical analysis of time-reassigned synchrosqueezing wavelet transform[END_REF] and TEWT [START_REF] Li | Time-extracting wavelet transform for characterizing impulsive-like signals and theoretical analysis[END_REF] are suitable for impulsive-like signals, and the proposed NTEWT is designed for frequency-varying linear chirp signal. Therefore, these TFA methods reassigning or extracting along the time direction are suitable for characterizing slowly or strongly varying signals in the frequency domain.

Differing from the squeezing manner of synchrosqueezing methods [START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF], [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF], [START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations[END_REF], [START_REF] He | Time-reassigned synchrosqueezing transform: the algorithm and its applications in mechanical signal processing[END_REF], [START_REF] Li | Theoretical analysis of time-reassigned synchrosqueezing wavelet transform[END_REF], [START_REF] Fourer | Second-order time-reassigned synchrosqueezing transform: application to draupner wave analysis[END_REF], the (3) According to Eq. ( 18), NTEWT maintains the reconstructed property of WT, and signal X(ω) can be recovered by WT along the GD trajectory, i.e.,

X( ω ψ a ) = aW ψ x (φ ( ω ψ a ), a)e -i ω ψ a b C x (φ ( ω ψ a , a) . ( 34
)
It motivates us to recover signal X(ω) from the NTEWT with a similar reconstruction expression:

X( ω ψ a ) = aN T e(φ ( ω ψ a ), a)e -i ω ψ a b C x (φ ( ω ψ a , a) . ( 35
)

Theoretical analysis of Newton time-extracting wavelet transform

In this section, we will provide a rigorous theoretical analysis of NTEWT that can identify and

d > 0, if each element X(ω) = K k=1 X k (ω) = K k=1 A k (ω)e -iφ k (ω) ∈ L ∞ (R), with A k (ω) ∈ C 1 (R) ∩ L ∞ (R) and φ k (ω) ∈ C 2 (R) satisfying: m k = sup ω∈R + |A k (ω)|, m k = sup ω∈R + |φ k (ω)|, |A k (ω)|, |φ k (ω)| ≤ |φ k (ω)|, ∀ω ∈ R + , φ k+1 (ω) -φ k (ω) ≥ d ω , ∀ω ∈ R + , k ∈ {1, 2, ..., K -1}. ( 36
)
Intuitively, for the function class H ,d , each component X k (ω) can be viewed as approximately a harmonic signal locally in the frequency domain, with slowly varying amplitude A k (ω) and GD φ k (ω).

In Theorem 3.1 for the TSWT [START_REF] Li | Theoretical analysis of time-reassigned synchrosqueezing wavelet transform[END_REF], it has been shown that X k (ω) and X k+1 (ω) satisfying well-separated

conditions ( φ k+1 (ω) -φ k (ω) ≥ d ω with d > 2ω ψ ∆) are separated. It has been also proven that the energy of WT W ψ X k (b, a) of the kth component X k (ω) is concentrated in the zone Z k := {(b, a); | φ k ( ω ψ a )-b a | < ∆},
X(ω) = K k=1 X k (ω) = K k=1
A k (ω)e -iφ k (ω) satisfies the following four conditions:

(1) Regularity conditions:

A k (ω) ∈ L ∞ (R) ∩ C 2 (R), φ k (ω) ∈ C 3 (R), k ∈ {1, 2, ..., K}.
(2) Boundedness conditions:

0 < inf ω∈R A k (ω) ≤ sup ω∈R A k (ω) < +∞, k ∈ {1, 2, ..., K} 0 < inf ω∈R φ k (ω) ≤ sup ω∈R φ k (ω) < +∞, sup ω∈R |φ k (ω)| < +∞.
(3) Growth conditions:

|A k (ω)|, |φ k (ω)| ≤ |φ k (ω)|, k ∈ {1, 2, ..., K}.
(4) Weakly-separated conditions: X k (ω) are separated with the distance d k,k+1 > 0, i.e. Our main result is then the following:

φ k+1 (ω) -φ k (ω) ≥ d k,k+1 ω , ω ∈ R + , k ∈ {1, 2, ..., K -1}. ( 37
Theorem 2. Let X(ω) = K k=1 X k (ω) ∈ C ,
d k,k+1 > max{ω ψ ∆ k , ω ψ ∆ k+1 }, ψ φ (t) ≤ G 0 , ψ φ (t) ≤ G 1 , ψ φ (t) ≤ G 2 , with G 0 , G 1 , G 2 being some
constants. Consider the WT W ψ X (b, a) of X(ω) with respect to ψ(ω), as well as the NTEWT N T e(b, a).

For k ∈ {1, 2, ..., K -1}, define a set T k := {(b, a); Π k (a) < b < Π k (a)}, with Π k (a) = max{φ k-1 ( ω ψ a ) + a∆ k-1 , φ k ( ω ψ a ) -a∆ k }, Π 1 (a) = φ 1 ( ω ψ a ) -∆ 1 , Π k (a) = min{φ k ( ω ψ a ) + a∆ k , φ k+1 ( ω ψ a ) -a∆ k+1 }, Π K (a) = φ K ( ω ψ a ) + ∆ K , (38)
Then, provided is sufficiently small, the following hold:

| tX (b, a) -φ k ( ω ψ a )| ≤ ˜ . (39) (c
| φ k ( ω ψ a )-b a | < ∆ k } and well-separated condition φ k+1 (ω) -φ k (ω) ≥ d k,k+1 ω with d k,k+1 > ω ψ (∆ k + ∆ k+1
). The set Z k we define in Lemma 3 is slightly different from the set Z k for TSWT [START_REF] Li | Theoretical analysis of time-reassigned synchrosqueezing wavelet transform[END_REF], we assume d k,k+1 and ∆ k vary with components, but d and ∆ in [START_REF] Li | Theoretical analysis of time-reassigned synchrosqueezing wavelet transform[END_REF] are constant.

The focus here is on proving that: [START_REF] Cohen | Time-frequency analysis: theory and applications[END_REF] the relationship between T k satisfying weakly-separated condition and Z k satisfying well-separated condition; (2) if the change-rate (in frequency) of the A k (ω) and φ k (ω) is small, compared with the change-rate of the GD φ k (ω) themselves, then the proposed method will identify and characterize the information of the signal X(ω) and their GDs. It can be concluded that, for frequency-varying signals with weakly-separated chirp-like components, the NTEWT can effectively represent them with satisfied TF energy concentration and reasonable reconstruction accuracy. Besides, in the future, we will try to follow the ideas and analysis methods of Daubechies [START_REF] Daubechies | ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform[END_REF] and Sejdić [START_REF] Sejdić | Quantitative Performance Analysis of Scalogram as Instantaneous Frequency Estimator[END_REF] about the theoretical analysis of noise, to further carry out a detailed and in-depth theoretical analysis of noise.

Numerical Analysis of the Behavior of NTEWT and Comparisons

In this section, we will compare the difference between the proposed NTEWT method and advanced TFA methods in addressing simulated and real signals, including SWT [START_REF] Daubechies | Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool[END_REF], MSWT [START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF], RS [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF], TSST [START_REF] He | Time-reassigned synchrosqueezing transform: the algorithm and its applications in mechanical signal processing[END_REF], and TSST2 [START_REF] He | Gaussian-modulated linear group delay model-application to second-order time-reassigned synchrosqueezing transform[END_REF], and we use three quantitative indicators, including the Rényi entropy [START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF], mean relative error (MRE) [START_REF] Li | Self-matched extracting wavelet transform and signal reconstruction[END_REF] and reconstruction quality factor (RQF) [START_REF] Fourer | Second-order time-reassigned synchrosqueezing transform: application to draupner wave analysis[END_REF], to evaluate the TF energy concentration, GD estimation and signal reconstruction performance, respectively.

Simulation study and performance analysis

The TF concentration is one of the outstanding features used for evaluating the TF performance of TFA methods, which can be quantified by Rényi entropy. The Rényi entropy of α order for TFR is defined as [START_REF] Wang | Matching synchrosqueezing transform: a useful tool for characterizing signals with fastvarying instantaneous frequency and application to machine fault diagnosis[END_REF]]

R α = 1 1 -α log 2 ( T F R(b, ω) T F R(b, ω)dbdω ) α dbdω. ( 40
)
where a smaller Rényi entropy value indicates that the TFA method can generate a more energyconcentrated TFR.

Herein, a simulated signal consisting of two components can be constructed in the frequency domain as:

X(ω) = X 1 (ω) + X 2 (ω) = A 1 (ω)e -iφ1(ω) + A 2 (ω)e -iφ2(ω) , (41)
where the higher the frequency, the closer the distance between the two components of X(ω): Rényi entropy, and suffers the least interference at high frequency, which means the NTEWT behaves the best TF concentration performance for X(ω). Meanwhile, Table 1 lists the numerical execution time of these TFA methods in addressing frequency-varying signals X(ω). It can be seen that under the framework of STFT, TSST and TSST2 can finish the analysis within one second. However, these TFA methods under the framework of WT need more time than TSST and TSST2, while NTEWT needs much less time than SWT, MSWT and RS.

A 1 (ω) = A 2 (ω) = e 0.001ω , φ 1 (ω) = -0.001ω 3 + 0.16ω 2 + 1.2ω,
In order to study the performance of TF concentration of those TFA methods in the presence of noise, we add Gaussian white noise with signal-to-noise ratio (SNR) of 0 dB to the test signal (see Fig. 3 and Fig. 4). Compared with other TF results, NTEWT achieves the highest resolution TFR and obtains the follows that among these TFA methods, NTEWT obtains the lowest Rényi entropy and behaves the best TF concentration performance whatever the noise level. Thus, NTEWT provides the most concentrated TFRs for analyzing noisy signals or noise-free signals.

Furthermore, we use the second quantified indicator, i.e., MRE, to evaluate the performance of N-TEWT on detecting the GD feature of the signal, MRE can be employed to measure the errors of estimated GD, which is defined by [START_REF] Li | Self-matched extracting wavelet transform and signal reconstruction[END_REF],

MRE = 1 N G GD -GD GD 1 , (43)
where • 1 denotes l 1 -norm, N G is the discrete length of the GD, GD and GD respectively represent the true GD and the estimated GD. Generally, a lower MRE value means a better GD estimation performance.

In Fig. 6, red lines denote true GD trajectory, blue lines denote estimated GD trajectories, and MRE values are written on both sides. As shown in Fig. 6, for two TFA post-processing methods along the frequency direction, the GD trajectory estimated by SWT deviates from the true GD trajectory most seriously, and as the second-order extension of SWT, MSWT reduces the deviation of both. However, for three TFA post-processing methods along the time direction, the GD trajectories estimated by TSST, TSST2 and NTEWT are closer to the true GD trajectory. Moreover, for RS along the frequency direction and the time direction, the GD trajectory estimated by RS is also very close to the true GD trajectory. It can be also seen that the GD trajectories estimated by the SWT, MSWT, TSST and TSST2 are interfered at high frequency, while the GD trajectories estimated by RS and NTEWT are almost unaffected at high frequency, and highly consistent with the true GD trajectory. And the MREs of GD estimated by RS and NTEWT are smaller than those of the other four TFA methods. Thus, RS and NTEWT shows better GD estimation performance.

Finally, the third quantified indicator RQF is used to evaluate the reconstruction performance, which is defined as [START_REF] Fourer | Second-order time-reassigned synchrosqueezing transform: application to draupner wave analysis[END_REF] RQF = 10 log 10

x(t) 2 x(t) -x r (t) 2 , (44)
where x(t) and x r (t) denote the original signal and the reconstructed signal, respectively. A higher RQF value indicates a better reconstruction performance. Figs. 7,8 show reconstructed results of five TFA methods and their respective RQF values, and we find that the RQF values of TSST, TSST2 and NTEWT are lower than those of SWT and MSWT, and NTEWT achieves the smallest reconstruction error and the highest RQF value, it means that NTEWT behaves better reconstruction performance than other four reconstruction methods. Therefore, it can be concluded that, the TFA methods for post-processing operations along the time direction has more obvious advantages, in analyzing the strongly-varying signals in the frequency domain, than the TFA methods for post-processing operations along the frequency direction, especially the NTEWT behaves the best TF concentration, GD estimation and reconstruction performance for strongly frequency-varying signal X(ω) through three quantified indicators.

Application

In this subsection, we employe two real-world signals to validate the effectiveness of the proposed method. The first case is a popular bat signal recorded by Rice University [START_REF] Yu | Synchroextracting transform[END_REF], its waveform and spectrum are displayed in Fig. 9 Next, we consider a bearing vibration signal with the outer race fault, which was recorded by the Case Western Reserve University [START_REF] He | Gaussian-modulated linear group delay model-application to second-order time-reassigned synchrosqueezing transform[END_REF], and its waveform and spectrum are displayed in Fig. 12. As shown in Fig. 12 (b), the frequency band of this vibration signal is mainly concentrated in the frequency range of 2.5-4 kHz, thus we compare the TF features of this frequency range in Figs. [START_REF] Kodera | A new method for the numerical analysis of nonstationary signals[END_REF][START_REF] Kodera | Analysis of time-varying signals with small BT values[END_REF]. Focus on the TF results and Rényi entropies, the NTEWT can provide a more concentrated TFR than other five TFA methods. This thus demonstrates the interest of the proposed new technique in real applications.

Conclusion

The main contribution of this paper is the introduction of Newton GD estimator and NTEWT, and its theoretical analysis. In this paper, we firstly convert the GD estimation problem into solving the fixed point problem, and combine with Newton algorithm to propose the Newton GD estimator for a direction, especially the NTEWT has a more concentrated TFR, better GD estimation performance, and good reconstruction performance. In the future, we will try to develop a rigorous theoretical analysis of noise, and further explore the potential of the proposed method in practical applications. Proof. (1) On the one hand, based on Eq. ( 38), Π k (a) and Π k+1 (a) can be written as

Π k (a) = min{φ k ( ω ψ a ) + a∆ k , φ k+1 ( ω ψ a ) -a∆ k+1 }, Π k+1 (a) = max{φ k ( ω ψ a ) + a∆ k , φ k+1 ( ω ψ a ) -a∆ k+1 }, (46)
it is obvious that Π k (a) ≤ Π k+1 (a), thus

T k ∩ T k+1 = ∅. Combined with Π k (a) ≤ φ k+1 ( ω ψ a ) -a∆ k+1 , Z k+1 = {(b, a); φ k+1 ( ω ψ a ) -a∆ k+1 < b < φ k+1 ( ω ψ a ) + a∆ k+1 }, (47)
it means that the lower boundary of the zones Z k+1 is above the upper boundary of the zones T k , i.e.

T k ∩ Z k+1 = ∅.
On the other hand, for any (b, a)

∈ T k , i.e. Π k (a) < b < Π k (a), combined with b > Π k (a) ≥ φ k ( ω ψ a ) -a∆ k , b < Π k (a) ≤ φ k ( ω ψ a ) + a∆ k , (48)
this gives (b, a) ∈ Z k , and this gives

T k ⊂ Z k . (2) Since X(ω) = K k=1 X k (ω) ∈ C ,d , d k,k+1 satisfies the separation condition φ k+1 ( ω ψ a ) -φ k ( ω ψ a ) ≥ ad k,k+1 ω ψ . When d k,k+1 > max{ω ψ ∆ k , ω ψ ∆ k+1 }, we have φ k+1 ( ω ψ a ) -φ k ( ω ψ a ) ≥ ad k,k+1 ω ψ > max{a∆ k , a∆ k+1 }, (49)
φ k+1 ( ω ψ a ) > φ k ( ω ψ a ) + a∆ k , φ k+1 ( ω ψ a ) -a∆ k+1 > φ k ( ω ψ a ). (50)
Considering

φ k ( ω ψ a ) > max{φ k-1 ( ω ψ a ) + a∆ k-1 , φ k ( ω ψ a ) -a∆ k } = Π k (a), (51)
and

φ k ( ω ψ a ) < min{φ k ( ω ψ a ) + a∆ k , φ k+1 ( ω ψ a ) -a∆ k+1 } = Π k (a), (52)
ω ψ a ) -φ k ( ω ψ a ) ≥ ad k,k+1 ω ψ > a(∆ k + ∆ k+1 ), which means that φ k+1 ( ω ψ a ) -a∆ k+1 > φ k ( ω ψ a ) + a∆ k . This leads to Π k (a) = max{φ k-1 ( ω ψ a ) + a∆ k-1 , φ k ( ω ψ a ) -a∆ k } = φ k ( ω ψ a ) -a∆ k , Π k (a) = min{φ k ( ω ψ a ) + a∆ k , φ k+1 ( ω ψ a ) -a∆ k+1 } = φ k ( ω ψ a ) + a∆ k . (53)
Thus, when d k,k+1 > ω ψ (∆ k + ∆ k+1 ), we have

T k = Z k .
The above lemma indicates that:

(1) Each zone T k does not intersect with each other, T k is a subset of Z k . Moreover, T k and Z ( = k) don't intersect, which means that the zone T k for component X k (ω) is separated from the zone Z for other component X (ω), and it has been proved that in [START_REF] Li | Theoretical analysis of time-reassigned synchrosqueezing wavelet transform[END_REF], Z provides an estimate of the main energy distribution region of W ψ X (b, a). In what follows, denote

M k = sup ω∈R A k (ω) < ∞, M k = max(sup ω∈R |A k (ω)|, sup ω∈R φ k (ω)) < ∞, M k = sup ω∈R |φ k (ω)| < ∞. If O(A k ) and O(φ k ) are neglected, the kth component X k (ω)

can be approximated by

Taylor expansion as

X k (ω) ≈ Xk (ω) = A k ( ω ψ a )e -i[φ k ( ω ψ a )+φ k ( ω ψ a )(ω- ω ψ a )+ 1 2 φ k ( ω ψ a )(ω- ω ψ a ) 2 ] . (54)
To simplify some notations, denote t k (b, a) =

φ k ( ω ψ a )-b a . If (b, a) ∈ T k , then |t k (b, a)| < ∆ k . Set Ψ k (t k (b, a)) = ψ * φ (t k (b,a)) ψ * φ (t k (b,a))
, and for (b, a) ∈ T k , assume

|ψ φ (t k (b, a))| ≤ Q k , |Ψ k (t k (b, a))| ≤ Q k , |Ψ k (t k (b, a))| ≤ Q k . (55)
Estimate 1. For each k ∈ {1, . . . , K}, we have

|A k (ω) -A k ( ω ψ a )| ≤ M k |ω - ω ψ a | |φ k (ω) -φ k ( ω ψ a ) -φ k ( ω ψ a )(ω - ω ψ a ) - 1 2 φ k ( ω ψ a )(ω - ω ψ a ) 2 | ≤ 1 6 M k |ω - ω ψ a | 3 (56)
Proof. When ω ≥ ω ψ a , noticing Definition 10, we have (the case ω < ω ψ a can be done in a similar way):

|A k (ω) -A k ( ω ψ a )| = | ω ω ψ a A k (x)dx| ≤ ω ω ψ a |A k (x)|dx ≤ ω ω ψ a |φ k (x)|dx ≤ M k |ω - ω ψ a |, (57)
On the other hand, it is easy to show that

|φ k (ω) -φ k ( ω ψ a ) -φ k ( ω ψ a )(ω - ω ψ a ) - 1 2 φ k ( ω ψ a )(ω - ω ψ a ) 2 | = | ω ω ψ a x ω ψ a [φ k (y) -φ k ( ω ψ a )]dydx| = | ω ω ψ a x ω ψ a y ω ψ a φ k (τ )dτ dydx| ≤ M k ω ω ψ a x ω ψ a |y - ω ψ a |dydx = 1 6 M k |ω - ω ψ a | 3 . (58)
Next, we will prove Theorem 2 (b) by approximating Xk (ω) to X k (ω), which first needs to calculate several WTs of Xk (ω) and their upper bounds.

where

Υ 0,0 (a) = M k Q k 2πa , Υ 0,1 (a) = M k Q k (∆ k + M k a 2 Q k ) 2πa , Υ 1,0 (a) = M k Q k (ω ψ + Q k ) 2πa 2 , Υ 1,1 (a) = M k Q k 2πa 2 (1 + ∆ k Q k + ∆ k ω ψ + M k a 2 (ω ψ Q k + Q k + Q 2 k )). (60)
(2) Furthermore, for (b, a) ∈ T k , for any 1

≤ = k ≤ K, denote β = sup (b,a)∈R×R + {|t (b, a)|}, we have |W ψ X (b, a)| ≤ J 0,0 (a), |W ψ X (b, a)| ≤ J 0,1 (a) |∂ b W ψ X (b, a)| ≤ J 1,0 (a), |∂ b W ψ X (b, a)| ≤ J 1,1 (a) (61)
where

J 0,0 (a) = M G 0 2πa , J 0,1 (a) = M 2πa 3 (a 2 G 0 β + G 1 M ), J 1,0 (a) = M 2πa 2 (G 0 ω ψ + G 1 ), J 1,1 (a) = M 2πa 2 (G 0 + ω ψ (G 0 β + G 1 M a 2 ) + β G 1 + M G 2 a 2 ) (62)
Proof.

(1) Based on Eqs. ( 5) and ( 54), the WT W ψ Xk (b, a) of Xk (ω) with ψ(ω) can be expressed as

W ψ Xk (b, a) = 1 2π R A k ( ω ψ a )e -i[φ k ( ω ψ a )+φ k ( ω ψ a )(ω- ω ψ a )+φ k ( ω ψ a ) (ω- ω ψ a ) 2 2 ] ψ * (aω)e iωb dω = 1 2πa A k ( ω ψ a )e -iφ k ( ω ψ a ) e i ω ψ a b R e -iφ k ( ω ψ a ) ω 2 2a 2 ψ * (ω + ω ψ )e -i(φ k ( ω ψ a )-b) ω a dω = 1 2πa X k ( ω ψ a )e i ω ψ a b R ψ * φ (ω)e -i(φ k ( ω ψ a )-b) ω a dω = 1 2πa X k ( ω ψ a )e i ω ψ a b ψ * φ (t k (b, a)). (63)
= it k (b, a) - φ k ( ω ψ a ) a 2 Ψ(t k (b, a)), ∂ b W ψ Xk (b, a) = Λ 2,k (a, b)W ψ Xk (b, a), Λ 2,k (a, b) = i ω ψ a - 1 a Ψ(t k (b, a)), ∂ b W ψ Xk (b, a) = Λ 3,k (a, b)W ψ Xk (b, a), Λ 3,k (a, b) = -( i + it k (b, a)Ψ(t k (b, a)) + ω ψ t k (b, a) a + φ k ( ω ψ a ) a 3 (iω ψ Ψ(t k (b, a)) -Ψ (t k (b, a)) -Ψ 2 (t k (b, a)))). (64)
From the above results, it is easy to obtain the bounds of these WTs.

(

φ (t), if |t (b, a)| ≥ ∆ with d , +1 > max{ω ψ ∆ , ω ψ ∆ +1 }, then |ψ φ (t (b, a))| ≤ G 0 , |ψ φ (t (b, a))| ≤ G 1 , |ψ φ (t (b, a))| ≤ G 2 .
Thus, by Eqs. ( 63) and (64), we obtain

|W ψ X (b, a)| = | 1 2πa X ( ω ψ a )e i ω ψ a b ψ * φ (t (b, a))| ≤ M G 0 2πa = J 0,0 (a), ( 65
) |W ψ X (b, a)| ≤ |it (b, a)W ψ X (b, a)| + | 1 2πa 3 X ( ω ψ a )e i ω ψ a b φ ( ω ψ a )ψ * φ (t (b, a))| ≤ β M G 0 2πa + M G 1 2πa 3 M = M 2πa 3 (a 2 G 0 β + G 1 M ) = J 0,1 (a), (66) |∂ b W ψ X (b, a)| ≤ |i ω ψ a W ψ X (b, a)| + | 1 2πa 2 X ( ω ψ a )e i ω ψ a b ψ * φ (t (b, a))| ≤ M 2πa 2 (G 0 ω ψ + G 1 ) = J 1,0 (a), (67) and
M G 0 2πa 2 + M β G 1 2πa 2 + M M G 2 2πa 4 = M 2πa 2 (G 0 + ω ψ (G 0 β + G 1 M a 2 ) + β G 1 + M G 2 a 2 ) = J 1,1 (a).
Next, we will consider the error analysis of the Newton GD estimator and the true GD. (91)

Proof of Theorem 2 (b)

The proof of Theorem 2 is finished.

  On the one hand, based on different features of signals, signals are divided into time-varying signals and frequency-varying signals. The TF energy distribution of time-varying signals spread around the IF trajectory, which can help to better identify the variation law of time-varying signals along the time direction, while the TF energy distribution of frequency-varying signals spread around the GD trajectory, which can effectively reflect the variation law of frequency-varying signals along the frequency direction.

  (a) tX (b, a) = b if and only if ∂ b ln M ψ X (b, a) = 0. (b) If tX (b, a) is a contraction operator for a given scale a, then tX (b, a) = b if and only if the point (b, a) is a ridge point. Proof. (a) For ψ(ω), since tψ(t) F T

( b )

 b Assume b * is a fixed point of tX (b, a), i.e. tX (b * , a) = b * . Then, for proving that (b * , a) is a ridge point, we only need to show ∂ 2 ∂b 2 ln M ψ X (b * , a) < 0, by Eq. (12) which is equivalent to prove ∂ b tX (b * , a) < 1. Consider a sequence {b n } +∞ n=1 with lim n→+∞ b n = b * , and denote F (b) = tX (b, a). Since tX (b, a) is a contraction operator, by Definition 6, there exists a constant k such that 0 < k < 1 and

  a) is a ridge point, then by (a), we have tX (b, a) = b. Hence, Theorem 1 (b) holds.

  (a) tX (b, a) is a contraction operator. (b) tX (b, a) is a biased estimate of the signal GD. (c) tX (b, a) = b if and only if the point (b, a) is at the GD trajectory of signal X(ω).

(

  d) tX (b, a) = b if and only if the point (b, a) is a ridge point of WT W ψ X (b, a).

( 1 )

 1 tX (b, a) is a biased estimate of the GD of X(ω); (2) the fixed point of time reassignment operator tX (b, a) is an accurate estimate of the GD; (3) the GD trajectory of signal X(ω) and the ridge trajectory of WT W ψ X (b, a) are the same.

  (a) tX (b, a) = b if and only if tX (b, a) = b. (b) tX (b, a) is an accurate GD estimation of the signal X(ω).

  (c) tX (b, a) = b if and only if the point (b, a) is at the GD trajectory of signal X(ω).

  (d) tX (b, a) = b if and only if the point (b, a) is a ridge point of WT W ψ X (b, a). Proof. (a) By Eq. (

  NTEWT only extracts the WT coefficients at the fixed point trajectory b = tX (b, a), while removing the TF coefficients beyond the fixed point trajectory b = tX (b, a). In particular, for the frequency-varying linear chirp signal X(ω), it can be concluded from the above discussion that: (1) by Lemma 2 (c), the NTEWT exactly extracts the WT coefficients at the GD trajectory of signal; (2) as shown in Lemma 2 (d), the fixed points of the NGD estimator tX (b, a) are the time ridge points, which means that NTEWT can remain the maximum energy of the WT modulus, then NTEWT can obtain a more energy concentrated TFR;

  characterize a function class C ,d . To better understand the function class C ,d for NTEWT, let's first review a function class H ,d for TSWT. Definition 9. [38] The class H ,d is said to be a set of all superposition of well-separated frequencydomain harmonic-like functions, with accuracy > 0 and separation

  and the energy distribution of W ψ X (b, a) is concentrated in the zones 1≤k≤K Z k . Thus, zone Z k provides an estimate of the main energy distribution region of W ψ X k (b, a). Compared with the class H ,d for TSWT, the class C ,d for NTEWT is not restricted to the class H ,d composed of several harmonic-like components, which is mathematically defined as: Definition 10. The class C ,d is said to be a set of all superposition of weakly-separated frequency-varying chirp-like components (WFCC) up to accuracy > 0 and with separation d k,k+1 > 0, if each element

  )Intuitively, the function class C ,d is composed of several oscillatory components X k (ω), where the change rate (in frequency) of A k (ω) and φ k (ω) is much smaller than the change rate of φ k (ω) itself, which means that each component X k (ω) can be viewed as approximately a chirp signal locally in the frequency domain, with slowly varying amplitude A k (ω) and φ k (ω). It is worth noting that, different from the constant separation constant d in the well-separated condition[START_REF] He | Time-reassigned synchrosqueezing transform: the algorithm and its applications in mechanical signal processing[END_REF], each separation constant d k,k+1 in the weakly-separated condition (37) here varies with components.

1 4.

 1 d , and set ˜ := Pick a wavelet ψ(ω) such that inverse FT ψ φ (t) of its variant ψφ (ω) = e iφ k ( ω ψ a ) ω 2 2a 2 ψ(ω + ω ψ ) satisfies the following conditions: if |t| > ∆ k with

  (a) For k ∈ {1, 2, ..., K -1}, and for each pair (b, a), we have T k = ∅ with (φ k ( ω ψ a ), a) ∈ T k , and T k ∩ T k+1 = ∅. (b) For each k ∈ {1, . . . , K -1}, and for each pair (b, a) ∈ T k for which holds |W ψ X (b, a)| > ˜ , then we have

Theorem 2

 2 ) Moreover, for each k ∈ {1, . . . , K}, and for each pair (b, a) ∈ T k , there exists a constant C > 0 such that, | 2πaN T e(φ k ( ω ψ a ), a)e -i ω ψ a φ k ( basically tells us that: (a) If d k,k+1 > max{ω ψ ∆ k , ω ψ ∆ k+1 }, then weakly-separated chirp-like components of C ,d may interfere with each other. However, each zone T k does not intersect with each other, the zone T k corresponds to the energy distribution region of W ψ X k (b, a) of the kth component X k (ω) at the point (φ k ( ω ψ a ), a) and its vicinity, and the zone T k for component X k (ω) is separated from the main energy distribution regions of WT of other components. (b) The NGD estimator tX (b, a) is indeed a high-precision GD estimator in the zone T k . (c) The signal can be reconstructed with reasonably high accuracy by NTEWT. The proof of this theorem is available in Appendix A. In the Lemma 3 of Appendix A, to further compare T k and Z k , we modify the set Z k for TSWT [38], i.e. Z k = {(b, a); | φ k ( ω ψ a )-b a | < ∆} and well-separated condition φ k+1 (ω) -φ k (ω) ≥ d ω with d > 2ω ψ ∆. In Lemma 3, we define a set Z k := {(b, a);

φ 2 ( 42 )Figs. 1 ,Figure 1 :Figure 2 :

 242112 Figs. 1, 2 show the TF results obtained by different TFA methods in addressing frequency-varying signals X(ω), where third-order Rényi entropy values are displayed in the bottom. It can be observed that some TFA post-processing methods along the frequency direction, such as SWT and MSWT, provide blurred TFRs for frequency-varying signals X(ω). However, some TFA post-processing methods along the time direction, such as TSST, TSST2 and NTEWT, improve the readability of TFRs of frequencyvarying signals X(ω). The RS along the time and frequency directions also obtain a high resolution TFR. Compared with the other five methods, NTEWT generates the most concentrated TFR, obtains the least

Figure 3 :

 3 Figure 3: TF results for X(ω) with SNR=0 dB Gaussian white noise: (a) SWT, (b) MSWT, (c) RS, (d) zoom of the SWT, (e) zoom of the MSWT, (f) zoom of the RS.

Figure 4 :

 4 Figure 4: TF results for X(ω) with SNR=0 dB Gaussian white noise: (a) TSST, (b) TSST2, (c) NTEWT, (d) zoom of the TSST, (e) zoom of the TSST2, (f) zoom of the NTEWT.

Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 5678 Figure 5: Under different noise levels (SNRs of 0-30 (dB)), Rényi entropy curves of TFRs generated by TFA methods.

  (a), (b). We first use NTEWT to reconstruct the signal, and its reconstruction results are shown in Fig. 9 (c), (d). It can be seen that NTEWT can reconstruct signals well with high RQF value and with small reconstruction error, thus NTEWT behaves the good reconstruction ability. Furthermore, in order to understand the non-linear behaviors of the bat signal precisely, we use SWT, MSWT, RS, TSST, TSST2 and NTEWT to analyze the bat signal, and corresponding TF results are shown in Figs. 10, 11. Both the quantitative measure by the Rényi entropy and the visual interpretation show that the NTEWT can accurately characterize the bat signal and essentially improve the TF energy concentration of the bat signal.

Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 9101112 Figure 9: (a) The waveform of the bat signal. (b) The spectrum of the bat signal. (c) The original signal and the NTEWT reconstructed signal. (d) The reconstructed signal and the NTEWT reconstructed error.

Figure 13 :

 13 Figure 13: TF results for vibration signal: (a) SWT, (b) zoom of the SWT, (c) MSWT, (d) zoom of the MSWT, (e) RS, (f) zoom of the RS.

Figure 14 :

 14 Figure 14: TF results for vibration signal: (a) TSST, (b) zoom of the TSST, (c) TSST2, (d) zoom of the TSST2, (e) NTEWT, (f) zoom of the NTEWT.

  leads to (φ k ( ω ψ a ), a) ∈ T k , which means that T k = ∅.(3) According to separation condition and d k,k+1 > ω ψ (∆ k + ∆ k+1 ), we have φ k+1 (

( 2 )

 2 It is noteworthy that when d k,k+1 > max{ω ψ ∆ k , ω ψ ∆ k+1 }, weakly-separated chirp-like components of C ,d may interfere with each other, and the zone T k corresponds to the energy distribution region of W ψ X k (b, a) of the kth component X k (ω) at the point (φ k ( ω ψ a ), a) and its vicinity. (3) If d k,k+1 > ω ψ (∆ k +∆ k+1 ), then we have T k = Z k , it follows that weakly-separated conditions are restricted to well-separated cases, and the class C ,d contains some signals with well-separated components.

Estimate 2 . ( 1 )

 21 For each k ∈ {1, . . . , K} and (b, a) ∈ T k , we have |W ψ Xk (b, a)| ≤ Υ 0,0 (a), |W ψ Xk (b, a)| ≤ Υ 0,1 (a), |∂ b W ψ Xk (b, a)| ≤ Υ 1,0 (a), |∂ b W ψ Xk (b, a)| ≤ Υ 1,1 (a),

Thus

  a)| ≤ 1 2πa M k Q k = Υ 0,0 (a).Furthermore, W ψ Xk (b, a), ∂ b W ψ Xk (b, a) and ∂ b W a) = Λ 1,k (a, b)W ψ Xk (b, a), Λ 1,k (a, b)

)

  According to Lemma 3, for (b, a) ∈ T k and k = , we have T k ∩ Z = ∅, and (b, a) / ∈ Z , i.e., |t (b, a)| ≥ ∆ . According to the condition on ψ

,

  ω ψ a b ψ * φ (t (b, a))| ≤ ω ψ M 2πa 4 (a 2 G 0 β + G 1 M ) +

Estimate 3 .X) 2 ]Estimate 4 . ( 1 )

 3241 For k ∈ {1, . . . , K}, and (b, a) ∈ R × R + , we have|W ψ X k (b, a) -W ψ Xk (b, a)| ≤ M k (k (b, a) -W ψ Xk (b, a)| ≤ M k (I 1 0,1 (a) + M k 6 I 3 0,1 (a)) |∂ b W ψ X k (b, a) -∂ b W ψ Xk (b, a)| ≤ M k (I 1 1,0 (a) + M k 6 I 3 1,0 (a)) |∂ b W ψ X k (b, a) -∂ b W ψ Xk (b, a)| ≤ M k (I 1 (a) = 1 2πa p R |ω -ω ψ a | r |(aω) p ψ(q) (aω)|dω, p, q = 0, 1, r = 1, 3. Proof. According to the WT definition, for (b, a) ∈ R × R + , we have |W ψ X k (b, a) -W ψ Xk (b, a)| ≤ | 1 2π R (A k (ω) -A k ( ω ψa))e -iφ k (ω) ψ * (aω)e iωb dω|+ | 1 2π R A k ( ω ψ a )(e -iφ k (ω) -e -i[φ k ( ω ψ a )+φ k ( ) ψ * (aω)e iωb dω| ≤ 1 2π R |A k (ω) -A k ( ω ψ a )|| ψ(aω)|dω + 1 2π R |A k ( ω ψ a )||φ k (ω) -φ k ( ω ψ a ) -φ k ( the differentialmean value theorem e ix -e i0 = ie iξ x for the second inequality. The proofs of the other estimates are analogous. For k ∈ {1, . . . , K}, and (b, a) ∈ R × R + , we have |W ψ X (b, a) -(a)), |∂ b W ψ X (b, a) -K k=1 ∂ b W ψ Xk (b, a)| ≤ K k=1 M k (I 1 1,0 (a) +

Estimate 6 . 2 Ψ( 2 ) 1 2 2 ) 1 2 H ( 1 4 2 M

 622121212 For k ∈ {1, . . . , K}, and (b, a) ∈ T k such that |W ψ X (b, a)| > ˜ , we have (1) Denote P k = inf a∈R + {|1+i φ k ( (t k (b, a)|}, if there exists some constant S k such that |Ψ(t k (b, a))t k (b, a)Ψ (t k (b, a))| ≤ S k , then | t Xk (b, a) -φ k ( Furthermore, we have | tX (b, a) -φ k ( Γ1(a)) H , with H = inf a∈R + {|1 -∂ b tX (b, a)|}.Proof. (1) By Eq. (32), we have| t Xk (b, a) -φ k ( ω ψ a )| = |b -b -φ k ( t k (b, a)) -t k (b, a)Ψ (t k (b, a)))According to the definition of tX (b, a), we have then| tX (b, a) -φ k ( ω ψ a )| = |b -b -tX (b, a) 1 -∂ b tX (b, a) -φ k ( ω ψ a )| = | b -tX (b, a) -(b -φ k ( ω ψ a ))(1 -∂ b tX (b, a)) 1 -∂ b tX (b, a) | ≤ | b -tX (b, a) -(b -t Xk (b, a)) 1 -∂ b tX (b, a) | + | (b -φ k ( ω ψ a ))(1 -∂ b tX (b, a)) -(b -φ k ( ω ψ a ))(1 -∂ b t Xk (b, a)) 1 -∂ b tX (b, a) | + | b -t Xk (b, a) -(b -φ k ( ω ψ a ))(1 -∂ b t Xk (b, a)) 1 -∂ b tX (b, a) | ≤ | tX (b, a) -t Xk (b, a) 1 -∂ b tX (b, a) | + | (b -φ k ( ω ψ a ))(∂ b tX (b, a) -∂ b t Xk (b, a)) 1 -∂ b tX (b, a) | + | ( t Xk (b, a) -φ k ( ω ψ a ))(1 -∂ b t Xk (b, a)) 1 -∂ b tX (b,a) | ≤ Γ 0 (a) + a∆ k Γ 1 (a) + 1 b tX (b,a) | = |1 + ∂ b tX (b,a)-∂ b t Xk (b,a) 1-∂ b tX (b,a)

  Assume

1 4 Γ 2

 12 (a) < 1, then we have| tX (b, a) -φ k ( ω ψ a )| ≤ ˜ . (88)Proof of Theorem 2 (c). For (b, a) ∈ T k , we evaluate the WT in Eq. (63) along the trajectory b = φ k ( |Ω 0,0 (a) = C˜ .

Table 1 :

 1 REQUIRED COMPUTATION TIME.

	TFA	SWT MSWT	RS	TSST TSST2 NTEWT
	Time(s) 4.749	7.928	5.255 0.134	0.211	2.663

Acknowledgements

This work was supported by National Natural Science Foundation of China under Grant U20B2075, supported by the ANR ASCETE project with Grant number ANR-19-CE48-0001-01, supported by the Fundamental Research Funds for the Central Universities under Grant G2021KY05103, and supported by the Natural Science Foundation of Shaanxi Province under Grant 2023-JC-QN-0739. The authors thank the editor and the anonymous reviewers for their valuable comments and suggestions that helped improve the quality of this paper.

Data availability statement

The authors do not have permission to share data.

Appendix A

For convenience, we divide the proof of Theorem 2 into several simple estimates and demonstrate one by one. The following lemma clearly explains the relationship between T k satisfying weakly-separated condition and Z k satisfying well-separated condition. Lemma 3. For k ∈ {1, 2, ..., K}, define a set Z k

then the following hold:

(1)

(2) Furthermore, for (b, a) ∈ T k , we have

where Ω p,q (a) = K k=1

(2) For (b, a) ∈ T k , according to Lemma 3, we have (b, a) /

Thus, by Estimate 2 and Estimate 3, it follows that

The proofs of the other estimates are analogous.

Proof.

(1) Based on Eqs. ( 19), ( 64) and (72), we obtain

(2) By Eq. ( 19), we have

this gives

then, we can get