A Fast and Robust Spectrogram Reassignment Method

: The improvement of the readability of time-frequency transforms is an important topic in the ﬁeld of fast-oscillating signal processing. The reassignment method is often used due to its adaptivity to different transforms and nice formal properties. However, it strongly depends on the selection of the analysis window and it requires the computation of the same transform using three different but well-deﬁned windows. The aim of this work is to provide a simple method for spectrogram reassignment, named FIRST (Fast Iterative and Robust Reassignment Thinning), with comparable or better precision than classical reassignment method, a reduced computational effort, and a near independence of the adopted analysis window. To this aim, the time-frequency evolution of a multicomponent signal is formally provided and, based on this law, only a subset of time-frequency points is used to improve spectrogram readability. Those points are the ones less inﬂuenced by interfering components. Preliminary results show that the proposed method can efﬁciently reassign spectrograms more accurately than the classical method in the case of interfering signal components, with a signiﬁcant gain in terms of required computational effort.


Introduction
A sparse representation of a function strongly depends on its properties.This task is made more difficult in the case of non-stationary time-frequency functions [1][2][3][4], as a single linear transform is not able to provide a really sparse representation of the function.That is why either nonlinear transforms or nonlinear operations in the transformed domain are necessary.
In this paper, time-varying oscillatory functions have been considered.More precisely, we focus on amplitude and frequency-modulated multicomponent signals which are defined as [1,4] where f k is the k-th mode, a k and φ k respectively are its amplitude and phase functions of the time variable t while M is the number of modes.Equation (1) represents chirp-like signals, as for example radar, audio, speech and seismic signals, gravitational waves, and so on.In several applications it is required to separate each single mode of such kind of signals to retrieve its phase.Specifically, the main goal is to estimate the instantaneous frequency (IF) of each mode, i.e., the positive value of the derivative with respect to the time variable of the phase φ k (t) [5][6][7].As a result, it is desirable to have a transform which can well separate the modes of the signal to allow a straightforward estimation of the corresponding instantaneous frequencies.The literature concerning linear and quadratic time-frequency analysis [2,4,[8][9][10][11][12][13] has proved that localized time-frequency analysis is able to provide interesting and quite efficient solutions to the problem.Some examples are the Short-Term Fourier Transform (STFT), Gabor Transform, Wavelet Transform, etc.All these transforms are based on the selection of a proper analysis window with time-frequency compaction properties whose limit is subjected to the Heisenberg principle.Unfortunately, although window-dependent, time-frequency transforms are not able to sparsify and separate each single mode with an acceptable precision due to interference patterns in the time-frequency plane.Additional operators have to be applied to the adopted time-frequency transform in order to further compact information; some interesting examples are the reassignment method [14][15][16][17][18] and the more recent synchrosqueezing transform [19][20][21][22][23][24][25] that can be applied to a wide class of transforms, hold for discrete signals and can be extended to filter banks [26,27].Their main idea is to move each point toward the center of mass of the energy distribution.The latter is evaluated in an interval whose amplitude is consistent with the support of the analysis window-see Figure 1.Unfortunately, modes separation is possible whenever the separability condition is satisfied, i.e., whenever the distance between the instantaneous frequencies of each pair of modes is at least equal to the bandwidth of the analysis window.If the separability condition is not satisfied, the reassignment method can fail since centroids cannot be correctly estimated (see Figure 2); for similar reasons, the synchrosqueezing does not separate the modes.Even though the use of a different window can sometimes help in better separating signal components in the time-frequency plane, it is not guaranteed the existence of an optimal window for separating all signal components; in addition, for those signals whose instantaneous frequency curves intersect, it does not exist a window able to separate them at each point.An accurate separation of signal modes at each point of the time-frequency plane allows for signal identification even from a limited region of the plane, resulting useful in the case of signals with close IFs in the observation time interval or finite-length signals.This paper aims at giving a step forward in this sense by providing a different pointwise reassignment strategy which is more robust to the lack of separability of signal modes, without using a priori estimation of corresponding IFs or the identification of non-separable regions.To this aim, FM signals with constant amplitude have been considered and FIRST (fast iterative and robust reassignment thinning) has been presented.FIRST mainly relies on the description of the modulus of signal STFT through a proper time-frequency evolution law.The latter allows us to characterize a limited set of points which represent good information to reassign even in the interference region.A local approximation of square root of the spectrogram is then employed to derive a fast algorithm for estimating the amount of reassignment for each selected point.As a result, for each point in the time-frequency plane, reassignment consists of an iterative procedure which gradually converges to the centroid of the energy distribution of the corresponding mode.Preliminary experimental results, achieved on multicomponent signals with constant amplitude, show that FIRST is able to: (i) provide comparable results with the reassignment method in regions where the separability condition is satisfied (non-interfering regions); (ii) improve reassignment result in correspondence to the non-separability region (i.e., the interference region depicted in Figure 2d); (iii) be computationally efficient.
The remainder of the paper is the following.Next section presents a brief formal description of standard reassignment method.Section 3 presents the proposed method and the details concerning the corresponding reassignment procedure.Section 4 contains some numerical results and comparative studies.Finally, the last section draws the conclusions.
where g belongs to the Schwartz class (i.e., the space of smooth functions with fast decaying derivatives of any order) and represents the analysis window- * denotes the complex conjugate.
The spectrogram is then defined as the square modulus of S g f (u, ξ), i.e., It is an element of the Cohen's class (i.e., the family of time-frequency energy distributions covariant by translations in time and frequency), since it is quadratic, time and frequency covariant, and preserves energy.
The spectrogram provides a sparse representation in the time-frequency plane of a multicomponent signal, as in Equation ( 1), if it consists of M distinct curves, each of them corresponding to the IF of each signal single mode.Unfortunately, this requirement cannot always be satisfied, as it happens for the signal shown in Figure 3; in this case the spectrogram spreads, providing a not enough compact signal representation in the time-frequency plane.
Spectrogram spreading is evident if one considers the following equivalent definition [4,14] where is the Wigner-Ville [4] distribution of the function f .The spectrogram of f is then the smoothed version of the Wigner-Ville distribution of f while the smoothing kernel is in turn the Wigner-Ville distribution of the analysis window.The Wigner-Ville distribution is a measure of signal energy density which derives from the correlation of the signal with itself.This is the reason why it guarantees high compaction properties which strictly depend on the time-frequency properties of the function itself.In particular, the following result holds [4,14].Proposition 1.If f (t) = a(t)e iφ(t) , then Hence, the center of mass along the frequency direction of the Wigner-Ville distribution of f is the instantaneous frequency of f .
Unfortunately, for a multicomponent signal the quadratic terms of the Wigner-Ville distribution cause some oscillatory interference patterns that alter the localization of each single IF in the time-frequency plane.As a result, the smoothing introduced by the spectrogram contributes to reduce the oscillatory interference between individual components introduced by the same WV f but, at the same time, it spreads energy density distribution.
The aim of reassignment is to compensate the spreading effects caused by the 2D smoothing introduced by the spectrogram.The main idea is to move the local energy to the centroid of the distribution, whose coordinates are The mean of the reassignment defines the reassigned spectrogram, i.e., Reassignment holds in a region which depends on the support of the analysis window.
Proposition 2. Let ϕ(u, ξ) denote the phase of the STFT S g f (u, ξ) of a function f , computed through the analysis windows g(t), and let S g f (u, ξ) and S tg f (u, ξ) be the STFT of f computed respectively through g (t), tg(t).Let us set where Re( * ) and Im( * ) respectively denote the real and the imaginary part of * ; then 1. Equations ( 4) and ( 5) are equivalent to 2. Equations ( 9) and ( 10) are equivalent to The proof of Proposition 2 can be found in [14,15,28].
It is worth observing that Equations ( 11) and ( 12) provide a straightforward procedure for the estimation of the centroid of the distribution.In fact, they consist of a proper combination of the STFTs of the analyzed function using three different but specific analysis windows, i.e., g, g , tg.On the other hand, the properties of the analysis windows strongly influence reassignment result, especially in the case of multicomponent signals.
One of the main constraints in a reassignment-like method is the separability condition on the individual signal components [15]: time-frequency points are correctly reassigned to their corresponding mode if the following condition is satisfied where ∆ω is the frequency bandwidth of the analysis window.If this property is not met, it is not possible to achieve the desired multicomponents decomposition, independently of the "effectiveness" of the adopted transform.For example, the synchrosqueezing transform is a frequency analysis method that can decompose complex signals into time-varying oscillatory components.It is a form of time-frequency reassignment that is both sparse and invertible, allowing for the recovery of the signal-see [15] for details.Unfortunately, the properties of such transform are guaranteed only if the separability condition is satisfied.The aim of this paper is to define a reassignment method that can characterize signal components with high precision even when the adopted time-frequency transform is not able to resolve signal components.The main idea mainly consists of weakening the separability condition.The goal is to describe the relation between points in the time-frequency plane and to use this relation for selecting those points where the separability condition is nearly verified; as a result, reassignment is only applied to those subsets of points.In other words, the separability condition is translated into the concept of interference between modes in the time-frequency plane; hence, only those points less influenced by interference terms are used for spectrogram reassignment.

Time-Frequency Evolution Law
Let f (t) = a cos(φ(t)) be a frequency-modulated signal with constant amplitude and let S g f (u, ξ) be the Short-Term Fourier Transform of f (STFT) computed through a real and symmetric window g, with support [− 1 2 , 1  2 ] and such that ĝ(0 where s is a fixed scaling factor while ε(u, ξ) is almost negligible.
For the sake of simplicity, in the sequel superscripts will be omitted and, without loss of generality, the scaling parameter s will be set equal to 1; hence, Using the methodologies in [6,29,30], the following Proposition can be proved.
where P u and P ξ are the partial derivatives of P with respect to u and ξ respectively.
Proof of Proposition 3. From Equation ( 14), it follows P(u, ξ) = a 2 4 ĝ2 (ξ − φ (u)) and then By comparing the partial derivatives of P with respect to ξ and u in Equations ( 16) and ( 17) we get the proof.
Equation (15) exactly describes the relation between spectrogram points in the time-frequency plane.In particular, these points belong to specific curves in the time-frequency plane, as proved by the following proposition.
Proof of Proposition 4. Let us consider the parametric curve with τ ∈ R, and P(u(τ), ξ(τ)).The derivative of P(u(τ), ξ(τ)) with respect to τ is where u = du dτ and ξ = dξ dτ .Using Equation ( 15), we have that Finally, the proof derives from the solution of the system of ordinary differential equations in Equation (22).
It is worth observing that the characteristic curves in Equation ( 18) are parallel to the instantaneous frequency curve, which in turn is a particular characteristic curve (c = 0)-see Figure 4.The following property for the spectrogram can be then derived.Proposition 5 implies that points with the same amplitude can be moved in the same way toward their center of mass.
Unfortunately, in case of multicomponent signals, Equation (15) does not hold ∀ (u, ξ) ∈ R 2 ; it turns out that we cannot derive characteristic curves by writing the evolution law in the time-frequency plane for a multicomponent signal.The main motivation is exactly the separability condition in Equation ( 13) that cannot be verified for all (u, ξ) ∈ R 2 .However, Equation ( 18) associates a group of curves C c,φ , depending on the parameter c, to the same component (mode); in particular, the curve with c = 0 is the ridge curve [4].As a consequence, if the ridge curve of a mode does not satisfy the separability condition with respect to another mode, it does not imply that any characteristic curve does-see Figure 5.More precisely, given two modes with IF respectively φ 1 (u) and φ 2 (u), the separability condition can be weakened as follows Hence, the farther the two modes, the more the points satisfying the separability condition in Definition 2. Curves satisfying constraints in Definition 2 can be characterized by the following proposition.Proof of Proposition 6.By setting ω = ξ 1 − φ 2 (u), the proof directly derives from the definition of window bandwidth [4], i.e., ĝ(ω) << 1 for |ω| ≥ ∆ω, when applied to the mode with φ 2 (u) as IF.
Remark 2. It is worth observing that the condition ĝ(ξ By observing that the spectrogram of a two-component signal is where More in general, for fixed u the following inequality holds: Let ξ 1 : = 1 H , and = 1 K , with ĝ symmetric and monotonically decreasing window function e., H < K, and then inequality in Equation ( 24) becomes By accounting for the definition of window bandwidth [4], previous inequality is not trivial for H < √ 2, i.e., for |ξ 1 − φ 1 (u)| < ∆ω 2 .Hence, time-frequency plane points (u, ξ) such that the ratio P(u, ξ) with high probability satisfy the separability condition, as in Definition 2, and then we expect that reassignment can be more correct if it is applied only to those points in the interference region.In our simulations we set ρ = 1 4 .

Fast Iterative and Robust Spectrogram Thinning (FIRST)
The profile of the window function can be used for a direct estimation of the reassignment shift.By considering the regularity of the window function, for a fixed u, we can compute the Taylor expansion around the ridge point, i.e., By neglecting the error term, we get and, by multiplying the second member by P(u,φ (u)) P(u,φ (u)) , we have The term P(u,φ (u)) P (u,φ (u)) only depends on the window function and its second derivative computed at the ridge point, i.e., ξ = φ (u), while P(u,ξ) P(u,φ (u)) is the value of the normalized spectrogram at point (u, ξ).
Even though previous estimation is almost acceptable whenever ξ approaches φ (u), for points far from φ (u), i.e., P(u,ξ) P(u,φ (u)) << 1, it does not provide a good approximation of the distance from the ridge.
As a matter of fact, by setting τ n = ξ n − φ (u) we can construct the following sequence that directly derives from a generalization of Equations ( 26)-( 28) written for the function Q(u, ξ) = P(u, ξ) for simplicity, where α ≤ 2 is a positive real parameter.
As proved in the following proposition, the sequence τ n monotonically converges to 0, i.e., ξ n → φ (u), and then it provides a reassignment method.Proposition 7. Let ĝ(ω) denote the modulus of the Fourier transform of a real, positive and symmetric analysis window and let Q(u, ξ) denote the square root of the spectrogram of a monocomponent function; then, for fixed u the sequence |Q (u,φ (u))| and R(τ) = Q(u,τ+φ (u)) Q(u,φ (u)) , monotonically converges to 0.
Equation ( 29) defines an iterative spectrogram reassignment method and it will be denoted as FIRST in the sequel.For monocomponent signals it represents a different way of reassigning time-frequency points and it exactly reassigns on the ridge curve as the standard reassignment does.On the contrary, in case of multicomponent signals it allows us to correctly reassign points in the interference region without changing the analysis window; hence, it contributes to reduce the region of incorrectly reassigned points thanks to the use of a weak separability condition.As a result, only a limited set of points is reassigned and each of them is forced to converge by sequential decreasing shifts to the ridge point of the mode it belongs to, i.e., the one with a dominant contribution at the considered location.
Although the better robustness to the non-separability condition, there exists a region where the weak separability condition in Definition 2 is not satisfied too.As a consequence, the shift cannot be correctly reassigned as it depends on the value of the spectrogram P(u, ξ n ) of the single mode which is unknown whenever the point ξ n does not satisfy the weak separability condition.However, the lower P(u, ξ n ) the higher the shift and the higher the probability to be far from the interference region.As a result, the larger n the higher the eventual reassignment error.Hence, in case of multicomponent signals, in correspondence to the interference region, we expect a correct reassignment for the very first iterations while convergence is not assured as n increases.
To force the convergence to the ridge point of each component, we can proceed as follows.From the contraction property it holds As a result, in case of multicomponent signals, ∀ n ≥ 1 we can modify Equation ( 29) as follows with τ 0 the initial point.The use of K instead of the spectrogram value for determining the elements of the sequence, except for the first one, makes the shift toward the ridge more accurate-as it will be shown in the next section (each point is forced to move as if it is a point satisfying the separability condition with respect to another mode).The reassignment as in Equation (33) will be denoted as FIRST_K.
Remark 3. K characterizes the analysis window ( τ satisfies Equation (32)) but its computation could be not straightforward-nonlinear iterative methods should be required.Hence, by considering that the closer , a more simple reassignment can be performed by simply setting K = 1 − α 2 and using the first few iterations of the resulting sequence, as it will be shown in the experimental results.

Computational Complexity
This section aims at estimating the number of operations per pixel required by the proposed reassignment method.Only multiplications, divisions, and comparisons will be considered.The estimated complexity is then compared with the one required by the standard reassignment.To this aim we will focus on the number of operations required by the computation of the centroid for each point of the time-frequency plane as it represents the actual difference with respect to the standard method.Equation (33) (FIRST_K) will be considered.
L denotes the signal length, L × F is the time-frequency dimension of the STFT, N is the number of adopted iterations and The estimation of the centroid using FIRST_K requires: • the computation of one STFT, i.e., (2 F log 2 (F) + F) for each signal sample; • F log 2 (F) comparisons of each u for the construction of the set Ω; • 1 multiplication for each iteration and each point in the time-frequency domain belonging to the set Ω and 1 square root for the first iteration.
Hence, the proposed reassignment method requires C proposed = |Ω|N + |Ω|β + (2 F log 2 (F) + F)L + L F log 2 (F) operations, where β is the number of operations required for the computation of the square root.
In the standard reassignment, the estimation of the centroid requires: • 2 multiplications and 1 division for each equation in (11) and ( 12) and each point in the time-frequency domain; • the computation of three STFTs, i.e., 3(2 F log 2 (F) + F) for each signal sample.
accurate than FIRST in correspondence to interference region since the sequence in (33) converges but its limit is not the ridge point of the signal mode.Results do not change if FIRST_K is used by setting , independently of the adopted analysis window-as discussed in Remark 3.However, in this case the sequence converges to a point which can be closer to the ridge point.
Figure 9 shows reassignment results provided by the proposed iterative procedure for the signals in Figures 1 and 2 where the number of iterations N has been adaptively estimated for each ratio Q(u,ξ) Q(u,φ (u)) and requiring a tolerance ε corresponding to the discretization step of the frequency axis adopted in the numerical computation of the STFT, i.e., N = log Figure 10 depicts the convergence rate for two different points in the spectrogram (ρ = 0.10 and ρ = 0.50-see Equation ( 34)).As it can be observed, for the signal with two separated modes, the proposed method is able to correctly reassign each point since the value of the spectrogram at τ 0 is exact.For the signal with two non-separated modes, the convergence is slightly altered (especially for points closer to the ridge, i.e., the ones for which the separability condition does not hold).It is due to the error on the spectrogram value caused by the interference of the two modes.In this case, FIRST_K provides a more robust solution since it only depends on the analysis window while it does not depend on the value of the spectrogram except for the initial point.
Finally, it is worth observing that the proposed method results more computationally efficient than the standard reassignment as it just requires the computation of one spectrogram; in addition, for points with the same spectrogram value, the sequence τ n in Equation ( 33) is the same-this allows us to drastically reduce the number of operations and the computing time of reassignment.

Conclusions
In this paper, an iterative spectrogram reassignment method for multicomponent signals has been presented.The method relies on a generalization of the separability condition in the case of multicomponent signals and then it mainly consists of successive shifts of selected points in the spectrogram.These points are those satisfying the separability condition with respect to signal modes, i.e., the ones less influenced by the interaction of signal modes in the time-frequency plane.As a result these points are moved as if the spectrogram is the one of a monocomponent signal.The sequential shifts are the ones estimated from a local Taylor approximation of the analysis window.Convergence of the method has been proved and discussions concerning its conditioning have been included.The latter allows us to define a robust and faster reassignment method which only requires the computation of one spectrogram.Extensive numerical studies show that the proposed method is more robust than standard reassignment method in the case of multicomponent signals with non-separable modes with constant amplitudes, especially in regions of the time-frequency plane where the modes strongly interfere.The convergence is achieved after two or three iterations making the method computationally advantageous.Future research will be devoted to an in-depth numerical study for further improving reassignment in the interference region.Particular attention will be devoted to the robustness to incorrect initial values and/or noisy data as well as to the extension to non-constant amplitude modes.

Proposition 4 .
The characteristic curves C c,φ of the partial differential equation in Equation (15) are

Proposition 5 .
The spectrogram is a constant function along each characteristic curve.Proof of Proposition 5.The proof straightforwardly follows from Equation(21).

Figure 5 .Proposition 6 .
Figure 5. Section at a fixed u in the time-frequency plane of the spectrogram of two different multicomponent signals-u = 200 (a), u = 210 (b).u belongs to an interference region.As it can be observed, there are several points that satisfy the separability condition even in the interference region-i.e., points where signal spectrogram (solid line) coincides with the spectrogram of one single mode (dotted line).

.
As it can be observed, the use of a different number of iterations for points with a different spectrogram value, allows us to correctly reassign points outside and inside the interference region.The same Figure shows the reassignment result which has been obtained using one iteration of the sequence Equation (33); results are almost equivalent, making FIRST_K a robust and computable version of FIRST in the case of lack of both weak and strong separability condition.

Figure 9 .Figure 10 .
Figure 9. (a) Reassigned spectrogram of the signal in Figure 1 using FIRST, where the number of iterations depends on the value of the spectrogram at the initial point.(b) Reassigned spectrogram of the same signal using FIRST_K.(c) Reassigned spectrogram of the signal in Figure 2 using FIRST, where the number of iterations depends on the value of the spectrogram at the initial point.(d) Reassigned spectrogram of the same signal using FIRST_K.As it can be observed, except for the central part, points in the interference regions are correctly reassigned even by FIRST_K, which is a less refined iterative method.

2. The Reassignment Method: Theoretical Model Definition 1. The
STFT at time u ∈ R and frequency ξ