Basis approach to estimate the instantaneous frequencies in multicomponent AM-FM signals

In this paper, an analytical approach to estimate the instantaneous frequencies of a multicomponent signal is presented. A non-stationary signal composed of oscillation modes or resonances is described by a multicomponent AM-FM model. The proposed method has two main stages. At first, the signal is decomposed into its oscillation components. Afterwards, the instantaneous frequency of each component is estimated. The decomposition stage is performed through the basis expansion exploiting orthogonal rational functions in the complex plane. Orthogonal rational bases are generalized to expand linear time-varying systems. To decompose the non-stationary signal, its equivalent time-varying system is sought. The time-varying poles of this system are required to construct appropriate basis functions. An adaptive data segmentation algorithm is provided for this purpose. The effect of noise is scrutinized analytically and evaluated experimentally to verify the robustness of the new method. The performance of this method in extraction of embedded instantaneous frequencies is asserted by simulations on both synthetic data and real-world audio signal.


Introduction
Non-stationary signals which are a compound of constituents with time-varying amplitudes and frequencies can be characterized by amplitude-modulated frequencymodulated (AM-FM) models. This modeling is attended to express genuine signals in communications [1,2], acoustic and speech processing [3][4][5][6], biomedical signal processing [7], and image processing [8]. To estimate instantaneous amplitudes (IAs) and instantaneous frequencies (IFs) embedded in a multicomponent signal [9], the first step is decomposing it into oscillation components. This procedure is termed 'demodulation', 'separation' or 'decomposition' in the literature. Each component should represent a valid oscillation for which the definition of instantaneous frequency is physically meaningful [10].
All methods of multicomponent AM-FM signal decomposition and parameter estimation can be categorized into non-parametric and parametric methods. One class of *Correspondence: hamidami@aut.ac.ir 1 Electrical Engineering Department, Amirkabir University of Technology, Tehran 15914, Iran Full list of author information is available at the end of the article non-parametric approaches is based on the joint timefrequency processing. Widespread time-frequency distributions (TFDs) such as short-time Fourier transform (STFT), Wigner-Ville distribution, and Choi-Williams distribution are employed [11]. These methods are limited by the well-known compromise between time and frequency resolution. Moreover, the cross-terms appear troublesome. The energy separation algorithm (ESA) which uses a non-linear differential operator, called Teager-Kaiser energy operator (TKEO), is another method [12]. The energy separation algorithm which tracks the energy of the source producing the signal is originally applicable for monocomponent AM-FM signals. Nevertheless, they are modified for application in multicomponent cases by designing a bank of filters. Multiband-ESA (MESA) that consists of bandpass filtering followed by monocomponent energy separation is introduced based on this concept [13]. The separation performed by bandpass filtering is proper when components are spectrally far enough.
Huang et al. [14] proposed an iterative technique known as the empirical mode decomposition (EMD). This technique is an algorithmic way to extract oscillation modes embedded in the signal, named intrinsic mode functions http://asmp.eurasipjournals.com/content/2014/1/8 (IMFs). Each IMF gives a valid IF which is estimated by applying the Hilbert transform (HT) [14]. The proposed algorithm for the implementation of EMD called the sifting process showed several drawbacks such as sensitivity to perturbation and mode mixing problem [15]. To overcome these difficulties, the original EMD is modified, and a new algorithm is developed which is named ensemble empirical mode decomposition (EEMD) [15]. This method is indeed the iteration of EMD for noise-added signals. In each iteration, controlled white noise is added to the data, and EMD is applied. Each individual trial may generate noisy results, but the noise is canceled out by taking the average of the results. Thus, true solution is the ensemble mean of enough trials. Although each trial produces a set of IMFs, the sum of IMFs is not necessarily an IMF. An empirical solution for this issue is suggested in [15]. EMD and EEMD methods suffer from the lack of analytic foundation. Some research has attempted to establish and improve analytical aspects of these empirical approaches [16]. In [17], an alternative algorithm for EMD is introduced based on iterating certain filters, such as Toeplitz filters. The results of iterative filtering are similar to those of the conventional sifting process. Although the authors of [17] do not claim superiority for their method, it lays down a mathematical framework for an alternative approach to EMD. The convergence of iterative filtering EMD is studied in [18]. A variant of EMD to decompose multiscale data is proposed in [19]. This work provides some theoretical understanding of EMD for a class of multiscale data and introduces two algorithms, Newton-Raphson-based EMD and ODE-based EMD, as the variations of the sifting process. The decomposition of multiscale data based on EMD is pursued in [20] inspired by the compressed sensing theory. The sparsest representation of multiscale data is sought within the largest possible dictionary constructed of IMFs. The problem is formulated as a non-linear L 1 optimization, and an iterative algorithm is proposed to solve it. Noise and perturbation in data may cause numerical instability in this method. Daubechies et al. developed a method which captures the philosophy of EMD and decomposes special functions in a defined class [21]. This method employs a combination of wavelet analysis and reallocation technique called synchrosqueezing transform which aim to sharpen a time-frequency representation. Synchrosqueezed wavelet transform is also investigated for signal sampling and denoising applications in multicomponent signal analysis [22]. In [23], an algorithm for AM-FM parameter estimation is proposed based on the iterated application of the Hilbert transform to amplitude envelopes obtained by adaptively low-pass filters. Furthermore, the IF of AM-FM components can be calculated by a posteriori adaptive segmentation of the acquired phase signal. Another iterative AM-FM decomposition is suggested in [24] using the quasi-harmonic model (QHM) for quasi-harmonic signals such as voiced speech.
There are various parametric approaches to extract IFs of multicomponent signals. One common approach is based on signal segmentation, while some simplifying assumptions such as constant frequency seem logical in each segment. Then, an estimator is designed to estimate the model parameters segment by segment. In [25], the maximum windowed likelihood (MWL) criterion is used to estimate the AM-FM components. The high nonlinearity of this method makes the necessary optimization difficult. Another parametric approach is based on the statistical modeling of the signal according to its statistical attributes and assumptions. Speech signals are statistically modeled as AM-FM signals, and the extended Kalman filter (EKF) is applied for demodulation [3]. The idea of EKF is also exploited in [26]. Polynomial phase signal (PPS) modeling is another parametric approach which is employed for AM-FM signals [27].
The interpretations and estimation of instantaneous frequencies embedded in a multicomponent signal have been controversial [28]. Three different approaches are proposed to estimate the IFs after the decomposition stage. In the first approach, the Hilbert transform is exploited to get the analytic signal whose phase is differentiated to find the IF [10]. The energy operator (TKEO) is utilized in the second approach [12], and the third one is based on TFD [11,29]. Different definitions of IF are considered in these approaches, and consequently, their results are not necessarily equivalent. In [28] and [30], the different definitions and estimation methods are compared and discussed. The main contribution of this paper is to develop a novel approach based on the expansion of timevarying systems by orthogonal rational functions. The method introduced in [31] is extended and improved to be applied as the essence of the new method for IF estimation. An adaptive segmentation procedure in the proposed algorithm allows us to estimate the IFs locally. The decomposition is performed using orthogonal rational functions.

Problem statement
Multicomponent signals are first introduced in [9]. A multicomponent AM-FM model describes a non-stationary signal as the combination of oscillation terms with timevarying amplitudes and frequencies: where N c is the number of components. A k (t) and θ k (t) are time-varying envelope and time-varying phase of the http://asmp.eurasipjournals.com/content/2014/1/8 kth component respectively, and the instantaneous frequency denoted by f k (t) is defined from θ k (t) : The general model in (1) can be interpreted as the signal expansion by a generalized complex exponential basis, which are exponential functions with time-varying amplitudes and frequencies. The decomposition of the multicomponent AM-FM signal is investigated through this point of view. Therefore, we are going to find an appropriate basis to expand the non-stationary signal x(t): The functions {g k ; k = 1, 2, · · · , K} should represent the oscillation modes in signal, for which the instantaneous frequency is meaningfully definable. Accordingly, each term represents a valid IF of the multicomponent signal.
The main idea to attain such decomposition is expanding the corresponding system of the AM-FM signal in the complex plane. Since the transfer function of a realistic linear system has a rational representation, it can be expanded by orthogonal rational functions in the complex z-plane. Returning back to the time-domain, each rational function is equivalent to a generalized exponential basis and represents one valid oscillation term or resonance. To perform this procedure, we should specify the generating system of the AM-FM signal. The corresponding system of a non-stationary signal is modeled by a linear time-varying (LTV) system [32]. LTV models have been applied to describe non-stationary signals [33,34]. Our proposed method is developed based on this approach of modeling. Let us consider the discrete-time AM-FM signal x[n] , obtained by time sampling of x(t) at the rate of f s . Its generating system is modeled as a LTV system, which can be described by a bivariate function to characterize the input-output linear relationship [32].
where {G k (m, z); k = 1, 2, · · · , K} is a rational basis. Fortunately, to find the orthogonal rational functions for system expansion, it is not necessary to have the system's transfer function with all the details. The knowledge of poles or logical assumptions about them are sufficient to extract a proper basis [35].

Orthogonal rational functions
The decomposition step in the proposed method is indeed an expansion of the AM-FM signal, which is accomplished through the incorporating expansion of the equivalent time-varying system by orthogonal rational functions.
Generally, the knowledge about the poles is sufficient as a priori information to describe the desired space being spanned by a rational basis. Let us consider a set of M time-varying poles, {ξ k [m] , k = 1, . . . , M}. We can make a first-order IIR transfer function by each pole. So, a basis set is constructed including all specified poles, but not orthogonal. The Blaschke products [35] formed by these poles are two-dimensional functions, Applying the Gram-Schmidt procedure on these rational functions with respect to z, two-dimensional functions are obtained: This is the same routine for finding Takenaka-Malmquist functions [36]. Now, these functions are generalized to two-dimensional functions for the LTV system expansion. The resultant functions in (6) are orthogonal with respect to z in the complex domain. The inner product of each pair of these functions is a function of time, m: Utilizing the Cauchy integral implies that d kl [m] is zero at each snapshot, m, for k = l. Taking the inverse Z-transform of G k (m, z) produces g k [m, n] . Since the Ztransform and its inverse are homomorphic transforms, these functions would preserve orthogonality with respect to n. Two problems remain to be studied. At first, the underlying time-varying poles of the corresponding LTV http://asmp.eurasipjournals.com/content/2014/1/8 system should be determined. Secondly, univariate terms should be extracted from bivariate functions, g k [m, n] , to express the oscillation modes of x[n] . These issues are resolved simultaneously by adaptive segmentation which is addressed in the following section.

Time-varying modeling
The concept of poles and zeros are also generalized for linear time-varying systems. Concerning the stability and behavior of the LTV systems, several definitions of poles or eigenvalues of such systems have been proposed [37], depending on the characterization method of the LTV system. The notion of time-varying poles in this paper is founded on the time-varying autoregressive model. Parametric models for LTI systems can be generalized for LTV ones by imposing time-varying parameters on the model. The AM-FM signal, x[n] , is modeled by a time-varying autoregressive (TVAR) of order M: {a m [n] , m = 1, . . . , M} are the time-varying parameters, and ν[n] is the zero-mean innovation process, also addressed as a modeling error. The most general case of this model is where the parameters are completely uncorrelated at each time sample. Therefore, each time sample of x[n] would be represented by M unknown coefficients; hence, it is not a practical approach. Based on a common practical assumption, the non-stationary signal is approximately regarded locally stationary or quasi-stationary. This assumption implies that the parameters of the TVAR model are correlated, and the coefficients are supposed to be constant in subintervals of the total time span, referred to as segments. This model is called a block stationary AR model [34]. For multicomponent AM-FM signals whose IAs and IFs are slowly time-varying or piecewise-constant, the segmentation strategy is applicable. By virtue of this assumption, a real multicomponent AM-FM signal over its support is considered as a superposition of temporarily more limited signals with constant frequencies. These intervals can generally have various lengths, and different methods from fixed-length windowing to adaptive segmentation algorithms are introduced to determine the borders of the segments [23]. In the proposed method, the segmentation is performed adaptively from the aspect of TVAR parameter estimation.

Segmentation procedure
The entire signal of N samples is segmented into L blocks with various lengths: where = 1, . . . , L and n 0 = 0. The TVAR coefficients are supposed to be constant in each segment. The mean square error in the th segment is given by where {a ,m , m = 1, . . . , M} are the TVAR coefficients of the th segment. The boundaries of each segment are determined such that the error J remains under a specified threshold. The segmentation algorithm operates as follows. At the start of each stage, the length of the current segment (say ) is considered the minimum possible length, equal to the order of the TVAR model, i.e., n = n −1 + M. The TVAR coefficients, a ,m , are estimated by the recursive least squares (RLS) technique, and the error in (10) is computed. If it is still greater than the prespecified threshold, the length of the segment increases by one sample, and the calculations are repeated. This procedure continues by one-sample increment in each stage until the error falls below the threshold. Now, the boundaries and the length of the current segment are determined, and the procedure starts over the next time sample for another segment establishment. This algorithm runs through the entire signal repeatedly and stops at the end of the data batch. The question arises here about the threshold setting and how it can affect the accuracy of the IF estimation. This issue is scrutinized in the succeeding subsection separately. Once the TVAR parameters are estimated, the corresponding time-varying poles, denoted by {ξ k [m] , k = 1, . . . , M}, are obtained by applying the Z-transform of (8) with respect to n.

Error analysis
It is noteworthy to mention the relation between the error caused by segmentation and the error of IF estimation. This analysis leads us to select a reliable error threshold in the adaptive segmentation procedure. Let us consider a discrete-time AM-FM component: The error in the instantaneous phase imposed by the TVAR modeling in each segment, denoted by θ , sets off an error in signal: For very small phase errors, the following approximation is considered using the Maclaurin series: where The error e[n] whose instantaneous amplitude is absolute error of modeling is also an AM-FM signal: The error in phase is now transduced to the error of amplitude.
So, the absolute error of phase depends on the signal envelope. This means that for a fixed threshold, where η[n] is constant over the entire signal, larger phase errors can occur when IA becomes smaller. Therefore, the threshold should vary adaptively, adjusted to the envelope of the observed signal. In other words, the locally normalized error for each segment is a proper threshold. Since the IA evolves slowly, its mean or minimum amount during the segment can be utilized for normalization. The normalized threshold is denoted byη for brevity: Thus, the inequality (17) is practically used as the following one: The square of |ε[n]| is obtained from (15): Exploiting this relation in the inequality (19) and performing some mathematical reformulations result in a bound for the phase error: Whenη → 0, the right-hand side of the above inequality is approximately equal toη 2 . Keeping the phase error ( θ ) under control, the error of IF is consequently controlled. By definition (2), IF is the derivative of instantaneous phase, which is a difference equation in the discrete-time situation: where ω[n] = 2πf [n] is the instantaneous frequency in radian per second, and T s denotes the sampling time. In the worst case, the maximum phase errors of two consecutive instants accumulate. Thus, the maximum error of IF is 2 θ f s . For example, ifη = 10 −3 , then from (21), the maximum phase error is almost 10 −3 , and the absolute error of IF is at most 0.2% of the sampling frequency. This error can be controlled by arbitrary selection ofη. A smaller threshold leads to wider segments, in which the assumption of constant frequency is no longer respected. Our experiments verified that the condition of piecewiseconstant frequency for slowly varying IFs is satisfied forη in the order of 10 −3 ∼ 10 −2 .

Estimation framework
The main algorithm of IF estimation takes the extracted time-varying poles to construct G k (m, z) in (6) where L is the number of total segments, and  (23), which controls the effect of borderline samples. The Hamming window is commonly utilized as an analysis window in audio and speech processing [24,27]. When the signal is contaminated by noise, the timevarying poles estimated from noisy observations are http://asmp.eurasipjournals.com/content/2014/1/8 misplaced. Thus, the estimated IF incurs more error due to the error in the estimation of poles imposed by the noise. This issue is alleviated by increasing the order of the TVAR model. Each resonance of a clean signal is represented by a pair of time-varying poles; hence, the order of the TVAR model (M) is twice the number of components (N c ). Nonetheless, to improve the estimation of time-varying poles in the presence of noise, we should have M > 2N c . Therefore, extraneous poles appear besides the valid poles. A minimum distance classifier [38] is applied to assign the poles of the resonances in each segment and distinguish them from the invalid poles. The perturbation of the poles due to the estimation error of TVAR coefficients is investigated mathematically in the succeeding section. The steps of the proposed algorithm are summarized as follows: In this novel method of IF extraction, the Hilbert transform that is a global operator is not employed. Additionally, a linear model is applied to the phase of components segment by segment in spite of differentiating throughout. It makes the proposed method less sensitive to phase changes. Therefore, the adaptive segmentation is advantageous for both decomposition and frequency estimation.

Pole perturbation
The coefficients of the TVAR model (8), estimated through the RLS technique, are affected by noise. The perturbation in these coefficients leads to the perturbation in the resultant time-varying poles. Let p( , z) be the polynomial of AR model in the th segment whose roots are time-varying poles: where the coefficients are normalized, i.e., a ,0 = 1. If the perturbation a ,i occurs in the ith coefficient, a ,i , the polynomial (24) changes tõ From now on, the kth pole and its perturbation are denoted by ξ k and ξ k , respectively, and their arguments are neglected for brevity. When ξ k → 0, the above equation is simplified by employing Taylor's expansion: where p ( , z) is the first derivative of p( , z) with respect to z. This equation expresses the linear relation between perturbation in poles and perturbation in coefficients: Considering a ,i and ξ k as random variables, their variances, respectively σ 2 a ,i and σ 2 ξ k , are related linearly, Obtaining p ( , ξ k ) from (24), the ratio of the variances is given by This ratio indicates the sensitivity of the poles to the perturbation of the coefficients. Smaller ratio represents less sensitivity or more robustness of pole. When all poles are inside the unit circle, i.e., |ξ k | < 1, k = 1, . . . , M, the ratio in (30) takes its maximum value for i = M. This means that the perturbation in the coefficient of the highest order in the polynomial has the most effect on misplacing the poles. This maximum value for each pole is defined as its variance ratio: γ k is the variance ratio of ξ k which indicates the robustness of this pole. This parameter only depends on the positions of the underlying poles, which are determined by AM-FM signal characteristics.

Experimental evaluation
The proposed method is implemented on both synthetic and real-world signals, and its performance is compared to the results of two previously introduced methods. http://asmp.eurasipjournals.com/content/2014/1/8

Synthetic data
A two-component AM-FM signal is considered as below: The parameters of the first and the second components are α 1 = −2, β 1,1 = 10 3 , β 1,2 = 5 × 10 3 , β 1,3 = −20 × 10 3 and α 2 = −5, β 2,1 = 5 × 10 3 , β 2,2 = −3 × 10 3 , β 2,3 = 0. The two following IFs are embedded in this AM-FM signal: The signal in (32) is sampled at T s = 50 μs intervals, so the sampling frequency is f s = 20 kHz. Since the absolute values of the instantaneous frequencies increase over time, a limited span of signal is observed to avoid aliasing. The algorithm is run over N = 2, 000 samples of the signal. The adaptive segmentation procedure divides the data batch into L = 45 segments with different lengths. The threshold of error in the segmentation procedure is set to 0.001. The coefficients of the TVAR model and, accordingly, the time-varying poles are estimated through the RLS algorithm with the forgetting factor of 0.98. Figure 1 demonstrates the estimated IFs denoted byf i ; i = 1, 2 besides the original IFs.f 1 andf 2 track the original IFs very closely. The step-like variation, produced by segmentation, is obvious. For quantitative evaluation, the relative mean absolute error (MAE) of IF estimation is computed forf 1 andf 2 , and illustrated in Tables 1 and 2, respectively. The proposed method is compared with two previous methods, EEMD-HT and QHM. The EEMD-HT is a nonparametric method which decomposes the components by the EEMD algorithm and estimates the IFs of resultant IMFs utilizing the Hilbert transform [14]. The QHM is a parametric method which has been appraised on speech signals [24]. The result of decomposition of the original signal by the EEMD procedure is displayed in Figure 2. The relative standard deviation of added noise is 0.2, and the ensemble number for each run is 100. Nine IMFs are derived while there are just two embedded oscillation modes. The first and the second IMFs are expected oscillation modes, and the others are false IMFs generated due to deficiencies in the iterative algorithm. EEMD initially extracts faster oscillations. Thus, the first IMF has higher frequency and represents our second component, and reversely, the second IMF corresponds to the first component. The valid IMFs are selected and assigned subjectively or based on the result of estimation [15]. In this simulation, the estimated IFs of the first and the second IMFs are closer to f 2 and f 1 , respectively. MAEs of the estimation of IFs are recorded in Tables 1 and 2. The estimation errors of QHM are also provided in these tables for comparison. The analysis window of 64 samples and hop step of one sample is considered for this algorithm.
The simulation is repeated in the presence of additive white-Gaussian noise, and the errors for different SNRs are depicted in the same tables. Each value is the average of 100 iterations. The proposed method outperforms the other algorithms, especially in stronger noise. These results verify the robustness of the proposed algorithm. Although the QHM has smaller error for clean signal, it is more sensitive to noise. The EEMD, which is more robust than EMD, is still affected by the perturbations on  the amplitude of the data. Therefore, EEMD-HT has the worst performance in the presence of noise. For the clean signal (SNR = ∞), the order of the TVAR model, M, is equal to 4. Since higher error is imposed on pole estimation for stronger noise, M should increase to alleviate this problem. We have M = 24 for SNR = 30 dB and SNR = 10 dB. Figure 3 depicts the MAE of IF estimation with respect to M for both components. It demonstrates that the error decreases remarkably as the order increases. However, there are more poles for higher orders which raise the error of classification routine and consequently the total error of estimation. Thus, there is a compromise to select the optimum order. In Figure 3, MAEs off 1 and f 2 are minimum at M = 24 and M = 28, respectively. The effect of noise is illustrated by the perturbation variance ratios, γ 1 , . . . , γ 4 (31) for each segment in Figure 4. The variance ratios of the poles in a conjugate pair are very close to each other or even equal in some segments. γ 3 and γ 4 , corresponding to ξ 3 and ξ 4 or the second resonance, are smaller than γ 1 and γ 2 . This means that the poles corresponding to the first oscillation component are more sensitive to the noise. This observation justifies why the MAE forf 1 (Table 1) deteriorates more than the MAE off 2 (Table 2) due to the noise.   The proposed method can distinguish components and estimate embedded close IFs, even when they have crossovers. To illustrate this capability, the simulation of a similar two-component signal is repeated, except that the frequency of the second component, f 2 , is reduced by 4.8 kHz, i.e., β 2,1 = 1, 200 in (32). In this case, two IFs are closer and intersect each other. Both IFs are depicted in Figure 5. The proposed method is implemented in the same condition, and the results of estimation are exhibited in Figure 5, besides the original IFs. This method can track both IFs, and the intersection in their trajectories is handled properly. The degradation of estimation, especially around the crossover junction, is obvious by comparing Figures 5 and 1, but the trajectories are not missed and the estimated IFs are improved with time. This advantage arises from the step of pole classification in the proposed algorithm.

Real-world signal
Acoustic studies have revealed that many natural acoustic signals such as the sounds of songbirds and oceanic mammals fit into the AM-FM model [27,39]. A duration of 68 ms of a song of two songbirds, a canary, and a kinglet, recorded at 22.05 kHz sampling frequency, is considered for the experiment. The spectrum of this signal estimated by STFT with sliding window of 128 samples, an overlap of 100 samples, and frequency resolution bin of 2 Hz is depicted in Figure 6. This signal is composed of two AM- FM components and real-world disturbing signals such as ambient and instrumental noises. Four complex poles and consequently four bases are extracted which present two embedded dominant oscillations. Through the adaptive segmentation procedure, the signal is divided into L = 76 unequal blocks. The threshold of the segmentation procedure,η, is set to 0.01. The extracted IFs are illustrated in the same figure of the STFT spectrum by black marks. The estimated IF tracks the frequency of dominant oscillations. In Figure 7, the IFs are demonstrated in addition to the results of EEMD-HT and QHM. Nine IMFs are extracted through EEMD, but the first and the second ones are valid resonances whose IFs are displayed. The relative standard deviation of added noise is set to 0.2, and the number of ensembles is 100. It takes several samples for EEMD-HT to resolve two IFs. Furthermore, the local variations and spikes in the results of EEMD-HT and QHM are noticeable. By contrast, the IFs of the proposed method are piecewise-constant, and no spike appears due to the segmentation. The lengths of the segments and consequently the variations of the resultant IFs (flatness or volatility) are controllable by the threshold selected in the segmentation procedure.

Conclusion
The analytical approach developed for extraction of oscillation modes of multicomponent non-stationary signals is founded on several facts which are summarized herein. Firstly, the decomposition of a multicomponent signal is considered as a signal expansion and investigated from this point of view. Secondly, the equivalent timevarying system of the non-stationary signal is modeled by TVAR and expanded by the orthogonal rational functions. Finally, oscillation modes are constructed, employing these bivariate rational functions, and their IFs are estimated. For this purpose, the evolution of the time-varying  Figure 7 The results of IF estimation on the selected piece of the birds' song. http://asmp.eurasipjournals.com/content/2014/1/8 poles is considered piecewise-constant which imposes the same approximation in frequency estimation. The threshold parameter in the adaptive segmentation algorithm controls this approximation and governs the flatness of estimated IF against volatility.
The order of the TVAR model or equivalently the number of time-varying poles is another parameter which is set properly to confront the noise. Because this method utilizes the poles of the generating system of the AM-FM signal, the distortions on the amplitude, like noise, affect the estimation results mildly in comparison to those of the empirical methods. Although the noise degrades the estimation of underlying poles, it can be controlled by the order of the TVAR model. Simulations reveal the superiority of this method in the presence of noise. The controlling parameters in the proposed method yield some degree of freedom to adjust it for different realistic applications. Its capability to extract the embedded IFs in audio signals is illustrated.