A pde-Based Analysis of the Spectrogram Image for Instantaneous Frequency Estimation

: Instantaneous frequency (IF) is a fundamental feature in multicomponent signals analysis and its estimation is required in many practical applications. This goal can be successfully reached for well separated components, while it still is an open problem in case of interfering modes. Most of the methods addressing this issue are parametric, that is, they apply to a speciﬁc IF class. Alternative approaches consist of non-parametric time ﬁltering-based procedures, which do not show robustness to destructive interference—the most critical scenario in crossing modes. In this paper, a method for IF curves estimation is proposed. The case of amplitude and frequency modulated two-component signals is addressed by introducing a spectrogram time-frequency evolution law, whose coefﬁcients depend on signal IFs time derivatives, that is, the chirp rates. The problem is then turned into the resolution of a two-dimensional linear system which provides signal chirp rates; IF curves are then obtained by a simple integration. The method is non-parametric and it results quite robust to destructive interference. An estimate of the estimation error, as well as a numerical study concerning method sensitivity and robustness to noise are also provided in the paper.


Introduction
Instantaneous frequency (IF) estimation is of great interest in many applications dealing with non-stationary signals, such as radar and Micro Doppler systems [1][2][3], seismic signals [4], some kinds of gravitational waves [5,6], audio [7] and human speech signals [8], animal sounds [9] and biomedical signals [10]. Practical applications deal with multicomponent signals (MCS), that is, the superposition of individual waveforms, characterized by specific time-dependent frequency content, that is, IF, as well as timedependent amplitude, namely the instantaneous amplitude (IA) [11]. The correct analysis of MCS requires signal modes separation, that also allows for IFs estimation. Conversely, in many decomposition schemes, IFs estimation is a required step for the recovery of the individual modes by the observed mixture.
It is worth pointing out that IFs detection becomes a very challenging task if signal modes interfere with each other. That is why most of IF estimators in the literature are limited to non interfering components. The existing strategies addressing this issue can be grouped into two approaches: the former analyzes the target signal directly in the time domain, such as Empirical Mode Decomposition and its improvements [12][13][14], or in the frequency domain [15,16]. As an example, de-chirping technique aims at removing the "non-stationary term" of a signal so that a narrowband filter can be used for extracting the target component [17][18][19]. In brief, the first approach consists of less or more advanced filtering procedures, performed in the time or the frequency domain. Conversely, the second approach analyzes the signal in the joint time-frequency (TF) plane, that is, it is based on TF analysis.
Even though methods belonging to the first class do not suffer from eventual distortions introduced by the adopted device (TF transform), they commonly show some sensitivity to noise and result computationally demanding.
For the same reason, in some methods TFD is regarded as an image, whose peaks are tracked and their local connectivity is exploited to estimate IF [34,35]. This is also the basic idea of Viterbi algorithm-based methods [36,37]. As main advantage, the formulation of these methods is simple but they are generally computationally expensive and sensitive to noise. Furthermore, peaks tracking-based procedures suffer from the so-called switch problem, that is, the assignment of a detected ridge point to the wrong component. It limits their application to well separated modes (i.e., non interfering).
In this framework, some recently published works address this issue by the introduction of a generalized Viterbi algorithm. The latter provides an optimization step aimed at promoting the direction of the estimated IF curve [38,39], or its local monotonicity [40]. A similar peaks-detection and tracking procedure can be found in [41]. Following the same idea, Ridge Path Regrouping Method (RPRM) offers a less expensive alternative to Viterbi algorithm, although more sensitive to noise [42,43]. It is worth noticing that interference can highly affect the observed peaks' position, resulting in incorrect estimation. For this reason, more advanced procedures try to compensate for interference effects by means of some optimization techniques [32,42].
A powerful IF estimator suitable for wideband non linearly FM signals has been proposed in [44]. In this work, the decomposition problem (and thus IFs estimation) is turned into best signal demodulation (i.e., the removal of signal non-stationarity), by minimizing the bandwidth of the demodulated signal. This method does not require TF analysis (it belongs to the first class of approaches) and it is suitable for close or even crossing modes. In addition, it is non parametric, that is, not limited to a specific signal class, as commonly done in the literature, where IFs separation is achieved by estimating the parameters of a predetermined IF model [3,9,18,26,[45][46][47][48][49][50][51][52][53][54]. Since parametric methods usually involve multi-dimensional searches in the parameter space, they result computationally expensive, that is why a non-parametric approach is more advantageous. The method introduced in [44] is very efficient in dealing with crossing modes that interfere with each other in additive way, that is, when the observed energy of the mixture increases. Unfortunately, it shows some sensitivity to the presence of destructive interference, that is, when the observed energy decreases where the two modes overlap.
It is worth observing that the above-mentioned limitation is shared by many others approaches, such as reassignment method [55], synchrosqueezing [56,57] and syncroextracting transform [58,59]. This happens because the greater the loss of resolution in a TFD, the worst the accuracy. Methods introduced in [60,61] aim at overcoming this limitation. The same shortcoming is expected in RPRM, as it assumes that the observed peaks in signal TFD reflect IF shape as a whole. Unfortunately, this is not always the case. It is also worth pointing out that interference kind (additive or destructive) depends on signal modes characteristics. The latter can be more or less emphasized in a TFD, but it remains signal dependent. As a matter of fact, the variational method proposed in [44], although TF analysis-free, suffers from the same limitation, as it will be shown later.
This paper aims at proposing an effective IF detector suitable for AM-FM twocomponent signals with close, overlapped and even crossing modes. In this last critical case, a certain robustness to destructive interference is achieved. First, a spectrogram evolution law, whose coefficients depend on signal amplitudes and phases derivatives, is introduced. The latter allows for the definition of a two-dimensional linear system, whose solution provides IFs derivatives, namely the chirp rates. Finally, IFs are achieved by integration, up to a constant.
Method sensitivity to the involved points needed to define the linear system, as well as robustness to moderate noise, are studied. A comparison to the methods introduced in [27,44] is also provided. The presented method is designed for AM-FM signals and it does not require additional assumptions concerning the signal class, resulting non parametric. It can be used as initialization of the method in [44], as well as a guide line for all peaks tracking-based procedures addressing overlapping components [42]. Finally, the proposed method may be extended to MCS with N > 2 modes if the interference regions are known in advance, as it will be discussed in the sequel.
The paper is organized as follows. Section 2 presents the mathematical background useful to introduce the problem and the proposed method; some experimental results and comparative studies are provided in Section 3 and, finally, Section 4 draws the conclusions.

Spectrogram of AM-FM Signals
According to [11], AM-FM MCS are defined by Equation (1): where N is the number of components, f k ∈ L 2 (R) denotes the k−th mode, characterized by its phase φ k , its amplitude a k (IA) and phase time-derivatives φ k (t), that is, mode IF. Signal phases, IAs and IFs are assumed to be smooth time-varying functions. The Short-Time Fourier Transform (STFT) of f , with respect to a real and symmetric analysis window g ∈ L 2 (R), is Spectrogram is defined as STFT squared modulus, that is, P(u, ξ) = |S g f (u, ξ)| 2 . According to [62], if the STFT of a single mode f (t) = a(t) cos φ(t) is computed by a normalized window g whose support is − 1 2 , 1 2 , then the STFT computed by the dilated window g s (t) = 1 √ s g t s , with s > 0, can be expressed by Equation (3): where * stands for the Fourier Transform of * and (u, ξ) denotes a corrective term, which is negligible if a(t) and φ (t) have small relative variations over the support of the window g. Under this assumption, a simple spectrogram expression can be derived from Equation (3), that is, Equation (4) shows that spectrogram reaches its maximum along the curve (u, φ (u)), u ∈ supp{ f }, defined as the ridge curve or IF curve. Equation (3) also provides a model for the spectrogram of MCS, as shown in the sequel. Without loss of generality, let us consider the spectrogram of a two-component signal f (t) = a 1 (t) f 1 (t) + a 2 (t) f 2 (t), as the one depicted in Figure 1a. By STFT linearity and taking into account Equations (3) and (4), spectrogram can be written as For a fixed u, modes are well separated if the following separability condition holds true, that is, where ∆ω is the window frequency bandwidth [62]. In this case, for each ξ belonging to P 1 support, that is, ξ : |ξ − φ 1 (u)| < ∆ω 2 , it follows |ξ − φ 2 (u)| > ∆ω 2 . Bandwidth definition thus yields toĝ(ξ − φ 2 (u)) ≈ 0 and then Equation (5) reduces to that is, mode f 1 does not affect f 2 in the TF domain, and viceversa, as in the case shown in Figure 1c.
Unfortunately, if the condition in Equation (6) is not met, the two modes have close or even overlapping supports both in time and frequency, resulting in non-negligible crossterms. In particular, for a fixed u belonging to the non-separability region, spectrogram may reach its maximum value at a frequency point different from both the ridge points of the single modes, as in the case shown in Figure 1d. For this reason, IF detectors/enhancers relying on spectrogram maxima are unreliable, especially in case of crossing IFs such that their TFDs show a lack of resolution at the non-separability region, as for the spectrogram shown in Figure 1. As it can be observed, the detected maxima in (e) do not well approximate IF curves in (b) at the region where they get closer and intersect.
If spectrogram energy increases at the non separability region, with respect to the overall energy of the two isolated contributions, that is, P 1 and P 2 , interference is said to be additive. It is referred as destructive in the opposite case. It is worth noticing that, although spectrogram is accounted for, this taxonomy depends on the cosine sign in Equation (5), that is, on signal characteristics.
Regardless of the interference effect, there are some curves other than the ridge one that convey the same information and that are less affected by cross-terms [45,60]. As proven in [61,63], the spectrogram P(u, ξ) of a monocomponent signal f (t) = a(t) cos φ(t) satisfies the following advection equation whose characteristic curves C c,φ are with c = ξ 0 − φ (u 0 ) and (u 0 , ξ 0 ) is a point in the TF plane. It is worth observing that Equation (8) is a linear advection equation, whose characteristic curves, as in Equation (9), are shifted copies of the ridge curve, which is also a characteristic curve, that is, C 0,φ . It means that IF information is not limited to the ridge curve, but it can be found in any characteristic curve. In addition, the introduced curves allow to formally characterize those points less influenced by interference, by means of the weakened separability condition (WSC), which is recalled below. Definition 1 ([60]). Two modes with IFs φ 1 (u) and φ 2 (u) are separated at time location u if there exists at least one curve in C c 1 ,φ 1 , that is, ξ 1 (u) = φ 1 (u) + c 1 , such that |ξ 1 (u) − φ 2 (u)| ≥ ∆ω; or viceversa. Remark 1. Two modes not separated at time location u in the classical sense, that is, |φ 1 (u) − φ 2 (u)| = ε ≤ ∆ω, can be weakly separated in the sense specified by Definition 1, see spectrogram external sides in Figure 1d. Indeed, by assuming, for instance, By looking again at Equation (5) and by the same argument which led to Equation (7), it follows that In case of a constant amplitude signal, that is, a (u) = 0, ∀ u ∈ supp{ f }, from Equation (8) it follows that spectrogram is a constant function along its characteristic curves. As a result, spectrogram isolevel curves lie on the characteristic curves and, operatively, any distribution thresholding provides IF, up to a constant. Figure 1f depicts some characteristic curves obtained by thresholding the spectrogram shown in Figure 1a. As it can be observed, higher levels are more influenced by cross-terms, while lower levels reveal IF shapes. Unfortunately, this argument can not be applied to time-varying signals, as the term (8) is non zero, as in the example provided in Figure 2.

The Proposed Method
This paper aims at estimating the IFs of AM-FM two-component signals. For constant amplitude signals, the latter may be estimated, up to a constant, from the characteristic curves of Equation (8), that can be easily derived from spectrogram isolevel curves. Unfortunately, in case of time-varying amplitudes, the characteristic curves of Equation (8) can not be obtained by a simple spectrogram thresholding. However, Equation (8) still is informative for our purposes. Indeed, it can be regarded as a linear equation into the It is worth noticing that the latter do not depend on variable ξ. By evaluating Equation (8) at ξ : , the ridge point, as a single mode is now considered), one can easily obtain the ratio a(u) . More in general, the following linear system can be defined: For each fixed u ∈ supp{ f }, if ξ 1 and ξ 2 are selected so that the corresponding vectors (P ξ (u, ξ j ), −P(u, ξ j )), j = 1, 2, are linearly independent, then the linear system in Equation (11) has a unique solution, whose first entry provides φ (u), that is, signal chirp rate (CR). The corresponding IF can be then estimated by integration, up to a constant.
Once the case of a single mode has been presented, let us address the two-component one in the following Proposition. Proposition 1. Let f (t) = a 1 (t) cos φ 1 (t) + a 2 (t) cos φ 2 (t) be a two-components signal and let us setĝ k =ĝ(s(ξ − φ k (u))), where g is the analysis window with length s > 0, a k = a k (u) and φ k = φ k (u), k = 1, 2, ∆φ = φ 1 − φ 2 , and ( * ) denotes the time derivative of ( * ). Then, the spectrogram P(u, ξ) satisfies the following evolution law Proof. The proof can be found in the Appendix A.
For each fixed u, Equation (12) is a non-linear equation in the unknown variables 1, 2 and k = 1, 2, that is, the signal amplitudes, phases and their derivatives. If enough frequency points are considered, the latter can be estimated by solving a non-linear system. However, this is not the strategy of the presented paper. The concept of weakened separability is exploited, instead. Precisely, if evaluated along the characteristic curves satisfying the WSC, Equation (12) reduces to Equation (11) (see Remark 1 and [60,61]), which is thus accounted for, also in the two-component case.
For the single case, by selecting two different ξ values, namely ξ 1 , ξ 2 at a fixed u, according to Equation (11), the following linear system allows for the estimation of both φ 1 (u) and where P 1ξ := ∂P 1 (u,ξ j ) ∂ξ and P 1u := A similar system can be written for φ 2 (u) and a 2 (u) . In particular, it holds By using the same equation even for two-component signals, the following estimation for φ 1 (u) would be provided where P has been considered in place of P 1 . It is worth observing that the latter assumption is valid for points (u, ξ j ) belonging to the weak separability region, as in Definition 1. Unfortunately, it does not hold true if, at a fixed time location u, that region results empty or whenever finite difference approximations are used for the partial derivatives of the spectrogram. As a result, the estimation error has to be computed in dependence on the frequency points ξ 1 and ξ 2 , as done in the next section.

Estimation Error
In this section we are interested in the estimation error for φ 1 (u) in the region Ω where even both the separability condition and the WSC do not hold. As it will be pointed out later, the same estimation still hold in the region of weakened separability whenever numerical approximation is used for spectrogram derivatives. On the contrary, as mentioned in the previous section, outside Ω the problem resembles the case of a monocomponent signal and then the error is negligible or zero.
It is worth observing that the error on the chirp rate, that is, φ 1 (u), strongly and directly depends on the error on P 1 (u, ξ) estimation through P(u, ξ), that is, ε(u, ξ) due to the presence of cross terms. In particular, the latter can present an oscillating behaviour, due to the presence of the term cos(∆φ(u)) in ε(u, ξ) definition.

Remark 2.
Spectrogram partial derivatives can be approximated from the observed data by using central finite differences, that is, where h u and h ξ denote the discretization steps, respectively in time and frequency. For each fixed u, two frequency points ξ j , j = 1, 2 are needed to define the linear system as in Equation (11). In the case of a single mode, as well as separated components, any pair of frequency points can be employed. On the contrary, a more careful selection is required for close and overlapping components. In order to prevent instabilities due to the presence of cross-terms, it is convenient to uniformly select ξ j , j = 1, 2 (i.e., they do not depend on the specific u) such that they satisfy the WSC, as in Definition 1. Taking into account Equation (10), P(u, ξ) as a function of variable ξ is close to the spectrogram section of the single component (cross-terms are negligible). As a result, the approximation of P kξ , k = 1, 2 by P ξ is more accurate, as in the case shown in Figure 3a. Conversely, it generally holds that is, P(u, ξ), as a function of variable u, does not satisfy any separability condition, as shown in Figure 3b. Then, the approximation of P ku , k = 1, 2, through finite differences is less robust (cross-terms are not negligible). It is worth pointing out that the lack of separability in Equation (20) may affect any derivative estimates, even in the WSC region and Proposition 2 still holds.

Chirp Rate Regularization
As discussed in the previous section and shown in Proposition 2, φ k (u), k = 1, 2 (signal CRs) estimation at the non-separability region is more delicate as the cross terms cannot be completely neglected. However, the estimation error can be reduced by considering the local mean of the estimated CRsφ k (u), k = 1, 2 by taking advantage of its oscillating behaviour. More precisely, CR can be defined as where s > 0 denotes the window size. On the contrary, when modes are well separated, CR estimate is more stable. Therefore, by denoting the crossing point as u cross , the following function is considered as global CR approximation, that is, that is regularized by convolution, as follows with ρ(u) = 1 √ 2π 2 exp(− u 2 2 ), > 0. Finally, IFs are estimated by integrating CR k , k = 1, 2, up to a constant. The first order rectangle quadrature scheme is adopted in this paper. The selection of the crossing point is discussed in the next section.

Crossing Point Detection
The so-called switch problem, that is, the assignment of an estimated IF point to the wrong IF curve can be prevented by the knowledge of crossing point position. The latter allows us to correctly switch modes role when estimating the chirp rates. The crossing point can be estimated by considering to characteristic curves at low level and by computing the distance between those curves for increasing time. More precisely, let us assume φ 1 (u 0 ) > φ 2 (u 0 ), where u 0 is an initial point and let us consider the set {(u, ξ) : P(u, ξ) = K, min P(u, ξ) < K < max P(u, ξ)}, that provides If signal amplitudes are slow-varying and K is sufficiently low, the term 2a (u) a(u) P(u, ξ) in Equation (8) can be neglected and c 1 , c 2 in Equations (24) and (25) do not depend on the time variable. As a result, Equations (24) and (25) indicate two characteristic curves on spectrogram lateral sides, thus less affected by eventual cross terms. Their distance is which reaches its minimum when φ 1 (u) = φ 2 (u), that is, at the crossing point.

Results and Discussion
The proposed method has been tested on several two-component signals with different frequency and amplitude modulations. This section presents the more significant results.
Synthetic signals of length n = 512 are considered and a Chebyshev analysis window of length s = 44 and side-lobe magnitude factor r ∈ [80, 100] dB is used for STFT computation, unless otherwise specified. In general, an analysis window with smooth decay to zero is recommended in order to prevent instabilities in the estimations.
As presented in the previous section, the proposed procedure allows for IF estimation, up to a constant. The latter can be derived directly from spectrogram, by selecting a peak located at the separability region. To this aim, the thresholded spectrogram is processed columnwise and the number of detected peaks, as a function of time, is computed. The initial point u 0 is chosen as the first instant providing two peaks and sufficiently apart from TF boundary, in order to prevent instabilities.
The first test is aimed at evaluating method effectiveness in case of well separated components. A two-component signal with polynomial phases and gaussian amplitudes is considered in Figure 4. The corresponding spectrogram is shown in (a) and the initial time point for integration is depicted in (b). (c) shows a generic spectrogram section and the points, on the external sides, selected to define the linear system in Equation (11). The estimated CRs are plotted in (d) and the corresponding estimated IFs can be found in (e).
Finally, (f) shows the IFs curves achieved by the classical peak detector algorithm. As it can be observed, the two methods provide similar results.
The second test, illustrated in Figure 5, concerns a FM two-component signal with intersecting modes, that interfere with each other additively. The spectrogram is shown in (a). As it can be observed, some boundary effects occur, therefore the very initial times have not been employed in the procedure, in order to prevent instabilities in the final estimate. The crossing point, detected by minimizing the distance function as in Equation (26), is depicted in (b). A generic spectrogram section is shown in (c), with emphasized points used for defining the linear system. The estimated CRs are depicted in (d) and the final result is shown in (e). As it can be noticed by comparison with (f), the proposed method outperforms the classical peak detector algorithm at the non-separability region. Figure 6 aims at evaluating the sensibility of the proposed method to the frequency points ξ 1 , ξ 2 selection. With reference to the spectrogram shown in Figure 5a, (a) emphasizes two points lying on a lower characteristic curve. This choice results in IFs estimation shown in (b). As it can be observed, the result is comparable to the one in Figure 5e, except for some instabilities at the left-hand boundary side. Indeed, the lower the characteristic the more the boundary effects. Figure 6c shows two frequency points on the lateral spectrogram sides which are closer to spectrogram maxima. The corresponding IFs estimation is provided in (d). As it can be observed, the result is less accurate. Figure 7 refers to a FM two-component signal with intersecting modes and destructive interference. As shown in (a), spectrogram resolution is significantly reduced at the nonseparability region. The distance function, as defined in Equation (26), and the crossing point estimated as its minimum are depicted in (b). A generic spectrogram section with emphasized points for defining the linear system is shown in (c). The estimated CRs are depicted in (d) and the final IFs estimation, provided by the proposed approach, is shown in (e). The latter can be compared to the result given by the classical peaks detector method in (d). As it can be observed, the proposed method better improves IF curves resolution at the non separability region.     Figure 5a. The star markers indicate spectrogram values at the frequency points ξ 1 , ξ 2 , close to the external bandwidth boundary, employed for defining the linear system, as in Equation (11), whose solution leads to IFs estimation (b). (c) The frequency points ξ 1 , ξ 2 are selected close to spectrogram maxima, giving the estimated IFs (d).   9 concern different combinations of the following modes: modulated by the amplitudes

Figure 8 refers to the signal
whose modes interfere with each other additively. The corresponding result provided by the proposed method is shown in (b). As it can be noticed, IFs estimation is quite accurate, although some instabilities occur at the right-hand boundary side. This is due to the presence of spectrogram values close to zero, which affects derivatives approximation and thus the final result. Anyway, the achieved result outperforms the one given by the peak detector (c), as well as the one provided by the Radon Spectrogram-based method introduced in [27]-it is worth observing that method in [27] is designed for constant amplitude FM signals. In Figure 9, destructive interference is addressed by considering the signal In this case, the achieved IFs estimation is less accurate, but still better than the ones given by the state-of-the-art methods shown in (c) and (d), especially if one focuses on the non-separability region, that is, close to the crossing point. Figure 10 considers different combination of the modes h 1 (t) = cos(0.5 n t 4 + 500t), h 2 (t) = cos(n(1 − t) 2 + 120(1 − t)), modulated by the amplitudes With reference to the same Figure, the left-hand side refers to while the right-hand side refers to As it can be noticed, the proposed method is able to recover IF curves at the nonseparability region. Some instabilities occur at TF boundaries, due to the presence of spectrogram value close to zero.
Some results concerning noisy signals are shown in Figure 11. The latter can be compared respectively to Figures 8b, 9c and 10d. As it can be observed, the proposed method is robust to moderate noise.   The following FM signals g 1 (t) = cos(0.5 n t 3 + 1000t) + cos(0.8 n (1 − t) 2 + 500(1 − t)) (31) g 2 (t) = cos(0.5 n t 4 + 600t) + cos(1.5 n (1 − t) 2 + 100(1 − t)), with t ∈ [0, 1], sampled at frequency n = 2000 are accounted for in Figure 12. A Chebyshev analysis window of length s = 116 and side-lobe magnitude factor r = 400 dB is used for STFT computation. Our results are compared to the ones obtained by applying the method introduced in [44,64]. As it can be observed, the proposed approach is clearly advantageous in dealing with destructive interference. Finally, Figure 13 refers to the complex signal where i denotes the imaginary unit. The signal spectrogram computed with a window of length s = 44 is shown in Figure 13a, while s = 22 has been used for Figure 13c. As it can be noticed, a narrow window in time is preferable for our purposes, as it provides a more spread spectrogram and then IFs estimation is more accurate.

1.
The proposed method requires two-component signals and prior knowledge of the presence of a non-separability region. A method for detecting the non-separability region, as the one proposed in [65], can be used as a preprocessing step, in case of constant amplitude signals. Furthermore, weakly amplitude modulation is assumed, as in most approaches.

2.
It is worth underlying that the proposed method could be easily generalized to MCS with N > 2 modes if the TF non-separability regions are known and they are sufficiently separated. Indeed, under these assumptions, the analysis of a more complex signal reduces to the analysis of a two-component signal locally. This point will be investigated in-depth in future studies. 3.
The proposed model allows for IFs estimation, up to an integration constant, whose accuracy can affect the final result, as previously shown. However, it is worth pointing out that this error often causes the recovery of a characteristic different from the ridge, but still a characteristic. As a result, it contains significant information concerning IF. The integration constant can be directly estimated from spectrogram ridges belonging to the separability region. Also in this case, a prior TF localization of the non-separability region would solve the problem. 4.
The proposed procedure requires two frequency points for defining the linear system as in Equation (11), for each fixed u. The sensitivity to their selection has been numerically investigated. As shown in Figures 5 and 6, involving too low characteristic curves can affect the final accuracy due to boundary effects, while characteristic at higher levels are generally more subjected to interference, resulting in inaccurate IFs estimate. For this reason, a good compromise can be achieved by selecting two frequency points close to the one where spectrogram concavity changes, as done in the presented simulations.

5.
As shown in the experimental results, the proposed method is robust to additive interference and cross-terms. The presence of spectrogram values close to zero, which occurs in case of strong amplitude modulation as well as destructive interference, could result in instabilities in the derivatives approximation, that can affect the final IFs estimation. However, it is worth highlighting that the proposed approach outperforms the state-of-the-art method in dealing with the critical case of destructive interference, as shown in Figures 7 and 12. In addition, the proposed method has shown robustness to moderate noise. 6.
Spectrogram is widely used in practical applications because of its simplicity and computational benefit. It is well-known that, in case of MCS, spectrogram of close modes suffers from the presence of cross-terms that can affect each estimation concerning the signal. For this reason, more advanced and adaptive kernels with attenuated cross-terms, such as S-Transform or Locally Adaptive TF distribution, are often preferred in the literature. However, it is worth noticing that many real-life measurements, such as the ones concerning human gait classification and detection, precisely deal with spectrograms. That is why spectrogram processing is still of interest, today. In addition, cross-terms do not represent a limitation for the presented method, but a tool for estimating IF.

Conclusions
This paper has proposed a non-linear time-frequency evolution law for the spectrogram of a frequency and amplitude modulated two-component signal, which is employed for instantaneous frequencies (IFs) estimation. Based on the concept of weakened separability, the latter equation can be reduced to the one referred to a single mode, whose coefficients linearly depend on signal IFs time derivatives, namely the chirp rates (CRs). Approximating spectrogram derivatives and replacing them in the evolution law yields to a linear system with signals CRs as unknown variables. As a result, CRs can be detected by simply solving a two-dimensional linear system and IFs can be estimated by numerical integration. The experimental results presented in the paper have confirmed method effectiveness in dealing with amplitude modulated two-component signals with interfering modes. A comparative study with some of the state-of-the-art methods has been also provided, showing the benefits of the proposed method in processing modes affected by destructive interference.
IF curves estimation provided by the proposed method can be adopted as starting point for all variational methods dealing with interfering modes and requiring an initialization step. In addition, since it provides CRs estimation, the introduced procedure can be useful in reallocation techniques that involve high-order phase time derivatives.
As main advantage, the presented approach is non-parametric. However, IF law may be obtained by fixing the IF class and by interpolating the result provided by the proposed method.
As future developments, the method could be employed also for instantaneous amplitudes estimation. In addition, future studies will be aimed at defining an automatic procedure for interference region detection, in order to extend the proposed method to general multicomponent signals, having more than two modes. The applicability of the presented approach to time-frequency distributions different from the spectrogram will be also investigated.

Conflicts of Interest:
The authors declare no conflict of interest.
Proof of Proposition 2. In the following, Equation (16) is proved. Equation (17) is then obtained by a straightforward integration. Without loss of generality, we set s = 1 in Equations (4) and (5).
The estimation still is correct, independently of ξ 1 and ξ 2 .