Exact Short-Time Identification of Rational or Polynomial Exponent Signals

This paper describes a novel model and short-time estimation method of nonstationary or generalized sinusoids and their parameters. The model expresses the time derivative of log-amplitude and phase (complex exponent) with a rational function of time. When all poles of the rational function are simple, it is integrated to obtain the sum of polynomial and logarithm functions. The former generates a polynomial exponent signal, which is multiplied by an intensely time-varying function generated from the latter. Its direct estimator based on the weighted integral method provides an exact solution in the noiseless case from a small number of finite (short-time) Fourier coefficients; thus, multiple sinusoids separated in either the time domain or the frequency domain can be estimated independently. Several experimental tests of the basic performances are shown under possible application scenarios including the analog or digital AM/FM wave demodulation, radar/sonar pulse detection and parameterization, and combined uses of the proposed method with pulse compression.


I. INTRODUCTION
M ODELING and estimation of nonstationary sinusoids are one of the oldest and still continuing research subjects in a wide range of signal processing fields including radar, sonar, communications, acoustics, and optics. Signals as the target of detection and analysis are often characterized by their own deterministic natures. Although there are advantages as well as disadvantages in parametric approaches, the aim of model-based approaches is to extract more widely and accurately these natures in order to understand and exploit the underlying phenomena using a shorter analysis window and an increased signal-to-noise ratio (SNR). Signals in communications and audio applications are complex mixture of amplitude-modulated (AM) and frequency-modulated (FM) components with rapidly varying amplitude and phase, which are sometimes too important to be neglected for their efficient, detailed, and fruitful analysis.
Polynomial phase (including log-amplitude) signal (PPS) models are the most fundamental and extensively studied ones. The estimation methods include recursive multiplication [1], phase unwrapping and linear regression [2], polynomial phase transform [3], [4], generalized ambiguity function [5], nonlinear instantaneous least squares estimate [6], cubic or higher-order phase functions [7], [8], [9], [10], and distribution derivative method for closed form solutions [11], [12]. The PPSs, however, are not always sufficient to describe sudden transitions of real signals such as speech and music. Extensions of sinusoidal models have been explored including those of cyclostationary polynomial signals with random amplitude [13], the polynomial amplitude complex exponential (PACE) as the general solution of autoregressive model [14], gamma-enveloped sinusoids [15], the phase vocoder approach [16], and the coupled PPS and sinusoidal FM models [17], [18]. Models should be simple, mathematically tractable, and general for describing diversified signal components and for estimating parameters easily and accurately. Further preferable and unified models in this sense are still open issues. Signals for radars, sonars, and diagnostic equipment are so designed to extract delay and Doppler information under constraints in signal power and varying propagation conditions. Theoretical analysis of the waveforms is based on the ambiguity functions. Chirp signals or PPSs with suitable envelopes are commonly used. For the decoupled detection of position and velocity, linear period modulation (LPM), also known as the hyperbolic chirp/FM, has been widely studied. Studies include those on the LPM as a Doppler-insensitive waveform [19], resemblance of it to bat signals [20], quadrature detection of range and velocity [21], use of multiple waveform sets for delay Doppler imaging [22], estimation of the product of hyperbolic FM and chirp signals [23], generalized analysis of coupled and uncoupled waveforms [24], and use of sub-band correlators [25]. Most of these methods rely on the matched filter array that requires the maximum search, but the interpolation in the finite search intervals often reduces the accuracies. An appropriate use of a model-based, direct estimation method will be desirable to systematically construct suitable waveforms and to extend their overall performances.
In recent years, a mathematical technique, the weighted integral method (WIM), for differential/difference model parameter identification from finite duration observation has been developed and applied to sinusoidal parameter estimation [26], [27], [28], localization of poles of meromorphic function [29], white light interferometry [30], [31], sound source localization [32], magnetic dipole localization [33], and algebraic solution of optical flow [34], [35], [36], [37]. The WIM relies on mostly linear signal generation models; thus, it is free from iterative nonlinear optimization owing to functional nonlinearities of signal waveforms. The length of the observation interval can This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ be arbitrarily chosen without any loss of exactness, and the edge effects at the boundaries are eliminated using model-originated extra unknowns or a properly designed windowing of model equations. The purpose of this study, which is based on the WIM, is to develop a simple but extended model of nonstationary sinusoids and the direct algebraic estimation method of its parameters. In the following sections, a differential equation (DE) model is introduced, its basic properties are described briefly, and the direct algebraic estimators of its coefficients are obtained. Then, experimental tests of the basic performances are shown under several application scenarios including the AM/FM demodulation, radar/sonar pulse detection and parameterization, mixture detection, and combined uses of the proposed method with pulse compression.

A. Nonstationary Sinusoidal Differential Equation (NSDE)
Let f (t) be a complex continuous signal that is nonzero in the interval of interest. Then, where r(t) is a complex continuous function with its derivativė r(t). By multiplying both sides of (1) by f (t), and introducing a complex driving term ξ(t), the nonstationary sinusoidal differential equation is expressed aṡ with the general solution where ξ(t) = 0 and f (t) = 0 as where A is a complex amplitude, and r(t) and r(t) are the time-varying log-amplitude and phase, respectively. The driving term ξ(t) is introduced to avoid a meaningless solution f (t) = f (t) = 0 and to express a deviation from the ideal waveform due to noise. The physical meaning and existence condition of ξ(t) as a driving term will be a subject of further consideration. In (3), the redundancies in A and r(t) are removed by assuming r(0) = 0 at t = 0.

B. NSDE for Rational Exponent Signals
As the most vital part of the proposed models, consider the case where the differential exponentṙ(t) is a rational (including a polynomial) function of t expressed aṡ where a 0 , a 1 , . . . , a K and b 0 , b 1 , . . . , b L are complex constants. Without loss of generality, we can assume a 0 = 1. Under this modeling, the NSDE is expressed as which is the rational exponent differential equation (REDE). Although the integral of (4) includes log or arctan functions generally, we call the general solutions of (5) the rational exponent signals. If all poles z i (i = 1, . . . , K) of the rational function of (4) are simple,ṙ(t) can be decomposed as where P (t) is the polynomial of order less than L − K, and β i (i = 1, . . . , K) are complex constants. After integration, the signal has the form which can be seen as a polynomial exponent signal (exp{ } term) modulated by a complex time-varying amplitude (product term). The product term in (7) or the fraction terms in (4) can describe localized and rapid changes of amplitude and phase near the poles. On the basis of the Weierstrass approximation theorem, the polynomial term with a sufficient order can approximate continuous distributions uniformly with a required accuracy. Hence, it is suitable for the parametric description of arbitrary waveforms in a finite interval. The rational modeling using (4) extends this capability with a smaller number of coefficients. Also, the rational exponent is a common property of various pulsed waves in radars, sonars, and diagnosis/inspection equipment. The product, time shift, and their time shrinkage/stretch are also rational exponent signals. Therefore, estimation of the coefficients of (4) is useful for the detection and analysis of those waves in a wide range of applications.

III. IDENTIFICATION OF RATIONAL EXPONENT SIGNALS
The problem here is the estimation of coefficients a i (i = 1, 2, . . . , K) and b i (i = 0, 1, . . . , L) from the observed signal f (t) in a finite observation interval [−T /2, T/2] so that (5) is satisfied in the least-squares sense. The observation interval is involved in an isolated duration of the signal.

A. Algebraic Equations of Fourier Coefficients
By introducing finite duration Fourier bases {e −jnΔ ω t }, where Δ ω ≡ 2π/T , as a complete set of weight functions, we can equivalently express the condition for satisfying (5) Let us express the weighted integrals of f (t) as Note that no smooth window is introduced for removing edge effects. Then, by the integral by parts, the weighted integrals oḟ f (t) are expressed as where are integral boundary terms. Fortunately, these terms do not dependent on the frequency order n and the power i of t except for the parity. Therefore, by denoting we obtain the equivalent algebraic equations as where η n ≡ is the Fourier coefficient of the driving term. In (11), g n (i = 0, 1, . . . , max(K, L)) are obtained from the wave data, a 0 = 1, and a 1 , a 2 , . . . , a K , b 0 , b 1 , . . . , b L are unknowns, but F 0 involves instantaneous quantities which are not on the sample times of the sampled data. η n will be zero when the waveform has sufficient purity in the observation interval. The role of η n in this case is in the description of observation noise in g [i] n . The use of hypothesized Schwartz functions as test functions [11] or as the window function of localized Fourier transforms [12] will eliminate F + and F − . It introduces, however, considerable changes in all g [i] n . The WIM provides the elimination schemes and resultant equations according to the DE model of each problem.

B. Extra Unknown Method
In this method, we treat all unknowns including F 0 equally and simultaneously, and obtain them using a sufficient number of equations with different n so that the sum of squares of the driving term η n is minimized. Assuming a 0 = 1, it is thus expressed as a linear optimization problem where n indicates the summation among all n used in the minimization. The number of complex unknowns for the minimization is K + L + 2, whereas the equation is complex. Hence, the required number of equations is larger than or equal to K + L + 2. Not all n are necessary. This is very convenient in applications. By restricting the use of g n contributing to the waveform of interest, we can reduce interferences from different components as well as the computing costs. We refer to the extra unknown (XU) method using N F equations as the N F × XU method.

C. Windowing of Equation Method
Among the unknowns in (11), F 0 depends on the boundary values and noises of the waveform, but others do not. Firstly eliminating F 0 with the rapidly varying properties will facilitate the stabilization of succeeding numerical algorithms. The elimination of F 0 is equivalent to the relaxation of the original weight function e −jnΔ ω t into p(t)e −jnΔ ω t , where p(t) is the window function that satisfies p(±T /2) = 0 at boundaries. The necessary condition for (8) is expressed as Since the integral boundary terms vanish, it follows that n ≡g n , By assuming a 0 = 1, we can express the minimization criterion of the driving term as where n ≡ u n , and n is the summation among all n used in the minimization. The problem is linear; thus, it has a closed-form solution with the normal equation ⎡ To solve Eq. (17), a number of equations larger than or equal to K + L + 1 are necessary. Note that, for calculatingg [i] n and h [i] n , two or more Fourier coefficients are required. As a class of window functions, we use [26] which is the weighted sum of M -tuple basis functions then, it follows that This shows thatg n are obtained by convolutions with weights P m and mP m (−[(M − 1)/2] ≤ m ≤ [M/2]) on the Fourier coefficients without windowing. An example of p(t) for M = 2 (P 0 = P 1 = 1) is a half-duration cosine function which eliminates F 0 by adding order n and n + 1 equations.
using a weighted sum of order n − 1, n, and n + 1 equations with weights 1/2, 1, and 1/2, respectively. We refer to the windowing methods with M = 2 and M = 3 as the 2FW and 3FW methods, respectively. By increasing M so that the higher-order differentials of p(t) are zero at the boundaries, we can obtain an enhanced interference reduction capability from adjacent sinusoids at a slight expense of statistical efficiency [26] and also with a close relationship with the distribution derivative methods [11], [12]. n should be chosen around a peak ofg n distribution. The 2FW and 3FW methods using N F frequency orders are denoted as the N F × 2FW and N F × 3FW methods, respectively. Remark: The statistical efficiency (smallness of error variance) is best in the maximum likelihood (ML) method, but it requires an iteration and is easily trapped into local minima. The proposed method, although it is not best in the ML sense, provides directly an algebraic solution close to the global minimum. A refined solution closer to the ML, if necessary, will be obtained shortly by starting the ML iteration from the proposed estimate.

D. Desired Conditions for Stable Solution
A general theory is difficult to obtain. Thus, we describe here some requirements for the stability of algorithms.
Let n 1 , n 2 , . . . , n N F be all n used in (16). Then, so that b 0 , b 1 , . . . , b L are definite, the vectors g Then, it follows generally that hence, the weighted integral has a form that which implies that from 0th-to Lth-order differentials of the signal spectrum at frequencies n 1 Δ ω , n 2 Δ ω , . . . , n N F Δ ω should have sufficiently rich variations for definiteness. This seems to be possible when the frequencies are chosen from around the largest peak in an extent of large variations of the spectral distribution of a time-varying sinusoid. When f (t), which is modeled as a rational exponent signal, actually has a polynomial exponent in the observation interval, there exist c 0 , c 1 , . . . , c M such thaṫ Then, it follows for all n that n (0 ≤ i ≤ L); hence, the estimation matrix for (16) degenerates. This wave is, however, expressed truly by the rational exponent model when a 1 = · · · = a K = 0 and b 0 = c 0 , . . . , b M = c M , b M +1 = · · · = b L = 0. Therefore, the use of suitable regularization for the minimum norm solution will provide such estimates and avoid the ill condition. An example is shown in Example 2 of Section IV-A.

Let the integral of rational exponent be
where r(t) here is calculated using the estimated a i (0 ≤ i ≤ K), b i (0 ≤ i ≤ L), the partial fraction decomposition, and integrals such that r(0) = 0. To estimate the complex constant A of f (t) = Ae r(t) , we evaluate numerically the Fourier coefficients of e r(t) as Using the results, we express the least-squares method aŝ with the solutionÂ For the windowing of equation method, R n and g n should be those calculated with the window p(t).

IV. APPLICATION ISSUES AND EXPERIMENTAL TESTS
In this section, the proposed method is described, assuming particular conditions of various signal processing applications. Detailed analysis or comparison with existing methods is not the aim here. The applications include 1) AM-FM wave demodulation of analog or digital signal waveforms and 2) detection and parametrization of polynomial or rational exponent pulses. In the following experiments, sampled data with the length N (integer) were used. All integrals from the sample sequence were approximated using the "half sample shift two-sample extension" method described in [26].

A. Moving Estimate of Modulation Parameters
By the truncated Taylor expansion at t = 0, the exponent ζ(t) + jφ(t) of time-varying sinusoid is expressed aṡ provides an approximate of (i + 1)th order differentials of the log-amplitude and phase at the center of the observation interval [−T /2, T/2]. To estimate them, the normal equation is composed of a block for b 0 , b 1 , . . . , b L in (17). For AM-only or FM-only signals, the solution can be constrained to be real or imaginary. The number and set of n chosen should be larger than or equal to L + 1 or (L + 1)/2 so that the spectral peak and the lobe of the modulated wave are enclosed. Example 1. Demodulation of FM wave: Parameters to be estimated are b 0 (instantaneous frequency), b 1 (frequency rate), and others at t = 0 (center of observation interval). The amplitude is assumed to be constant.
The number of n for summation N F must be N F ≥ (L + 1)/2. Fig. 1 shows a numerical simulation result obtained using the 3rd order model and 2FW method (T = N = 64, N F = 5).
The FM wave (complex) is modulated by a constant frequency sinusoid with a linearly increasing amplitude, and 20% white Gaussian noise is added. In all coefficients, the rms estimation error indicated by the error bar is about three times larger than the CRLB (whole data usage in the observation interval). The relative error increases in proportion to the coefficient order. No dependences are seen on the frequency of its changes.

Example 2. Demodulation of digital QAM wave:
The quadrature amplitude modulated (QAM) wave of a digital code sequence differs significantly from those of continuous analog signals. The wave segment for each code is a pure sinusoid, whereas steep discontinuities are present at their boundaries. Detection of code boundaries is important for the synchronous recovery of the carrier and digital codes.
Description of such a discontinuity requires higher-order terms. The use of a rational exponent model is beneficial in this respect since, for |a 1 t| < 1, it equivalently has a polynomial exponent which involves an infinite number of higher-order terms. The log-amplitude rate (LAR) and instantaneous frequency (IF) are expressed directly by b 0 . The higher-order terms compose a geometrical series with the common ratio and coefficient of −a 1 and b 1 − a 1 b 0 , respectively. The estimates of a 1 , b 0 , and b 1 are provided by (17) with K = L = 1. Fig. 2 shows a result for a 16QAM wave obtained using the rational-exponent 2FW method (T = N = 64, N F = 5). In the bottom graph of a raw waveform (complex, 2% white Gaussian noise), periodic discontinuities are evident between codes. The estimates of b 0 (LAR) and b 0 (IF) are mostly zero and constant, respectively. Sharp antisymmetric spikes are evident in a 1 and b 1 at code boundaries. The top graphs show the demodulated I/Q components and the determinant of estimation matrix for a 1 , b 0 , and b 1 . The determinant reduces significantly in the code segment because the pure sinusoidal waveform there causes the degeneration. For recovering the I/Q components, the linearly increasing phase of the carrier must be suitably restored across the observation intervals.

B. Parameterization of Nonstationary Sinusoidal Pulses
The nonstationary sinusoidal pulses are widely used in radar, sonar, and measurement systems for detecting various conditions during its propagation. Most typically, they undergo the propagation delay τ and the time scaling s owing to the Doppler effect. Using (6), we can express such a wave as which is still the rational exponent signal with the same order (K, L). It also shows increased and reduced dependences on s of the polynomial and the fraction terms, respectively. τ and s will be obtainable by comparing the known and estimated coefficients. The cross-correlation between 2nd-order polynomial exponent signals such as the Gabor pulse and chirped Gaussian pulse (CGP) is also the 2nd-order polynomial exponent signal. Therefore, the proposed estimation can be used after pulse compression. These imply an extended range of applications of the parameter estimation approach. The detection/parametrization will begin with the preliminary search and localization of the target pulse.
Example 3. Moving estimation of LPM wave: This wave is also called the hyperbolic chirp/FM. Let the time origin and phase coefficient be t 0 and β, respectively. Then, the wave f (t) = Ae jφ(t) = Ae jβ log(t−t 0 ) has the "period" 2π/φ(t) = 2π(t − t 0 )/β, which is proportional to time from the origin. It is thus time-scale-invariant except for its phase constant. The REDE for this wave is expressed as where (5). The normal equation (real) for minimizing the driving term is obtained as n |u [1] n | 2 n u Fig. 3 shows a result obtained using the 2FW method (N = 128, N F = 3). In the LPM signal (complex, 10% white Gaussian noise in both real and imaginary parts), the sign of β is reversed at t < t 0 so that the phase rotation is positive. Except for the central interval with an excessively high frequency, the estimates of a 1 and b 1 and the calculated t 0 (time to the LPM origin from the center of observation interval) and β are appropriate. The rms errors of t 0 and β are about three times larger than those of the CRLB indicated by the gray zone around the true value.
Example 4. Estimation of LPM with Gaussian envelope: The Gaussian-enveloped LPM (GLPM) is used for the ranging and velocimetry of reflective objects. Let the parameters of the LPM phase be t 0 and β, and the origin and variance of the Gaussian envelope be t 1 t 0 and σ 2 , respectively. Then, the waveform is expressed as thus, the time differential of the exponent has the forṁ which corresponds to the case in (4) when The number of unknowns has increased from four to five for constructing the functional form of the rational exponent.  When all the waveform parameters are unknown, they are obtained using the relations 4 shows a moving estimation of GLPM with the 2FW method (N = 128, N F = 7). The middle low graphs are the estimated values of a 1 = −1/t 0 , b 0 = t 1 /σ 2 , and b 0 = −β/t 0 . The error variances are smallest for b 0 and b 0 (0th-order coefficients), and slightly larger for a 1 (1st-order coefficient). Since a 1 and b 0 are proportional to the coefficient β, which is time-scale-invariant, t 0 = β/ b 0 provides a better estimate when β is known. The top graph shows the calculated results of GLPM parameters when all of them are unknown. Localization of the waveform origin t 1 is possible using t 0 plus a known offset t 1 − t 0 (phase-based estimate) or the zero cross of t 1 (amplitude-based). The estimates of t 0 and t 1 have similar accuracies under this condition of waveform. The bottom graph is an example of a noisy wave (complex, origin is at the envelope peak), the middle low and high graphs are the estimated values of a 1 , b 0 , b 0 , and b 1 , b2, respectively, and the top graphs are the estimated t 0 , t 1 (left-axis scale) and α ≡ 1/2σ 2 and β (right axis scale, α is magnified 10 6 times). The vertical bars indicate the rms of estimation errors, and the thin gray zones indicate the CRLB.
They will vary according to the actual values, prior knowledge of GLPM parameters, and the presence/absence of a time scale change. After the moving estimates, extending the observation interval can further increase the accuracy.
In Figs. 5 and 6, the estimation error variances of t 0 = −1/ a 1 (LPM origin) at the waveform center t = t 1 = 0 are plotted in relation to the noise variance and the number of frequencies for the estimate, respectively. Fig. 5 shows that the 2FW method is stable under all conditions and has higher accuracy except in cases with small data length (N = 64 and N = 128). Its error variances are about 5 times larger than those of the CRLB in almost all cases. Note that it uses only N F N frequencies although the CRLB assumes all data use in the observation interval. Differences among the three methods are caused by their different schemes of eliminating the integral boundary term. The extra unknown method shows higher accuracies for short data, but it shows an accuracy saturation under high-SNR conditions. Slightly increased error variances in the 3FW method are due to its smoother windowing.   Fig. 6 shows that the estimation error variances reach their minima at small N F s and then increase gradually with the increase in N F s. A method to determine an appropriate N F is the calculation of the signal bandwidth estimated from the span of IFs in the observation interval as Adding 4 to BW/Δ ω for including the spectral sidelobe will be adequate for N F . The values are indicated in Fig. 6.

C. Decision of Detection Based on Estimated Parameters
Usually, the target pulse is detected by matched filters including a possible variation of waveforms in templates followed by the maximum selection and thresholding. The threshold is determined from an expected noise level. In the proposed method, these operations correspond to the comparisons of 1) estimated parameters with their predetermined ranges and 2) the residual error of the estimated waveform with the noise level. For 2), the residual error in (24) expressed as (30) can be compared with the expected value 2N F T σ 2 f , where 2σ 2 f and 2T σ 2 f are the noise variance of f and g n , respectively, and N F is the number of n used in the minimization.

D. Use Under Multipath Condition With Pulse Compression
The sum of rational exponent signals no longer satisfies their own REDEs. Thus, its estimation requires the isolated waveforms in the frequency domain or the time domain. For a delayed sum due to multipaths or clutters, satisfaction of this condition is very difficult. The matched filter, in contrast, is linear and maintains the additivity of signals; hence, it is extremely advantageous for use under such a condition. To extend the use of the parametric approach, we consider here some techniques for 1) superposition detection of nonstationary sinusoidal pulses and 2) combined use with pulse compression for superposition reduction.
Example 5. Measures of superposition of CGPs: The CGP is widely used for ranging and Doppler velocimetry. Detection of the superposed condition is important for the reliable use of the proposed method. The CGP waveform is expressed as thus, the time-differential of exponent has a forṁ which corresponds to the case in (4) when with the center freqnency ω 0 = 2β(t 1 − t 0 ) at t = t 1 .
Possible measures to detect the superposition will be the residual error in the amplitude/phase estimation (30) and the power of higher-order coefficients beyond the expected one. Fig. 7 shows the residual error and the higher-order polynomial coefficients of the 1st CGP at the center related to the displacement and relative amplitude of the 2nd CGP. The 2FW method (N = 128, N F = 9) is used and 10% white Gaussian noise is added. The anomalous response of b 1 (envelope curvature and chirp rate) is significant in the 0 to 2σ range of delay. Both the residual error and the excess-order coefficient increase their values in this delay range and become largest at about σ. For the application of the proposed method to the superposed CGPs in this example, this shows that at least an ∼ 2σ separation between them is required.
Example 6. Detection of CGPs with pulse compression: With pulse compression, the CGP obtains an enhanced signal energy and temporal resolution from distributed frequency components in the long envelope. In addition, the cross-correlation of different CGPs is the CGP, and autocorrelation of CGP reduces to the Gabor function, both of which can be identified directly using the proposed method. Fig. 8 shows the moving parameter estimation of CGPs (up-chirp) with a multipath distortion and 20% white Gaussian noise using the 2FW method (K = 0, L = 1, N = 128, N F = 9). Two CGPs have 1.6σ separation (σ = 100) and equal amplitudes. The middle low and high graphs are the estimated values of b 0 (IF), b 1 (frequency rate), b 0 (LAR), and b 1 (log-amplitude curvature), and the top graph is the measures of superposition used in Example 6. The rms estimation errors of the estimates and the scatters of the superposition measures are indicated by the vertical bars. The CRLB is indicated by the gray zones. The results show that both the zero crosses of b 0 Fig. 8. Moving estimate of CGP with multipath distortion using 2FW method (K = 0, L = 1, N = 128, N F = 9, 20% white Gaussian noise on the wave). The bottom graph is a noisy wave (complex) of superposed CGPs with 1.6σ separation (σ = 100) and equal amplitudes. The middle low graphs are the estimated b 0 (IF, black line) and b 1 (frequency rate, red line), the middle high graphs are the estimated b 0 (log-amplitude rate, black line) and b 1 (log-amplitude curvature, red line), and the top graph is the squared residual error in the amplitude/phase estimation. The vertical bars across them indicate the rms of estimation errors or scatters, and the gray zones indicate the CRLB.
and ω 0 crosses of b 0 of the 1st and 2nd CGPs become unclear and biased, and the estimates become more anomalous near the center of two CGPs. These situations are also indicated by the increase in residual error in the amplitude/phase estimates as a superposition measure. Fig. 9 shows the results after matched filtering. The compressed wave shown in the bottom graph provides enhanced separation and SNR. As shown in middle low and high graphs, the values of b 0 (IF) and b 1 (frequency rate) become constant (ω 0 and 0, respectively), whereas those of b 0 (LAR) and b 1 (log-amplitude curvature) increase and shorten the interference at the center; hence, the determination of zero crosses and increased curvatures at these crosses can be more accurate. In addition to the noise reduction via the matched filtering, significant reduction in estimation error of waveform parameters will be evident. The resolution of pulses and the improvement of the situation are also indicated by the drops of superposition measures at the matched-filter-output peaks in the top graph of Fig. 9.

V. SUMMARY AND CONCLUSION
In this paper, the rational exponent differential equation was introduced as a nonstationary sinusoidal model. It was then applied the weighted integral method to obtain the direct algebraic method to estimate its parameters. Experimental tests were performed to evaluate the basic performances of the proposed method under several possible application-specific conditions including AM/FM demodulation, radar/sonar pulse detection and parameterization, and combined uses of the proposed method with the multipath detection and pulse compression. They all showed extended and promising performances.
As a future work, it is desired to obtain relevant theories and algorithms for reliable applications of the proposed method to actual problems. They will include some modification to decrease the error variances toward the CRLB and methods for determining suitable model orders K and L for unknown time-varying sinusoids and an adequate observation length N and a number of frequency orders N F used for the estimation.