A new method for parameter estimation of high-order polynomial-phase signals

Signal Processing 142 (2018) 212–222

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

A new method for parameter estimation of high-order polynomial-phase signals Runqing Cao a,b, Ming Li a,b,∗, Lei Zuo a,b, Zeyu Wang a,b, Yunlong Lu a,b a b

National Lab for Radar Signal Processing, Xidian University, Xi’an, China Collaborative innovation Center of Information Sensing and Understanding, Xidian University, Xi’an, China

a r t i c l e

i n f o

Article history: Received 28 January 2017 Revised 3 June 2017 Accepted 10 June 2017 Available online 24 June 2017 Keywords: Polynomial phase signal High-order ambiguity function (HAF) Cubic phase function (CPF) Parameter estimation

a b s t r a c t Parameter estimation of a high-order polynomial phase signal (PPS) is considered in this paper. We propose a method to estimate phase parameters more efficiently and accurately. In the proposed method, we define an operator referred to as non-uniform sampled reducing-order operator (NURO) to reduce the order of a polynomial phase by half, when the order of PPS is even. By combined using NURO and phase differentiation (PD) operators, the PPS order is reduced to one, i.e., the PPS degenerates into a complex sinusoid. Then, the parameter estimation can be done by jointly using fast Fourier transform (FFT) and one-dimensional search. Compared with the traditional methods, the reducing-order procedure in the proposed method has lower-order nonlinearities. Simulation results show that the proposed method outperforms the hybrid CPF-HAF and HAF in both mean square error (MSE) and the threshold when the PPS order is higher than 5. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Polynomial phase signal (PPS) has widely application in the parameter estimation and target detection in radar, sonar, communication, and various nature signal analysis [1–7]. To estimate the parameters of PPS, the maximum likelihood (ML) estimation was proposed in 1960s. Although it performs well in the sense of mean square error (MSE), it suffers a heavy computational burden due to the multidimensional search. Therefore, methods with low computational burden attract more and more attention. The highorder ambiguity function (HAF) [2], which is realized by a onedimensional search, develops a phase differentiation (PD) technique to reduce the PPS order. By applying PD once, PPS order can be decreased by one. However, as to high-order PPS, many times of using PD leads to high-order nonlinearities, which increases the number of noise influenced terms. Therefore, the HAF algorithm has a poor performance at low signal-to-noise ratio (SNR). To estimate the third order phase parameter, some timefrequency rate representations are proposed, for example, the cubic phase function (CPF) [3], the high resolution time-frequency rate representation (HR-TFRR) and their extensions [4,5]. By combining the CPF and HAF, the hybrid CPF-HAF is proposed in [6] to estimate parameters of high-order PPSs. Compared with the HAF, the hybrid CPF-HAF algorithm requires two PDs less which leads to rel∗

Corresponding author. E-mail address: [email protected] (M. Li).

http://dx.doi.org/10.1016/j.sigpro.2017.06.011 0165-1684/© 2017 Elsevier B.V. All rights reserved.

atively lower MSE and lower SNR threshold. Although the hybrid CPF-HAF improves the estimate performance, the order of nonlinearities grows exponentially with the increase of the PPS order. Expect the hybrid CPF-HAF algorithm, the quasi-maximumlikelihood (QML) estimator and its extensions [8–12] are proposed to improve the performance at low SNR. By utilizing the robustness of the short-time Fourier transform (STFT) to noise influence, the QML based algorithm achieves a lower SNR threshold and a lower MSE than the HAF. However, the improvement is achieved at the cost of greatly increased computational complexity. In this paper, we propose a new method to estimate the phase parameters of a high-order PPS more efficiently and accurately. In the proposed method, we define the kernel in non-uniform sampled cubic phase function (NU-CPF) [13,14] as a reducing-order operator referred to as non-uniform sampled reducing-order operator (NURO). When the PPS order is even, NURO can decrease the PPS order by half. Utilizing this property, we perform a multistep procedure combining the PD and NURO operator to reduce the order of the PPS. Instead of using PD repeatedly in the HAF based approach, this procedure selects and uses either PD or NURO operator in each stage. When the PPS order is odd, we perform PD and reduce the PPS order by one. On the other hand, if the PPS order is even, we perform NURO and reduce the PPS order by half. In this manner, by using the PD or NURO repeatedly, we can transform a high-order PPS into a complex sinusoid with fewer bilinear operations i.e., involving lower-order kernel nonlinearities. Since the frequency of the sinusoid is proportional to the highest-order phase

R. Cao et al. / Signal Processing 142 (2018) 212–222

213

parameter, the estimate of highest-order phase parameter can be done efficiently by jointly using fast Fourier transform (FFT) and one dimensional search. The proposed algorithm outperforms traditional estimate algorithms for high-order PPSs in MSE and signalto-noise (SNR) threshold.

kernel [14]. When estimating the parameters of high-order PPSs, the order of kernel nonlinearities is very high for HAF algorithm, which leads to a high SNR threshold [12,15].

2. PD and NURO operator

For third order PPS, a technique called the CPF is proposed in [3]. It is defined by

2.1. PD operator

2.2. NURO operator

CP F (n, ) =



x(n ) = s(n ) + w(n ), − N/2 ≤ n ≤  N/2  P  i s(n ) = b0 exp ( jφ (n ) ) = b0 exp j ai (n) ,

(1)

i=0

where φ (n) is the signal phase, P is the order of the polynomial in the signal phase,  is the sampling interval, w is complex zero-mean white Gaussian noise with variance σ 2 , and {b0 , a0 , a1 , a2 , ..., aP } are real-valued parameters of the PPS. The length of the sequence is N + 1, where N is even. The SNR is defined as 2 SNR = b0 /σ 2 . In this paper, we assume b0 and a0 are known constant. The problem is to estimate parameters {a1 , a2 , ..., aP } from x(n). The PD operator is used in the HAF to decrease the order of the polynomial in signal phase. It is defined in [2] as

P Dτ [s(n )] = s∗ (n − τ )s(n + τ )



j 2



= b20 exp j

P−1 

np

(P−p+1  )/2

p=0 l=1 ϕ0 (n ) + θ1 nP−1 ,

( p+2l−1 )! p+2l−1 2l−1  τ a p+2l−1 p!(2l−1 )!



(2) 2 P P τ a

where θ1 = P , ϕ 0 (n) is a polynomial including lower-power terms of n, and  ·  denotes the rounding up operation. From (2) we can see that the resulting signal phase is a polynomial with respect to n, and the order of the PPS is decreased by one. Knowing that the original coefficient of the highest-power term nP in the phase of s(n) is θ0 = P aP , the coefficient is changed by PD operator according to

θ1 = 2 P τ θ0 .

(3)

The CPF has an interpretation as a time-frequency rate representation. Thus, the third and second-order parameters of the PPS can be estimated by evaluating the function at two different instants:

aˆ3 =

arg max |CP F (n1 , )| − arg max |CP F (0, )| 

P Dτ1 [x(n )] = x (n − τ1 )x(n + τ1 ) ∗



∗

1

.. .





CP [s(n )] = b0 exp

aˆP =

2P−1 P P !

P−1

i=1

j

P 

×b0 exp

j

b20

exp

j

=

b20

exp

j

=

exp

P 

The next lower-order parameter aP−1 can be estimated by performing the same procedure on the demodulated signal x(n ) exp(− jaP P nP ). From (4) it can be seen that for a P th-order PPS, the HAF algorithm requires P − 1 times of PD, i.e., the HAF has a kernel of 2P−1 th-order nonlinearities, where the order of nonlinearities equals to the number of signal factors in the product

a i 2 i 2 ( n − m )

2ai1 

i2

P/2  l=0

j

P/2 

i1

i1 /2  



i ni1 −2l m2l 2l

l=0

P/2 P  

l=0

(5)

l=0 i1 =2l

b20

P 

i1 =0



τi

ai1  ( n + m )

i2 =0

= b20 exp 2 j

.

i1

i1

i1 =0

Then the highest-order parameter can be estimated by

 P−1 arg max P Dτ1 τ2 ...τP−1 [x(n )] exp(− jn )  n

(9)

Substituting a Pth-order PPS described in (1) into (9) yields



(4)

(8)

CP [s(n )] = s(n − m )s(n + m ).

=

P−2 P−2 P DP−1 τ1 τ2 ...τP−1 [x (n )] = P Dτ1 τ2 ...τP−2 [x (n − τP−1 )] P Dτ1 τ2 ...τP−2 ×[x(n + τP−1 )].

(7)

To improve the performance of the non-linear technique based algorithm at low SNR, the idea of reducing the order of nonlinearities is brought up in [18] and [19]. Basing on CPF and PD operator, the hybrid CPF-HAF is also proposed in [6]. This algorithm applies P − 3 times of PD to reduce the order of polynomial in phase to 3, and then CPF is used to estimate the parameters. The resulting order of kernel nonlinearities is 2P−2 . Compared with HAF, it uses two PDs less, therefore, the lower-order kernel nonlinearities leads to a better performance in the MSE and SNR threshold [6]. However, the order of kernel nonlinearities still grows exponentially with the increase of the PPS order. Thus, the SNR threshold is still very high when estimating parameters of high-order PPS. It can be seen from (6) that, in the CPF, a bilinear operator is applied to the PPS, which can be written as

1

∗

,



P Dτ1 τ2 [x(n )] = P Dτ1 [x(n − τ2 )] P Dτ1 [x(n + τ2 )] 2



2n1 

aˆ2 = arg max |CP F (0, )|.

In the HAF algorithm, the property of PD is utilized to reduce the order of the polynomial in the signal phase. The HAF algorithm can be described as follows. First, perform PD repeatedly until the PPS is transformed into a complex sinusoid: 1

(6)

k

Consider a discrete PPS defined as

= b20 exp

x(n + k )x(n − k ) exp(− jk2 ).

i1 ! a i1 ni1 −2l m2l ( i 1 − 2 l )! ( 2 l )! i 1

P 2 m2l  i1 ! a i1 ni1 −2l ( 2 l )! ( i 1 − 2 l )! i 1 i1 =2l



2 φ 2l ( n ) 2l m , ( 2 l )!

(10)

To reduce the computational complexity of CPF, the NU-CPF is proposed in [13], where the bilinear operator is modified so that the CPF can be performed by using FFT. We name the modified kernel in the NU-CPF as the NURO operator, which is given by

NUCP [s(n )] = s(n −

√ √ cm )s(n + cm ).

(11)

Here, we exploit its usage of reducing PSS order. Similar to the derivation in (10), applying the NURO operator to a P th-order PPS

214

R. Cao et al. / Signal Processing 142 (2018) 212–222

described in (1) yields



NUCP [s(n )] = b20 exp

j

P/2  l=0

2 φ 2l ( n ) (cm )l , ( 2 l )!

(12)

where  ·  denotes the rounding down operation. When P is even, we have



j 2φ (n ) + φ 2 (n )cm

NUCP [s(n )] = b20 exp +

2φ 4 ( n ) 2 2 2φ P (n ) P/2 P/2 c m + ... + c m 4! P!





= b20 exp j

 ϕ0 (m ) + θ1 mP/2 ,

(13)

where θ1 = p is the coefficient before the highest-power term mP/2 , ϕ 0 (m) is a polynomial including lower-power terms of m. We can see from (13) that the resulting signal phase is a polynomial with respect to m, and the PPS order is decreased to P/2. Knowing that in the phase of s(n) given by (1), the original coefficient of the highest-power term nP is θ0 = P aP , NURO operator changes the coefficient θ 0 according to 2P cP/2 a

θ1 = 2cP/2 θ0 .

(14)

From above we can see that different from PD operator, which reduces the order of the polynomial in the signal phase one by one, the NURO operator can reduce the order by half. That means when P is even, the NURO operator can be used to reduce the PPS order more efficiently.

Less times of using bilinear operator implies lower-order nonlinearities, which improves the performance of the estimator at low SNR [6,18,19]. With this goal, we develop a new technique to reduce the PPS order to one by combined using PD and NURO operator. The reducing order approach uniquely depends on the PPS order P, which may consist of several PD or NURO operators. Here, we use a 0–1 sequence {Q1 , Q2 , ..., QR } to record the operators we need to perform, where Qr = 0 means the r th operator is a PD operator, and Qr = 1 means the r th operator is a NURO operator. R is the total number of operators we need. Moreover, we use sequence {P1 , P2 , ..., PR } to record the PPS order every time before we perform an operator. Thus, when P > 1, the sequence and number of PD and NURO operators can be determined by following procedure. First, we judge if P is even or odd and let P1 = P . If P is odd, that means we need to perform a PD operator, so we let Q1 = 0, and P decreases by 1. If P is even, that means we need to perform a NURO operator, so we let Q1 = 1, and P decreases by half. Then, we go on judging P and perform cycle operation, until P = 1. In this manner, we can determine the values of Q1 , Q2 , ..., QR and P1 , P2 , ..., PR in turn. Then we can construct a reducing-order kernel for a Pth order PPS, which is written as

where

OPQr [x(n )] =



P D [ x ( n )]

Qr = 0

NUCP [x(n )]

Qr = 1

,

(15)

r = 1, 2, ..., R.



θR = 2 c M

M−1 

cj × 2

H

H 

Pi τi

θ 0 = 2 R P c

i=1

j=1

H 

(Pi τi )aP ,

(17)

i=1

where Pi , i = 1, ..., H is the PPS order before we apply the ith PD operator. Pi can be easily obtained from sequence {P1 , P2 , ..., PR } and {Q1 , Q2 , ..., QR }. It follows that the DFT of KP [x(n)], which can be written as

f () =

U 

KP [x(n )] × exp(− jm ),

(18)

m=0

reaches its maximum at 0 = θR = 2R P c the upper limit on m.

H

i=1

(Pi τi )aP , where U is

3. Estimate method

2.3. Reducing-order kernel consisting of PD and NURO operators

KP [x(n )] = OPQR [...OPQ2 [OPQ1 [x(n )]]...],

operator, c1 , c2 , ...cM−1 in other operators are all set to be one, because the sampling density-control factors before constant parameters {m1 , m2 , ..., mM−1 } make no sense. For the special case where P = 1, we define KP [ · ] as an identity transform, i.e., KP [x(n )] = x ( n ). The reducing-order kernel KP [ · ] can decrease the PPS order to one, i.e., transform a P th-order PPS into a complex sinusoid. As we discussed in Part 1 and 2 Section 2, applying PD and NURO operator changes the coefficient before the maximum power term in the signal phase according to (3) and (14) respectively. Therefore, after applying H and M times of PD and NURO operators, the final coefficient of the maximum power term, i.e. the frequency of sinusoid KP [x(n)] is

(16)

The reducing out kernel we obtain is uniquely determined by PPS order P. Here specific expansions of KP [x(n)] for P = 5, 6, 7, 8 are given in Table. 1, where H and M are number of PD and NURO operators, and {τ 1 , τ 2 , ...τ H }, {m1 , m2 , ...mM−1 } are constant parameters in PD and NURO operators respectively. Except for the sampling density-control factor c before variable m in the last NURO

For a given Pth-order PPS, the proposed method of estimating the parameter a1 , a2 , ..., aP is given in the form of pseudo code. Input: P, P th-order noisy PPS x(n ) While P > 0 do Choose the reducing-order kernel KP [·] according to P If P > 1, then Resample the signal at non-integer time instants we need in order to evaluate KP [x(n )] end if Evaluate KP [x(n )] Evaluate DFT of KP [x(n )] Maximize the mode of DFT of KP [x(n )] with respect to discrete frequency . max Estimate aP by aˆP = H

2 R P c

i=1

Pi τi

P := P − 1 x(n ) := x(n ) × exp(− j aˆP P nP ) End Output: {aˆ1 , aˆ2 , ..., aˆP }

The estimate algorithm can also be interpreted into three steps in details as follows: Step 1 Apply the reducing-order kernel KP [ · ] described in Part 3 Section 2 on the PPS, i.e., evaluate KP [x(n)]. When evaluating KP [x(n)], we need samples at non-integer time instants. In the proposed method, these non-uniform signal samples are obtained from original integer samples via an interpolation process. The interpolation process includes two substeps [14]: First, interpolate x(n) by a factor 2 or 4 using a standard interpolation to obtain a new signal x0 (n). Then calculate the approximate signal value at non-integer time instants n0 by

x0 ( n0 ) =

n0 − f loor (n0 ) x0 (ceil (n0 )) ceil (n0 ) − f loor (n0 ) ceil (n0 ) − n0 + x0 ( f loor (n0 )), ceil (n0 ) − f loor (n0 )

(19)

where floor(n0 ) and ceil(n0 ) are two sample times closest to n0 , ceil(n0 ) is the upper time value, while floor(n0 ) is the lower one.

R. Cao et al. / Signal Processing 142 (2018) 212–222

215

Table 1 Expressions of reducing-order kernel and parameters for P = 5,6,7,8. P

Expanded expression of reducing-order kernel

{Q1 , ..., QR }, {P 1 , ..., P R }

5

K5 [x(n )]= NUCP[NUCP [P D[x(n )] ]]   √ √ √ √ ∗ ∗ =s(n +  m1 +√ cm + τ1 )s (n + m1 +√ cm − τ1 )s(n − m1 +√ cm + τ1 )s (n − m1 +√ cm − τ1 ) s(n + m1 − cm + τ1 )s∗ (n + m1 − cm − τ1 )s(n − m1 − cm + τ1 )s∗ (n − m1 − cm − τ1 )

{0,1,1} {5,4,2}

6

K6 [x(n )] = NUCP [P D[NUCP [x(n )] ]]   √ √ √ √ + m1 + cm + τ1 )s(n  − m1 + cm + τ1 )s∗ (n  + m1 + cm − τ1 )s∗ (n  − m1 + cm − τ1 ) = s (n  √ √ √ √ s(n + m1 − cm + τ1 )s(n − m1 − cm + τ1 )s∗ (n + m1 − cm − τ1 )s∗ (n − m1 − cm − τ1 )

{1,0,1} {6,3,2}

7

= NUCP [P D[NUCP [P D[x(n )]]]] K7 [x(n )]   √ √ √ √ + m1 + cm + τ2 − τ1 )s(n − m1 + cm + τ2 + τ1 )s∗ (n  − m1 + cm + τ2 − τ1 ) = s(n + m1 + cm + τ2 + τ1 )s∗ (n  √ √ √ √ s∗ (n + m1 + cm − τ2 + τ1 )s(n + m1 + cm − τ2 − τ1 )s∗ (n − m1 + cm − τ2 + τ1 )s(n − m1 + cm − τ2 − τ1 ) √ √ √ √ s(n + m1 − cm + τ2 + τ1 )s∗ (n + m1 − cm + τ2 − τ1 )s(n − m1 − cm + τ2 + τ1 )s∗ (n − m1 − cm + τ2 − τ1 ) √ √ √ √ s∗ (n + m1 − cm − τ2 + τ1 )s(n + m1 − cm − τ2 − τ1 )s∗ (n − m1 − cm − τ2 + τ1 )s(n − m1 − cm − τ2 − τ1 )

{0,1,0,1} {7,6,3,2}

8

= NUCP[NUCP[NUCP[x(n )] ]] K8 [x(n )]       √ √ √ √ = s(n + m1 + m2 + cm)s(n − m1 + m2 + cm)s(n + m1 − m2 + cm)s(n − m1 − m2 + cm )         √ √ √ √ s(n + m1 + m2 − cm)s(n − m1 + m2 − cm )s(n + m1 − m2 − cm)s(n − m1 − m2 − cm )

{1,1,1} {8,4,2}

Step 2 Calculate the function f(), which is defined as

f () =

U 

KP [x(n )] × exp(− jm ),

(20)

τi = ωi (gU − gD )/2, ωi ∈ (0, 1).

m=0

where U is the upper limit on m. It can be seen from (20) that f() is the discrete Fourier transform (DFT) of KP [x(n)]. Thus, we calculate f() by using FFT. Since KP [x(n)] is a sinusoid, according to the analysis in Part 3 Section 2 f() reaches its maximum at

0 = 2R P c Hi=1 (Pi τi )aP . Step 3 Estimate the parameter aP by

aˆP =

1 arg max | f ()|, H 

2R P c Pi τi

0, the range of the variable n in the new PPS x (n) is updated to [gD + τ1 , gU − τ1 ]. τ 1 can be determined within the range (0, (gU − gD )/2). Therefore, let

(21)

i=1

where τ i is the lag parameter of the ith PD. The estimate method is for the highest-order parameter aP in the signal phase. Once aP is estimated, the lower-order parameter aP−1 can be estimated by performing the same procedure on demodulated signal x(n ) exp(− jaˆP P nP ). The procedure is repeated until all the parameters are estimated. Here we summarize the proposed method as follows. Instead of using PD repeatedly in the HAF based approach, the proposed method selects and uses one of PD and NURO operator depending on P in each stage. When P is odd, we perform PD and reduce the PPS order by one. On the other hand, if P is even, we perform NURO and reduce the PPS order by half. In this manner, by using PD or NURO repeatedly, one can transform a high-order PPS into a complex sinusoid with fewer times of bilinear operations, i.e., involving lower-order kernel nonlinearities. 4. Parameter definition In Section 3, in order to calculate KP [x(n)], we need to determine the values of n, c, U, lag parameters {τ 1 , τ 2 , ...τ H } and {m1 , m2 , ...mM−1 }. The reducing-order kernel KP [ · ] consists of a series of operators. Since the range of the parameter in the r th operator OPQr depends on previous ones from the first to the (r − 1 ) th operator OPQ1 , OPQ2 , ..., OPQr−1 , we need to determine the parameters from the operator OPQ1 to OPQR in turn. Let us begin from a P th-order PPS described in (1). The variable of the PPS is n. We assume the range of n is [gD , gU ], where the upper and lower limits gD and gU are interaction variable, respectively. Their initial values are gD = −N/2, and gU = N/2. If a PD operator is applied on the PPS, the renewed PPS is x (n ) = x(n + τ1 )x∗ (n − τ1 ) in terms of (2). Assuming τ 1 >

(22)

If a NURO operator is applied on the PPS, the renewed PPS √ √ is x (m1 ) = x(n + c1 m1 )x(n − c1 m1 ) in terms of (11). Assuming m1 > 0 and c1 = 1, the range of the variable m1 of the new PPS 2 x (m1 ) is updated to [0, (min(gU − n, n − gD ) ) ], where min(gU − n, n − gD ) denotes the minor one between gU − n and n − gD . Now the variable n in the original PPS should be determined as a constant parameter. In order to achieve a better estimate performance, we usually tend to sample the signal in a wider time range. Therefore, in order to make the range of variable m1 as wide as possible, let

n = (gU + gD )/2.

(23) 2

Thus the range of the new variable m1 is [0, (gU + gD ) /4]. As the sequence of the operators is known, we perform the cycle operation, renew the PPS and update the variable range. When performing PD operator, we determine τ i by (22). On the other hand, when performing NURO operator, we determine n and following mj by (23). It can be seen from (22) and (23) that these parameters n, {τ 1 , τ 2 , ...τ H } and {m1 , m2 , ...mM−1 } depend on range limits. According to the interaction procedure described above, the upper and lower range limits gU and gD of the variable in the r th operator can be calculated interactively by following function

g(r, g0 ) = fQr (· · · fQ2 ( fQ1 (g0 ))),

(24)

where gU and gD are the two elements of vector g, and vector g0 = (−N/2, N/2 ). The output of the function fQr is also a two-element vector. fQ is defined as r



f Qr ( g ) =

(gD + τu , gU − τu ) ( 0,

(gU −gD )

2

4

)

Qr = 0 Qr = 1

,

r = 1, . . . , r − 1, (25)

where τ u represents the parameters of PDs before the r th param eter OPQr . u can be calculated by u = r − rl=1 Ql . Thus, the general formula for calculating τ i and mj can be expressed as

τi =

ωi (gU (I (i ), g0 ) − gD (I (i ), g0 ) ) 2

, ωi ∈ (0, 1 ),

(26)

and

mj =

gU (J ( j ), g0 ) + gD (J ( j ), g0 ) , 2

(27)

216

R. Cao et al. / Signal Processing 142 (2018) 212–222

Fig. 1. MSE of a5 and a6 estimate with respect to U, when SNR = 30 dB (a) for a5 . (b) for a6 .

where I(i) is the index of the I th operator in operator sequence OPQ1 , OPQ2 , ..., OPQR , which is the i th PD; and J(j) is the index of the J th operator in operator sequence OPQ1 , OPQ2 , ..., OPQR , which is the j th NURO operator. These two mapping relations can be easily obtained from the 0–1 sequence {Q1 , Q2 , ..., QR } that is derived in Part 3 Section 2. The value of ωi effects the MSE of the estimator. The optimal value of ωi varies with the change of Ps. We optimize the value of ωi with respect to the theoretical MSE for P = 5, 6, and 7 in Appendix B. The result is as follows: for P = 5, ω1opt = 1/5; for P = 6, ω1opt = 1/3; for P = 7, ω1opt = 1/7, ω2opt = 1/3; The optimal values of ωi for other Ps can be obtained by following the procedure in Appendix B. Next, we will define the value of another parameter c and U. Parameter c controls the sampling density of KP [x(n)]. When the range of variable cm is constant (determined by parameters n, {τ 1 , τ 2 , ...τ H } and {m1 , m2 , ...mM−1 } we discussed above), c is inverse proportional to number of samples U. According to [14], if the signal is resampled too densely, there will be redundancy in presentation of the information. Therefore, in order to make it easy to process and analyze, the sampling density should be sparse enough so that all the samples which contribute to form the kernel are white. Therefore, it requires that the number of samples after resampling procedure not more than original number of samples 2N+ 1, which can be written as

2RU ≤ N + 1.

(28)

Under this condition, signal should be sampled as dense as possible to achieve a better MSE performance. The MSE of a5 and a6 estimation with respect to U is shown in Fig. 1, where the MSE is derived by 10 0 0 times of Monte Carlo tests, parameters {τ 1 , τ 2 , ...τ H } and {m1 , m2 , ...mM−1 } are determined by (26) and (27), N = 256, U = 4, 5, ..., 64 and the SNR is 30 dB. Fig. 1 shows that the MSE gets lower as U increases and it decreases asymptotically to a limited extend. If we continue increasing U beyond the threshold in Fig. 1, the improvement of MSE is limited (about 0.8 dB), while the cost is expensive, because a large U implies a immensely increased amount of the calculation. Therefore, we set U to be its upper limit in (28):





U = N/2R ,

(29)

and c can be determined by

c=

gU (R, d0 ) − gD (R, d0 ) . U

(30)

5. Analysis and discussion 5.1. Order of kernel nonlinearities Take a sixth-order PPS for example, i.e., P = 6. By implementing Step 1 in the estimate algorithm we obtain H = 1, M = 2, R = 3, and sequence {Q1 , Q2 , Q3 } = {1, 0, 1}. Thus, we can apply a NURO operator reducing the order of the PPS to 3, then apply a PD operator reducing the order of the PPS to 2, finally a NURO operator, making the signal a sinusoid. Thus the kernel for estimating a6 can be expressed as

K6 [s(n )] = NUC Pm [P Dτ1 [NUC Pm1 [s(n )]]]   √ √ = s∗ (n − m1 − τ1 + cm )s∗ (n + m1 − τ1 + cm )   √ √ s(n − m1 + τ1 + cm )s(n + m1 + τ1 + cm )   √ √ s∗ (n − m1 − τ1 − cm )s∗ (n + m1 − τ1 − cm )   √ √ s(n − m1 + τ1 − cm )s(n + m1 + τ1 − cm ). (31) It can be seen from (31) that the designed kernel has eighthorder nonlinearities. On the other hand, HAF for sixth-order PPSs has thirty second-order nonlinearities, and the hybrid CPF-HAF has sixteenth-order nonlinearities. Low-order nonlinearities imply less cross-term of noise and signal component, therefore the estimation is less effected by noise. So the proposed method has lower SNR threshold. To estimate parameter aP in a P th-order PPS, the proposed method requires 2R th order kernel nonlinearities, where R is an P integer satisfying R < 2log2 . On the other hand, the HAF has 2P−1 th-order nonlinearities, and the CPF-HAF has 2P−2 th − order nonlinearities. As shown in Fig. 2, compared with the HAF and the hybrid CPF-HAF algorithm, the proposed method has a lower-order nonlinearities kernel when P > 5. As P grows, the order nonlinearities used in the proposed method grows much more slowly than that in the other two algorithms. So the proposed method performs better than traditional algorithms especially for high-order PPSs. 5.2. Theoretical MSE We derive the MSE of aP estimate for P th-order PPS in Appendix A. Calculating τ i , mj , U, c by (26), (27) (29) and (30), and substituting them into (A.23), we derive the MSEs for different Ps. The results along with theoretical MSEs of the hybrid CPFHAF [6] are shown in Table 2. The calculation of Cramer–Rao low bounds (CRLB) can be found in [16]. Here we define λ as the ratio

R. Cao et al. / Signal Processing 142 (2018) 212–222

217

order of nonlinearities

600 HAF Hybrid CPF-HAF Proposed method

500 40 400 20

0 4

300

5

6

200 100 0

4

5

6

7 P

8

9

10

Fig. 2. Orders of nonlinearities of the kernel in HAF, hybrid CPF-HAF, and the proposed method.

Fig. 3. MSE of a6 estimated by the HAF, hybrid CPF-HAF and the proposed method. Table 2 Theoretical MSEs of aP estimated by the hybrid CPF-HAF and the proposed method.

λ

P

Asymptotic MSE of the hybrid CPF-HAF

Asymptotic MSE of the proposed method

CRLB

5

10.16×105 SNRN 11 10

625162 SNRN 11 10

349272 SNRN 11 10

1.6252

6

2.4763×107 SNRN 13 12

1.0848×107 SNRN 13 12

5549544 SNRN 13 12

2.2827

7

6.6078×108 SNRN 15 14

1.2502×108 SNRN 15 14

8.8339×107 SNRN 15 14

5.2854

9

1.4079×109 SNRN 17 16

3.5821

8

10

1.8948×10 SNRN 17 16

5.2897×10 SNRN 17 16

of the asymptotic MSE of the hybrid CPF-HAF to that of the proposed method. As shown in Table 2 the asymptotic MSEs of the hybrid CPF-HAF and the proposed method are all close to the CRLBs. For P = 5, 6, 7, 8, the asymptotic MSE of the hybrid CPF-HAF algorithm is respectively 1.6525, 2.2827, 5.2854, 3.5821 times as much as that of the proposed method. The proposed method has better MSE performance in comparison with the hybrid CPF-HAF algorithm. 5.3. Calculation complexity

Fig. 4. MSEs of aP estimated by HAF, hybrid CPF-HAF, QML, and the proposed method (a) P = 5, (b) P = 7, (c) P = 8.

Assuming the interpolation used in Step 1 of the estimate method is performed by zero padding in the frequency domain, the overall calculation complexity of the interpolation procedure is O(Nlog N) [13]. Calculating KP [x(n)] with the non-uniformly spaced samples needs 2N times of multiplication, whose calculation complexity is O(N). Calculating the FFT of KP [x(n)] requires Ulog2 U times of addition and (1/2)Ulog2 U times of multiplication. Knowing that U = 2N/2R , the calculation complexity of calculating the FFT of KP [x(n)] is O(Nlog N). Estimating parameters a1 , a2 ..., aP needs P times of the amount of calculation we discussed above. However, P is independent of N and P  N. Therefore, the overall calcula-

tion complexity of the proposed method is O(Nlog N). The calculation complexity of the HAF [2], the hybrid CPF-HAF [13], the nonuniform sampled CPF-HAF [13], and the QML [8] algorithm along with proposed algorithm is listed in Table 3. It can be seen that the calculation complexity of proposed method is lower than that of the hybrid CPF-HAF and the QML, and of the same level with the HAF and the non-uniform sampled CPF-HAF. It can be seen from the analysis that the proposed method achieves a lower nonlinearity and MSE when compared with other nonlinear transform based algorithms including the HAF and the

218

R. Cao et al. / Signal Processing 142 (2018) 212–222 Table 3 Calculation complexity of the HAF, hybrid CPF-HAF, the non-uniform sampled CPF-HAF, the QML, and the proposed method. Algorithm

HAF

Hybrid CPF-HAF

Non-uniform sampled CPF-HAF

Proposed method

QML

Calculation complexity

O(Nlog N)

O(N 2 )

O(Nlog N)

O(Nlog N)

O(N 3 )

Fig. 5. MSEs of parameters estimation of a 6-th order PPS by using HAF, hybrid CPF-HAF, QML, and the proposed method (a–f): MSEs of a6 –a1 estimation.

hybrid CPF-HAF algorithm. Because the estimation can be performed by using FFT, the proposed method has a much lower calculation complexity than the QML algorithm. 6. Simulation results To evaluate the proposed algorithm, first we estimate the parameter a6 of a sixth-order PPS

s6 (t ) = exp( j (1.2n + 1.3 × 10−3 (n )2 +2 × 10−4 (n )3 + 4.96 × 10−6 (n )4 +5.5 × 10−9 (n )5 + 1.5 × 10−9 (n )6 )) by the HAF, CPF-HAF, QML, and proposed algorithm respectively. The number of the signal samples is 257, n ∈ [−128, 128], and

the sampling interval  = 0.1. By using (26), (27), (29) and (30), we obtain U = 23, c = 228 /207, τ1 = 8192/3, m1 = 64. The window width in the QML is chosen to be 32n. MSEs of a6 estimation obtained by 10 0 0 times of Monte Carlo tests are shown in Fig. 3. As shown in Fig. 3, the MSE of the proposed method is close to the CRLBs above the threshold and consistent with theoretical MSE. Due to the lower-order nonlinearities, the SNR threshold of proposed method is lower than that of the hybrid CPF-HAF and HAF by about 4 dB and 14 dB, respectively. The MSE of the proposed method is about 2 dB lower than hybrid CPF-HAF and 6 dB lower than HAF. Compared with the three kinds of non-linear based estimators, the QML has better estimating performance at all SNRs.

R. Cao et al. / Signal Processing 142 (2018) 212–222

Fig. 6. MSEs of parameters estimation of a 7th-order PPS by using HAF, hybrid CPF-HAF, QML, and the proposed method (a–g): MSEs of a7 –a1 estimation.

219

220

R. Cao et al. / Signal Processing 142 (2018) 212–222

We also estimate ap of P th-order PPS, for P = 5, 7, 8. The signal expressions are



s5 (t ) = exp j 1.2n + 1.3 × 10−3 (n )2 + 2 × 10−4 (n )3 +4.96 × 10−6 (n )4 + 5.5 × 10−9 (n )5





s7 (t ) = exp j 1.2n + 1.3 × 10−3 (n )2 + 2 × 10−4 (n )3 +4.96 × 10−6 (n )4 + 5.5 × 10−9 (n )5 +1.5 × 10−9 (n )6 + 1.6 × 10−9 (n )7





s8 (t ) = exp j 1.2n + 1.3 × 10−3 (n )2 + 2 × 10−4 (n )3 +4.96 × 10−6 (n )4 +5.5 × 10−9 (n )5 +1.5 × 10−9 (n )6



+1.6 × 10−9 (n )7 + 1.6 × 10−8 (n )8 . The number of the signal samples is 257, n ∈ [−128, 128], and the sampling interval  = 0.1. MSEs of aP estimation are obtained by 500 times of Monte Carlo tests. For P = 5, U = 29, c = 234 /18125, τ1 = 128/5, m1 = 131072/5. For P = 7, U = 25, c = (9 × 232 )/60025, τ1 = 128/7, τ2 = 98304/49, m1 = 294912/49. For P = 8, U = 9, c = 250 /9, m1 = 8192, m2 = 33445532. As shown in Fig. 4, the MSEs of the proposed method for all Ps are close to the CRLB above the threshold. The MSE of the proposed method is about 3–6 dB lower than the hybrid CPF-HAF for different Ps. The threshold is about 8 dB and 18 dB lower than the hybrid CPF-HAF for P = 7 and P = 8 respectively. We can see that the advantages of proposed method in SNR threshold become greater as the order of the PPS increases. Then we do another two experiment to evaluate the performance of the proposed algorithm when estimating lower-order parameters a1 , a2 , ..., aP−1 . The PPS orders are P = 6 and 7. Signal expressions are respectively



s6 (t ) = exp j 1.1n + 1.38 × 10−3 (n )2 +5.5 × 10



(n ) + 1.5 × 10 (n ) −9

5

6



s7 (t ) = exp j 1.1n + 1.38 × 10−3 (n )2 +1.89 × 10−4 (n )3 + 5.06 × 10−6 (n )4 +5.5 × 10−9 (n )5 + 1.5 × 10−9 (n )6 +1.6 × 10−9 (n )

 7

Acknowledgment This work is supported by the National Defense Foundation of China under grant 61424010302162401002, Shanghai Aerospace Science and Technology Innovation Funds (SAST) under grant 90716160 0 09, National Natural Science Foundation of China under grant 61501342, and in part by the Fundamental Research Funds for the Central Universities (JB160213). The authors would like to thank the Associate Editor and the Reviewers for their constructive comments that helped us improve the paper. Appendix A In Appendix A, we derive the asymptotic MSE of aP estimate for a P th-order PPS by following the first order perturbation analysis in [17]. As we discussed in Section 3, aP is estimated by

aˆP = 2R



1 M

P i=1

Pi τi c

 U arg max KP [x(n )] × exp(− jm ) .  m=0

(A.1)

Here we define the function

f () =

U 

KP [s(n )] × exp(− jm ).

(A.2)

m=0

+1.89 × 10−4 (n )3 + 5.06 × 10−6 (n )4 −9

reduces the PPS order with lower-order nonlinearities. Therefore, it provides more accurate estimate and lower SNR threshold than traditional methods. Since the method performs less bilinear operations and it can be implemented by FFT, it has low computational complexity. The results are supported by theoretical analysis and confirmed by numerical examples.

.

The number of signal samples is N = 513, the sampling interval is  = 0.05. We use the HAF, the hybrid CPF-HAF, the QML, and the proposed method to estimate parameters {a1 , ..., aP }. The window width of the STFT in the QML algorithm is 32n. In Fig. 5, the performance of proposed method for lowerorder parameters a1 , a2 , ..., aP−1 is similar to the estimation of the highest-order parameter aP . The SNR threshold of the proposed method is 8 dB, which is respectively 7 dB and 15 dB lower than the hybrid CPF-HAF and the HAF algorithm, and 8 dB higher than the threshold of the QML algorithm. As to the MSE, the QML has the best performance. However, as we discussed in Part 3 Section 5, compared with other nonlinear transform based algorithms, the QML based algorithm has a greater computational complexity. Fig. 6 shows similar results with Fig. 5. For P = 7, the MSE of the proposed method is lower than those of the hybrid CPFHAF and HAF for about 0–6 dB and 0–10 dB respectively. The SNR threshold of the proposed algorithm is about 12 dB for a7 , a5 , a3 , a1 , and 8 dB for a6 , a4 , a2 . The thresholds of the proposed algorithm are lower than that of the HAF and hybrid CPF-HAF, and higher than that of the QML.

As we discussed in Part 3, Section 2, KP [s(n)] degenerated into

a complex sinusoid whose frequency is θR = 2R P c H i=1 (Pi τi )aP . R

2 Thus we have KP [s(n )] = b0 ρ × exp(2R P c H i=1 (Pi τi )aP m ), where ρ is a phase factor independent of m. Then by substituting

KP [s(n)] and 0 = 2R P c H i=1 (Pi τi )aP into (A.2), we derive

f ( 0 ) =

U 

b0

ρ × exp j2  c R

m=1

P

H 

(Pi τi )aP m

i=1

× exp(− j2R P c

H 

(Pi τi )aP m )

i=1

=

U 

b0

2R

ρ

m=1

= b0

2R

ρU.

(A.3)

Similarly, we have U ∂ f ( 0 )  2R = b0 ρ (− jm ) ∂ m=1

=

1 2R 1 2R b0 ρ (− j )(U + 1 )U ≈ b0 ρ (− j )U 2 , 2 2

(A.4)

U ∂ 2 f ( 0 )  2R = b0 ρ (−m2 ) 2 ∂ m=1

1 2R 1 2R = − b0 ρU (U + 1 )(2U + 1 ) ≈ − b0 ρU 3 , 6 3

7. Conclusion In this paper, we propose a method for parameter estimation of high-order PPSs. By combining PD and NURO operator, the method



2R

δ f ( 0 ) =

U  m=0

(KP [x(n )] − KP [s(n )]) exp(− j0 m ),

(A.5)

(A.6)

R. Cao et al. / Signal Processing 142 (2018) 212–222

we can obtain

U  ∂δ f (0 ) = −j m(KP [x(n )] − KP [s(n )]) exp(− j0 m ). ∂ m=0

(A.7)

In (A.4) and (A.5), the approximations are derived by omitting the lower power terms of U. Based on the first order perturbation analysis [17], the perturbation of the maximum position of function f() is

A B

δ = − , where

(A.8)





∂δ f ∗ (0 ) ∂ f (0 ) ∗ + δ f ( 0 ) , ∂ ∂   ∂ 2 f ∗ ( 0 ) ∂ f ( 0 ) ∂ f ∗ ( 0 ) B = 2Re f (0 ) + . ∂ ∂ ∂ 2 A = 2Re f (0 )

(A.9)

  U  1 m ∗ lim ρ − (E {KP [x(n )]} −KP [s(n )] ) exp(− j0 m ) = 0 3 2 N→∞ U 2U m=1 (A.18) Then, from (A.15) we obtain

lim E {δ } = 0

The perturbation of maximum of function f() is shown to be zero asymptotically. Because of (A.1), the aP estimate is (to first order approximation) asymptotically unbiased [19,20]. Based on the first order perturbation analysis [17], the expectation of the squared δ  is given by

(A.10) E



δ

2



 

=

E A2

(A.11)

A ≈ −2b20 U[Im{η}],

(A.12)

R

m−

m=0

U 2



× (KP [x(n )] − KP [s(n )] )∗ exp(− j0 m ).



E {ηη

δ  = − 12b2ImR U{η} 3 0  U  ∗ 12 m 1 = − 2R Im ρ ( U 3 − 2U 2 )(KP [x(n )] − KP [s(n )] ) exp(− j0 m ) . b m=1

E { δ } = −

2R

b0



Im

ρ

 U   1 m − (E {KP [x(n )]} − KP [s(n )] )∗ U3 2U 2

m=1



× exp(− j0 m )

(A.15)

KP [x(n)] is the noisy signal processed by reducing-order kernel. Obviously, |E{KP [x(n)]}| has an upper limit. We define the upper limit as L. Thus, we have

  U  1 m ∗ − ( E K [ x ( n ) ] − K [ s ( n ) ] ) exp ( − j  m ) { } ρ P P 0 3 2 U 2 U m=1   U  m 1 ∗ ≤ 3 − 2 (E {KP [x(n )]} − KP [s(n )] ) exp(− j0 m ) m=1

=

U

2U

 U   1 m 3 − 2 × |(E {KP [x(n )]} − KP [s(n )] )| U 2U

 2 η = 0.

 U   m 1 R + × (L + b20 ) U3 2U 2

m=1

Substituting (A.11)–(A.13), (A.16) and (A.17) into (A.15), we get the theoretical MSE of the parameter aP estimate for Pth-order PPS, which is

E



δ aP

2



L + b20 U +1 + 2 2U 2U

(A.16)

Knowing (A.16) and



lim

U→∞

R b20

L+ U +1 + 2U 2U 2

  E δ 2

=

(  2R

Pc

6 2

H

i=1

l=1

≈

(Pi τi ) )

2R 

U3

  2R l

(  2R

Pc

SN Rl H

i=1

2

.

(A.23)

(Pi τi ) )

Appendix B In Appendix B, we optimize parameter ωi with respect to the theoretical MSE we derived in Appendix A for P = 5, 6, 7. The MSE in (A.23) can be written as

E



δ aP

2



6

where

l=1

=



ξ = U3 c

 

2R 



ξ

H 

2R

2R l



P

SN Rl H

i=1

2 ,

(B.1)

Pi

2 τi

.

(B.2)

i=1

R

=

(A.22)

When P = 5, the form of the PPS is given by (1). By following the procedure described in Part 3 Section 2, we derive the 0–1 sequence {Q1 , Q2 , Q3 } = {0, 1, 1}. That means we need to perform a PD operator followed by two NURO operators to reduce the order of the given PPS to one. By using (26), (29), and (30) we can obtain

m=1

≤

(A.21)

2R ) denotes the binomial coefficient. l

(A.14) It follows that

 

R

l

E

Substituting (A.11)–(A.13) into (A.8) yields

12

(A.20)

2  1 2R 2R+1 −2l 2l b0 σ , } ≈ U3 l 12

∗

where (

(A.13)

0

.

B2

From (A.13), we get

1 R+1 B ≈ − U 4 b20 , 6

U  

(A.19)

N→∞

Substituting (A.3)–(A.7) into (A.9) and (A.10), we obtain

η=ρ

221

τ1 =

ω1 N 2

,

(B.3)

and

 = 0,

(A.17)

c = N 4 ( 1 − ω1 )4 / ( 26 U ).

(B.4)

222

R. Cao et al. / Signal Processing 142 (2018) 212–222

Substituting (B.3), (B.4), and (29) into (B.2) yields

ξ=

ω12 (1 − ω1 )8 N

  10 N 2R

214

.

It follows that

ω1 (1 − ω1 )7 (1 − 5ω1 )N10 ∂ξ = ∂ ω1 213

(B.5)

N 8

.

(B.6)

To minimize the MSE, we need to maximize ξ with respect to ∂ξ ∂ ω = 0. Knowing that ω 1 ∈ (0, 1), we obtain

ω1 . Therefore, we let ω1opt = 15 .

1

Following the same procedure, we obtain ω1opt = 13 for P = 6. When P = 7, there are two parameters ω1 and ω2 need to be optimized. By (26), (29), and (30) we obtain τ1 = ω1 N/2, τ2 = N2 ω2 (1 − ω1 )2 /8, c = N4 (1 − ω2 )2 (1 − ω1 )4 /(64U ). Substituting τ 1 , τ 2 , c and (29) into (B.2) yields

ξ=

N 16

N 14

220

ω12 (1 − ω1 )12 ω2 2 (1 − ω2 )4 .

∂ξ Let ∂∂ξ ω1 = ∂ ω2 = 0. Then we obtain ω1opt = P = 7.

(B.7) 1 7

and ω2opt =

1 3

for

References [1] B. Porat, Digital Processing of Random signals: Theory and Methods, Prentice-Hall, Englewood Cliffs, NJ, 1994. [2] B. Porat, B. Friedlander, Asymptotic statistical analysis of the high order ambiguity function for parameter estimation of polynomial-phase signals, IEEE Trans. Inf. Theory 42 (3) (1996) 995–1001. [3] P. O’Shea, A new technique for estimating instantaneous frequency rate, IEEE Signal Process. Lett. 9 (2002) 251–252. [4] L. Zuo, M. Li, Z. Liu, L. Ma, A high resolution time-frequency rate representation and the cross term suppression, IEEE Trans. Signal Process. 64 (10) (2016) 2463–2474.

[5] L. Zuo, M. Li, Xiang-Gen Xia, New smoothed time-frequency rate representations for suppressing cross terms, IEEE Trans. Signal Process. 65 (3) (2017) 733–747. [6] I. Djurovic, M. Simeunovic, S. Djukanovic, P. Wang, ‘A hybrid CPF–HAF estimation of polynomial-phase signals: detailed statistical analysis’, IEEE Trans. Signal Process. 60 (10) (2012) 5010–5023. [7] Y. Gao, H. Li, B. Himed, Knowledge-aided range-spread target detection for distributed MIMO radar in non-homogeneous environments, IEEE Trans. Signal Process. 65 (3) (2017) 617–627. [8] I. Djurovic, L. Stankovic, Quasi-maximum-likelihood estimator of polynomial phase signals, IET Signal Process 8 (4) (2014) 347–359. [9] I. Djurovic, L. Stankovic, STFT-based estimator of polynomial phase signals, Signal Process. 92 (11) (2012) 2769–2774. [10] I. Djurovic, On parameters of the QML PPS estimator, Signal Process. 116 (2015) 1–6. [11] I. Djurovic, High precision technique for PPS estimation in impulsive noise environment, Signal Process. 127 (2016) 151–155. [12] I. Djurovic, QML-RANSAC: PPS and FM signals estimation in heavy noise environments, Signal Process. 130 (2017) 142–151. [13] M. Simeunovic, I. Djurovic, Non-uniform sampled cubic phase function, Signal Process. 101 (2014) 99–103. [14] P. O’Shea, Improving polynomial phase parameter estimation by using nonuniformly spaced signal sample methods, IEEE Trans. Signal Process. 60 (7) (2012) 3405–3414. [15] S. Barbarossa, A. Scaglione, G.B. Giannakis, Product high-order ambiguity function for multicomponent polynomial-phase signal modeling, IEEE Trans. Signal Process. 46 (3) (1998) 691–708. [16] B. Ristic, B. Boashash, Comments on ‘The Cramer-Rao lower bounds for signals with constant amplitude and polynomial phase, IEEE Trans. Signal Process. 46 (6) (1998) 1708–1709. [17] S. Peleg, B. Porat, Linear FM signal parameter estimation from discrete-time observations, IEEE Trans. Aerosp. Electron. Syst. 27 (1991) 607–616. [18] I. Djurovic´ , M. Simeunovic´ , Parameter estimation of non-uniform sampled polynomial-phase signals using the HOCPF-WD, Signal Process. 106 (January) (2015) 253–258. [19] P. O’Shea, A fast algorithm for estimating the parameters of a quadratic FM signal, IEEE Trans. Signal Process. 52 (2) (2004) 385–393. [20] P. Wang, I. Djurovic, J. Yang, ‘Generalized high-order phase function for parameter estimation of polynomial phase signal’, IEEE Trans. Signal Process. 56 (7) (2008) 3023–3028.