Polynomial Phase Estimation Based on Adaptive Short-Time Fourier Transform

Polynomial phase signals (PPSs) have numerous applications in many fields including radar, sonar, geophysics, and radio communication systems. Therefore, estimation of PPS coefficients is very important. In this paper, a novel approach for PPS parameters estimation based on adaptive short-time Fourier transform (ASTFT), called the PPS-ASTFT estimator, is proposed. Using the PPS-ASTFT estimator, both one-dimensional and multi-dimensional searches and error propagation problems, which widely exist in PPSs field, are avoided. In the proposed algorithm, the instantaneous frequency (IF) is estimated by S-transform (ST), which can preserve information on signal phase and provide a variable resolution similar to the wavelet transform (WT). The width of the ASTFT analysis window is equal to the local stationary length, which is measured by the instantaneous frequency gradient (IFG). The IFG is calculated by the principal component analysis (PCA), which is robust to the noise. Moreover, to improve estimation accuracy, a refinement strategy is presented to estimate signal parameters. Since the PPS-ASTFT avoids parameter search, the proposed algorithm can be computed in a reasonable amount of time. The estimation performance, computational cost, and implementation of the PPS-ASTFT are also analyzed. The conducted numerical simulations support our theoretical results and demonstrate an excellent statistical performance of the proposed algorithm.


Introduction
In the past few decades, the non-stationary signal processing has been a very important research field [1][2][3][4], especially estimating polynomial phase signals (PPSs). The PPSs arise in many fields including radar, sonar, radio communication, biology, and geophysics systems [5][6][7]. The Mth order PPS can be described by the following model: where A and a m denote the amplitude and phase parameters, respectively, T is the signal duration, and j is the imaginary unit. Recently, numerous methods for PPS parameters estimation have been developed. The most popular methods can be roughly grouped into three classes: estimators based on phase unwrapping (PU), estimators based on phase-differentiation (PD), and estimators based on time-frequency representation (TFR). The PU-based estimators mainly include the Kitchen's unwrapping estimator [8], the SEparate-EStimate (SEES) estimator [9], and the least square unwrapping (LSU) estimator [10]. The denotes the phase of a complex number; unwrap(·) denotes unwrap operator; |·| denotes the modulus of the internal entity; · denotes the floor operator; (·) T and (·) −1 represent the transpose and inverse, respectively.

QML Estimator
In this section, a brief introduction to QML and concepts used in our method are given. Suppose H is a set of window widths.
Firstly, we estimate the IF of the signal by detecting the ridge of STFT, which is defined by: where h ∈ H andω h (n) denotes the IF when the window width is h. Secondly, we perform the polynomial regression [18,22] on the IF to coarsely estimate the coefficients {a k , k = 1, · · · , M} in the signal phase. The polynomial regression is defined by: whereâ h = [â 1,h ,â 2,h , · · ·,â M,h ] T , N h is the number of effective signal length, and X is represented by: Thirdly, the parameter refinement procedure based on the O'Shea refinement strategy [36] is employed to improve the accuracy of the obtained coarse results, and the final estimation of the phase Fourthly, we choose another window width from H and repeat the above procedure until all widths are used. Then, from all estimations obtained with different window widths, we select the best one as an output result. However, the QML needs to perform a one-dimensional search over the set of window widths to obtain the optimal result and the computation cost high.

Proposed Estimator
In this paper, a low complexity method based on the ASTFT is proposed. The proposed method avoids both one-dimensional and multi-dimensional searches, and the error propagation problems that widely exist in the PPSs field. The proposed method mainly consists of five parts: firstly, IF of PPS is coarsely estimated by S-transform; secondly, IFG is estimated by the PCA; thirdly, window width of each time instant is estimated according to the relationship between IFG and window width; fourthly, the phase parameters of PPS are coarsely estimated by adaptive STFT and polynomial regression; fifthly, O'Shea refinement strategy is used for refining on the obtained coarse results.

Instantaneous Frequency Estimation by S-Transform
It is well known that the STFT uses a fixed sliding window to determine the tiling of time-frequency plane. Since the window is invariable, it suffers from a poor time-frequency resolution.
An improved WT method uses an adaptive window whose width changes with frequency, thus a progressive resolution is provided. However, the time-scale plots produced by the WT are unsuitable for intuitive visual analysis and they suffer from a large mathematical burden. The S-transform is a combination of STFT and WT. The advantages of S-transform are that it can preserve the information on signal phase and provide a variable resolution similar to that of the WT. Moreover, by employing a scalable and variable Gaussian window and using the Fourier kernel, the S-transform provides the phase information referenced the time origin. The standard continuous S-transform of a signal x (t) is defined by: where t denotes the time, f is the frequency, w (τ − t, f ) is a scalable Gaussian window, which is represented by: and the standard deviation σ ( f ) is a function of frequency f , which is defined by: In previous equations, it can be easily seen that window width in the time domain is inversely proportional to the frequency, which further means that, in the time domain, the S-transform window is wider for lower frequency and narrower for higher frequency.
Then, the IF can be estimated by detecting the ridge of ST, which is represented by: where f inst_ST (t) denotes the IF.
The main purpose of our method is not to obtain the exact IF, namely, a small IF estimation error is considered as tolerable. During estimation of the IF, similarly to the Fourier transform, the aliasing issue influences estimation performance. In order to enhance algorithm adaptability, we can use the resolving aliasing method presented in [19].

Instantaneous Frequency Gradient Estimation
According to the above analysis, we know that a wide window should be employed when IF of the signal varies smoothly, and a narrow window should be employed when the IF varies sharply. Accordingly, the window width depends on the IFG of the signal, which is defined by: where f inst_ST (t) denotes the IFG. In [37], the difference operator is presented for estimating the IFG. For a discrete signal with a sampling interval ∆t, the discrete IFG at the kth time sampling point is represented by: where f inst_ST [k] relates to the kth IF measurement. It can be easily seen that the difference operator in Equation (12) is sensitive to the IF estimation error, especially in noisy environments because only two IF measuremets are used to evaluate the IFG. According to [33][34][35], the PCA is a standard tool in multivariate data analysis to reduce the number of dimensions, while retaining as much as possible data variation. For a time-frequency signal, there are two features including time and IF. By performing the PCA operation on the time-frequency signal, we can obtain the principal components, and the slope of the first principal component vector approximates to the IFG [38]. For example, to get the f inst_ST [k] of kth point, we utilize (2K + 1) estimated IFs, which are obtained by Section 3.1, Then, the covariance matrix of the two-dimensional measurements is given by: When we perform the eigenvalue decomposition operator on C, the eigenvalue and eigenvector are obtained as: e = e 11 e 12 e 21 e 22 and λ= [λ 1 , Suppose λ 1 is the larger eigenvalue, then the relevant eigenvector e 1 = [e 11 , e 12 ] T is the principal component vector. Thus, the IFG f inst_ST [k] is approximated by the slope of e 1 , which is defined by: The difference operator in Equation (12) can be viewed as a special case of the PCA with using only two IF measurements to calculate the IFG. Compared with the difference operator, the explained strategy can obtain better anti-noise performance by using more IF measurements (large K). It is noteworthy that, for different PPSs, the optimal K values are also different. The optimal K value is large for low order and slowly changing IF of PPS, and small for high order and quickly changing IF of PPS. If there is no prior knowledge about the PPS, K is set to 2 based on the experience. In that case, the anti-noise performance loss is about 2 dB compared to that for an optimal K value. However, the performance of the proposed method still outperforms most of the PPS algorithms. Finally, in order to improve the anti-noise ability of the proposed method, the classical polynomial regression [18,22] is applied to IFG.

Adaptive Window Width Estimation
Suppose ∆ t and ∆ ω are the time resolution and frequency resolution, respectively. According to the Heisenberg uncertainty principle [39], we have: When the value of ∆ t ∆ ω is equal to 1/2, it is called the Heisenberg box. Since the Heisenberg box is constant, time and frequency resolution cannot be maximized at the same time. Fortunately, a reasonable harmonization between time and frequency resolution by optimizing the window width using the signal characteristics is presented in [40].
In [40], for purely frequency modulated signals, the optimal window width relates to the IFG. If the purely frequency modulated signal is defined by: where ϕ (t) denotes phase function. Then, the theoretical optimal window width T t can be calculated by: where f inst_ST (t) is the IFG.
In the proposed method, the STFT is defined by: where ω σ (t) is the Gaussian window kernel function, which is defined by: where σ is the dilation parameter. According to [41], the Gaussian window width is defined as its full width at a half maximum value (FWHM) value, which is represented by: Thus, T G = 2T t . By substituting Equation (21) into T G = 2T t , we get: As a result, the actual adaptive Gaussian window width can be obtained by σ.

Coarse Estimator
We perform the ASTFT operator on the PPS by using an adaptive window width. The IF can be estimated by: Although this estimator is biased, the IF contains phase coefficients {a k , k = 1, · · · , M} information. When the classical polynomial regression [18,22] is applied to f inst_ASTFT (t), the coarsely estimated coefficients {â k , k = 1, · · · , M} can be obtained. The polynomial regression is defined by Equations (3)- (6). It is noticed thatΩ denotes angular frequency in Equation (3).
According to [18,22], the polynomial regression is a biased estimator, thus it is very difficult to derive a closed form expression for estimator variance. Due to the bias, the obtained estimation can be treated as a coarse result. Thus, we should use a fine estimation strategy to improve estimation accuracy.

Refinement Strategy
In this section, we use the O'Shea refinement strategy proposed in [36] to refine the obtained coarse results. This O'Shea refinement is divided into three steps.
Firstly, the PPS is dechirped by the coarse results obtained in Section 3.4, and then the low-pass filtering operation is performed on the residual signal. The processing is represented by: where L denotes the filter length, Q = N/L , is the floor operator, and t l = (l + Q/2) L∆t − N∆t/2, l ∈ [−Q/2, Q/2]. The filtering in this step is used to increase the SNR explained in [36].
Secondly, the unwrapping operation is performed on the phase ofx (t l ), thus: Then, V is modeled as a new PPS in the noise with phase parameters a v = [a 0 δa 1 · · · δa M ] T , where δa k = a k −â k , k = 1, · · · , M. Then, these new phase parameters are estimated according to: where G is the rectangular Vandermonde matrix, which is represented as: Thirdly, the final estimates of phase parameters are obtained using the following relations: The PPS-ASTFT method can be described by the steps given in Algorithm 1.

Algorithm 1
Input: The PPS signal is x (t).
Step 1: Estimate PPS IF f inst_ST by detecting the range of S-transform.
Step 2: Perform Equations (13) and (14) on f inst_ST to obtain the eigenvector e. Then, the IFG f inst_ST is estimated by Equation (15).
Step 3: The theoretical optimal window width T t is determined by Equation (18). Using the relationship given in Equation (22), σ = Tt √ 2 ln 2 , the actual Gaussian window width for each time instant is estimated.
Step 4: Based on the adaptive Gaussian window width for each time instant, instantaneous frequency f inst_ASTFT can be estimated by Equation (23). Then, the coarsely estimated results {â k , k = 1, · · · , M} are obtained by Equation (3).
In addition, when the order of PPS is unknown, we can use the strategy presented in [20]. For this case, we should search over the set of potential signal orders. The algorithm can be divided into three steps. Firstly, suppose M is the set of potential signal orders, we perform Step 1 to Step 5 in Algorithm 1 on the PPS for each order m ∈ M, and we can obtain the corresponding estimation coefficients â f inal k,m , 1 ≤ k ≤ m . Secondly, evaluate the quasi maximum likelihood function presented in [18][19][20], which is defined by: whereâ f inal k,m denotes the kth (1 ≤ k ≤ m) order coefficient. Thirdly, choose the value of signal order, which maximizes the quasi maximum likelihood function, as the final estimation of the PPS order, and the corresponding coefficients are the final estimation of the PPS coefficients.

Numerical Simulations
In this study, we consider a signal of the following form: where v (t) denotes the additive Gaussian noise with SNR ∈ [−5, 25] dB, T is the signal duration, and M = 5, 6 and 7. In order to prove the effectiveness of the proposed algorithm, the parameters selection strategy is similar to that in [19]. The amplitude and initial phase are A = 1, a 0 = 20, respectively. The length of signals is N = 256, and sampling interval is ∆t = 1/256. The phase parameters are given in Table 1. According to the analysis in [13,17,18], the PHAF and the PHAF-CPF show much better identifiability and noise rejection capability than the other PD-based methods. Moreover, the QML outperforms other TFR-based methods. Therefore, we chose these three methods as references.

Simulation Time
The main implementation procedures of the proposed method include the ST operator, the PCA operator, the STFT operator and the O'Shea refinement strategy. Assume the search space is N. Due to the computational efficiency of the Fast Fourier transform (FFT) used to calculate the ST, the total number of operations is O (N + Nlog 2 N) [42]. The PCA function requires O N (2K + 1)F 2 + F 3 , where F is the number of features. In this paper, the features contain only time and frequency, thus F is equal to two. Since the matrices X T X −1 X TΩ in Equation (3) and G T G −1 G T V in Equation (26) can be calculated in advance and stored in memory, the most complicated procedure in parameter estimation is the STFT evaluation with O N 2 log 2 N operations. Thus, the total cost of the proposed method is O N 2 log 2 N . Table 2 lists the average simulation time of PHAF, PHAF-CPF, QML, and PPS-ASTFT. In this paper, the number of lag sets of the PHAF and PHAF-CPF is set to 3. According to [43], the number of QML windows is set to N/8. The simulations are completed on a computer with an Intel Core G2030 (3.00 GHz) and 12 GB memory. In Table 2, it can be seen that the CPU time of PPS-ASTFT is significantly less than that of QML, and close to that of PHAF-CPF. Among these methods, the PHAF is the fastest one. It is noteworthy that simulation time of both PHAF and PCPF-HAF increases with the signal order. According to the analysis in [13,17], as the signal order increases, more PD operations are used, and those lead to increasing of calculation.

Performance Analysis
According to the metrics used for other estimation algorithms [18][19][20]22,36], we utilize the mean-squared error (MSE) defined as follow to evaluate accuracy of the proposed algorithm: where N trials denotes the number of trials, and here N trials = 500, andâ k,r is the estimated value of a k .
According to the analysis in [18], the nonlinear techniques for higher order PPSs produce results above the CRLB for a couple of dBs. Therefore, in order to prove effectiveness of the proposed method, the same refinement strategy was used for both PHAF and PCPF-HAF to improve their accuracies. In Figure 1, the MSEs of the four highest order characteristic parameters of 5th order PPS are shown, wherein it can be seen that the TFR-based estimators, QML and PPS-ASTFT, outperform the PD-based estimators, PHAF and PCPF-HAF. Based on the analyses in [13,17], the PD-based estimator suffers from error propagation and PD operation, which lead to a high SNR threshold. The SNR thresholds of PHAF and PCPF-HAF are 8 dB and 6 dB for M = 5, respectively. Compared to the PHAF, the PCPF-HAF performs less PD operation, thus its SNR threshold is 2 dB lower. Although the simulation time of PHAF and PCPF-HAF of 5th order PPS is less than that of proposed method, the SNR thresholds are much higher. In the TFR-based estimator, all parameters are estimated simultaneously, therefore, there is no error propagation problem. Moreover, as explained in [36], the filtering tends to cause the SNR to increase, and the O'Shea refinement strategy tends to improve estimation accuracy. Hence, the QML and PPS-ASTFT can obtain better estimation performance than the PHAF and PCPF-HAF methods. In Figure 1, the SNR thresholds of QML and PPS-ASTFT are 0 dB and 2 dB, respectively, and they are much lower than those of PHAF and PCPF-HAF. Although the SNR threshold of the proposed method is slightly higher than that of the QML method, the simulation time of the proposed method given in Table 1 is much less because the QML method needs to repeat the searching and estimating operation, as described in Section 2.  Figures 2 and 3 show the same conclusions. In addition, it is noteworthy that the SNR threshold of both PHAF and PCPF-HAF becomes higher with the increase of signal order. According to the analysis in [13,17,19], with the signal order increasing, the error propagation and PD operation affect the performance more seriously. Moreover, since simulation time of both PHAF and PCPF-HAF increases with the signal order, as seen in Table 1, therefore the performance of PHAF and PCPF-HAF is worse when the PPS order is higher. However, it is easily seen form Table 1 and Figures 1-3, the performance of both PPS-ASTFT and QML is almost unaffected by the signal order. According to the analysis on computational cost and estimation performance, the proposed method is more suitable for practical applications than other methods.

Conclusions
In this paper, a low complexity method named PPS-ASTFT is proposed to estimate the PPS parameters. In the proposed algorithm, the IF is estimated by S-transform, which can preserve information on signal phase and provide a variable resolution. The width of the ASTFT analysis window is calculated by the PCA, which is robust to the noise. In the parameters estimation, by using the refinement strategy, the estimation accuracy is improved. Through the simulations, it is shown that the proposed method costs a little more time than PCPF-HAF, but the SNR threshold of proposed method is significantly lower than that of PCPF-HAF. Although the PHAF costs less time than the proposed method, the performance of PHAF is worse. In addition, the SNR threshold of the proposed method is slightly higher than that of the QML method, but the simulation time is much less than that of QML because a one-dimensional search over the STFT window widths is avoided. Therefore, through comprehensive analyses, the PPS-ASTFT is more suitable for practical applications than other methods. In the future, our work will be mainly focused on the follow issues: (i) since the S-transform suffers from poor energy concentration in some situations, while recent development in this field has enhanced the energy concentration, but the computational cost is still unsatisfactory, thus we will focus on enhancing the energy concentration in the time-frequency domain with a low complexity; (ii) since the value of K of the PCA affects the estimation performance, we will study the relationship between that value and estimation performance; and, (iii) although the proposed method is an adaptive method, the adaptability is not great. To improve the adaptability, we will study the adaptive method, such as the method presented in [29].