- Email: [email protected]
Non-uniform sampled cubic phase function
Signal Processing 101 (2014) 99–103
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Fast communication
Non-uniform sampled cubic phase function Marko Simeunović n, Igor Djurović University of Montenegro, Electrical Engineering Department, Cetinjski put bb, 81 000 Podgorica, Montenegro
a r t i c l e i n f o
abstract
Article history: Received 14 December 2013 Received in revised form 24 January 2014 Accepted 5 February 2014 Available online 13 February 2014
Parameter estimation of polynomial phase signals (PPSs) based on the cubic phase function (CPF) and its extensions cannot be performed by using the fast Fourier transform (FT) algorithm. Therefore, in order to express the CPF by means of the FT, in this paper we propose a scheme for the CPF evaluation based on non-uniform sampling. Calculation complexity of the estimation procedure is significantly reduced, whereas the accuracy is the same or better compared to the original algorithm. & 2014 Elsevier B.V. All rights reserved.
Keywords: Polynomial-phase signal Cubic phase function Fourier transform
1. Introduction The cubic phase function (CPF) has been proposed for the estimation of third order polynomial phase signals (PPSs) [1]. Several generalizations of the CPF to PPSs of higher orders exist in the literature, including the recently proposed hybrid CPF–high-order ambiguity function (CPF–HAF) [2,3]. However, since the CPF-based techniques cannot be evaluated using the fast Fourier transform (FFT), they are characterized by larger computational complexity than the FT-based techniques such as the HAF and product HAF (PHAF) [4,5]. With the requirement of OðN2 Þ complex multiplications and additions, where N is the number of signal samples, the CPF can be unfeasible in real-time applications, for example, in radar and sonar applications [3]. A procedure for parameter estimation of higher order PPSs based on the non-uniform signal sampling has been proposed in [6,7]. By sampling the signal at nonequidistant time instants, the procedure lowers the nonlinearity of the estimator function, which, in turn, improves its performance, signal to noise ratio (SNR) threshold and mean squared error (MSE). In this paper, we use a similar non-uniform sampling scheme to
n
Corresponding author. E-mail address: [email protected] (M. Simeunović).
http://dx.doi.org/10.1016/j.sigpro.2014.02.005 0165-1684 & 2014 Elsevier B.V. All rights reserved.
transform the CPF in a form suitable for the FFT algorithm. The modified CPF requires OðN log NÞ operations and has lower MSE than the standard CPF. It can be generalized for other CPF-based techniques. Here, we implement the modified CPF in the evaluation of the CPF–HAF. The rest of paper is organized as follows. In Section 2, we present the signal model and overview some of the most popular PPS estimators. The proposed method is presented in Section 3. In this section, the performance study is given as well. Numerical examples supporting the theoretical analysis are given in Section 4 followed by the conclusions in Section 5.
2. CPF and related estimators Consider the following signal model: P
!
xðnÞ ¼ sðnÞ þ νðnÞ ¼ A exp j ∑ ai ðnΔÞi þ νðnÞ; i¼0
n A ½ N=2; N=2;
ð1Þ where νðnÞ is complex zero-mean white Gaussian noise with variance s2 and s(n) is a P-th order PPS with the amplitude A and phase parameters ai ; i ¼ 0; …; P. The number of signal samples is N þ1 and Δ is the sampling rate. Here, N is even positive integer and Δ satisfies the
100
M. Simeunović, I. Djurović / Signal Processing 101 (2014) 99–103
Nyquist–Shannon sampling theorem. Our interest is to estimate fA; a0 ; …; aP g from x(n). The maximum likelihood (ML) estimation of parameters fA; a0 ; …; aP g of higher order PPSs is computationally complex since it requires maximization of a P-dimensional function [8]. Therefore, the phase differentiation (PD)-based techniques are used to reduce the search space [4,5,9]. The PD is performed recursively by the auto-correlation function
The main problem in the CPF evaluation is transforming the underlying signal by quadratic phase matching sequence ∑m ½expð jΩðmΔÞ2 Þ that cannot be evaluated using the FFT. Therefore, the complexity of the CPF-based estimation procedure performed over a search grid of size O(N) is OðN 2 Þ, whereas the HAF-based procedure evaluated using the FFT requires OðN log NÞ operations. This is the reason for proposing the non-uniform sampled CPF in the next section.
1 PDKxðnÞ ½n; τ1 ; τ2 ; …; τK ¼ PDKxðnÞ ½n þτK ; τ1 ; τ2 ; …; τK 1 n
1 fPDKxðnÞ ½n τK ; τ1 ; τ2 ; …; τK 1 g ;
PD0xðnÞ ½n ¼ xðnÞ;
ð2Þ
where K is the number of PDs, τ1 ; τ2 ; …; τK are the lag parameters and PDKxðnÞ ½n; τ1 ; τ2 ; …; τK is the PD operator applied on x(n). The highest order phase parameter aP can be estimated from the HAF by performing ðP 1Þ PDs and periodogram maximization: a^ P ¼
arg maxjHAFðωÞj ω
2P 1 P!ΔP 1 ∏Pi ¼ 11 τi
;
1 ½n; τ1 ; …; τP 1 expð jωΔnÞ: HAFðωÞ ¼ ∑PDPxðnÞ
ð3Þ
n
Once aP is estimated, lower order parameters and amplitude can be obtained from the dechirped signal xd ðnÞ ¼ xðnÞexpð ja^ P ðnΔÞP Þ by repeating the procedure [4]. The PD-based estimation becomes less accurate as the PPS order increases. Each PD generates additional interference terms caused by noise and additional cross-terms in case of several signal components. Moreover, the dechirping procedure causes the error propagation from higher to lower order phase parameters. This effect becomes more emphasized with larger P. The crossterms could be reduced using the product form of the HAF (PHAF) [5], but the main problems associated with the HAF still remain. Therefore, by reducing the number of PDs, issues associated with the PD implementation are also reduced. However, lower number of PDs increases the computational complexity of estimation due to increased dimensionality of the search space. The CPF is introduced for the estimation of cubic phase signals (P ¼3) as [1] CPFðn; ΩÞ ¼ ∑xðn þ mÞxðn mÞexpð jΩðmΔÞ2 Þ:
ð4Þ
m
In the absence of noise, the CPF peaks at the second order phase derivative ΩðnÞ ¼ 2ð3a3 nΔ þ a2 Þ and two highest order phase parameters can be estimated from (4) calculated at n ¼0 and n ¼ n1 . The CPF requires lower number of the PDs with respect to the HAF for cubic phase signals and is correspondingly more accurate than the HAF. The CPF is extended for the estimation of higher order PPSs as [3] CPFHAFðn; ΩÞ ¼ ∑PDP 3 ½n þ m; τ1 ; …; τP 3 m
PDP 3 ½n m; τ1 ; …; τP 3 expð jΩðmΔÞ2 Þ:
3. Non-uniform sampled CPF The parameter estimation of higher order PPSs using a non-uniform sampling scheme has been recently proposed in [6]. Here, we apply the same scheme in the CPF definition in order to enable its evaluation by the FFT algorithms. Maximization of the CPF function is usually performed in two steps. In the first step, a coarse estimate is obtained and that estimate is refined in the second step by searching over the predefined grid around the coarse estimate. The second step can be very computationally demanding. Therefore, several refine search strategies have been proposed for that purpose [10,11]. These strategies are able to obtain very precise estimate by calculating the objective function at several points only. However, these strategies require the evaluation of the objective function by the FFT. Therefore, the CPF evaluation using the FFT would additionally reduce the complexity of the CPF-based estimation procedure. In this paper, we propose pffiffiffiffiffiffi to modify the CPF (4) by substituting m with m ¼ Ck. In that case, the autocorrelation contained in (4) has the following form: pffiffiffiffiffiffi pffiffiffiffiffiffi x1 ðkÞ ¼ xðn þ CkÞxðn CkÞ ¼ A2 expfjð6a3 nΔ3 Ck þ2a2 Δ2 Ckþ 2a3 ðnΔÞ3 þ2a2 ðnΔÞ2 þ a1 nΔ þ 2a0 Þg þνx ðnÞ
ð6Þ
and the corresponding CPF can be written as NUCPFx ðn; ΩÞ ¼ ∑x1 ðkÞexpð jΩCΔ2 kÞ ¼ FTfx1 ðkÞg:
ð7Þ
k
Relationship (7) gives us the CPF representation by means of the FT. Function (7) will be referred to as the nonuniform sampled CPF (NU-CPF). Note that the NU-CPF peaks at ΩðnÞ ¼ 2ð3a3 nΔ þ a2 Þ, same as the CPF. In order to maintain the same nature of noise after resampling and to achieve high resampling factor, C has to be chosen as C N=2 jnj [6]. The NU-CPF requires the evaluation p offfiffiffiffiffi the atffi ffi signal x(n) pffiffiffiffiffi non-integer time instants, i.e., n þ Ck and n Ck. Unknown signal values can be obtained using the interpolation procedure which can be described by two steps [6, Section III C.]:
ð5Þ
This approach is referred as the CPF–HAF, and it uses ðP 3Þ PDs to transform a P-th order PPS to a cubic phase signal, whose parameters are in turn estimated by the CPF.
(a) interpolate x(n) by a factor 2 or 4 using some standard interpolation technique to obtain xi(n); (b) calculate unknown signal value at arbitrary time
M. Simeunović, I. Djurović / Signal Processing 101 (2014) 99–103
instant n0 using the following formula: x i ðn0 Þ ¼
n0 ⌊n0 c ⌈n0 ⌉ n0 x ð⌈n0 ⌉Þ þ x ð⌊n0 cÞ; ⌈n0 ⌉ ⌊n0 c i ⌈n0 ⌉ ⌊n0 c i
ð8Þ
where ⌈n0 ⌉ and ⌊n0 c are two sample times closest to n0. ⌈n0 ⌉ is the upper time value, while ⌊n0 c is the lower one.
Assuming that the interpolation in step (a) is performed by zero padding in the frequency domain, the evaluation of xi(n) for all values of n requires OðNlog NÞ operations. On the other hand, the number of evaluations of (8) in step (b) depends on n used in the NU-CPF and it can be shown that it equals O(N). Therefore, the overall complexity of the interpolation procedure is OðN log NÞ operations. Taking this into account and knowing that maximization of (7) can be performed by the FFT algorithm, the NU-CPF-based estimation procedure requires O(N log N) operations. Similar procedure can be applied for all generalizations of the CPF to higher order PPSs. Here, we consider the CPF–HAF approach. If we, instead of PDP 3 ½n 7 m; pffiffiffiffiffiffi τ1 ; …; τP 3 , evaluate PDP 3 ½n 7 Ck; τ1 ; …; τP 3 in (5), the corresponding non-uniform sampled CPF–HAF (NU-CPF–HAF) has the following form: h i pffiffiffiffiffiffi NUCPFHAFx ðn; ΩÞ ¼ ∑PDP 3 n þ Ck; τ1 ; …; τP 3 k
h i pffiffiffiffiffiffi PDP 3 n Ck; τ1 ; …; τP 3 expð jΩCΔ2 kÞ: The remaining procedure is the same as for (5). Statistical properties of the NU-CPF have been determined in Appendix by the first-order perturbation method [1,3]. The analysis have shown that the NU-CPF is asymptotically unbiased estimator, i.e. Efδa^ p g 0; p ¼ 2; 3; with
MSE [dB]
In this section, we evaluate the performance of the proposed method on PPSs of the third and fifth orders ! Pk
xk ðnÞ ¼ exp j ∑ ai ðnΔÞi þ νðnÞ;
k ¼ 1; 2;
where Δ ¼ 1, N ¼514, P 1 ¼ 3, P 2 ¼ 5 and ai is the i-th element of vector G ¼ f0:9; 0:44; 2:15 10 4 ; 6:42 6 ; 4:12 10 8 ; 1:51 10 10 g. The CPF and NU-CPF are used in the parameter estimation of x1 ðnÞ, while the CPF–HAF and NU-CPF–HAF are used to estimate the parameters of x2 ðnÞ. The HAF estimator is also considered. The HAF-, CPFand CPF–HAF-based estimations are performed following the instructions from [1,3,4]. Monte-Carlo simulations with 300 trials are used. MSEs of the two highest order phase parameters are shown in Fig. 1 where they are compared with the corresponding Cramér–Rao lower bounds (CRLBs). Fig. 1 (a) and (b) corresponds to the parameters of signal x2 ðnÞ, whereas Fig. 1(c) and 1(d) corresponds to the parameters
PPS P=5, estimation of a
4
NU-CPF-HAF CPF-HAF HAF CRLB
-160 -180 -200
-260
SNR [dB] 0
5
MSE [dB]
10
15
20
25
30
5
10
15
20
25
30
PPS P=3, estimation of a2 NU-CPF CPF HAF CRLB
-60 -80
-140
-100
-160 -180 -10
0
MSE [dB] NU-CPF CPF HAF CRLB
-120
SNR [dB]
-220
PPS P=3, estimation of a3
-100
SNR [dB] -5
0
5
10
15
ð10Þ
i¼0
-140
-240
-80
4. Numerical example
MSE [dB]
NU-CPF-HAF CPF-HAF HAF CRLB
1 826:61 2 þ SNR Efðδa^ 3 Þ2 g : N 7 Δ6 SNR ð9Þ
In (9), the MSE E fðδa^ 3 Þ2 g is calculated for n1 ⌊0:0955Nc that minimizes the MSE. By comparing the obtained results with that of the CPF [1], it can be concluded that, for high SNR values, the NU-CPF and the CPF have approximately the same accuracy for a2, whereas the MSE of the NU-CPF for a3 is lower than that of the CPF by about 1 dB. To conclude, with non-uniform sampling we accomplish two aims: reduced calculation complexity and improved accuracy.
5
-220
-280
the MSEs given by 1 48 2 þ SNR Efðδa^ 2 Þ2 g ; Δ4 N5 SNR
PPS P=5, estimation of a
-200
101
SNR [dB]
-120 -10
-5
0
5
10
15
Fig. 1. MSEs of the two highest order phase parameters estimated by the CPF, NU-CPF, CPF–HAF, NU-CPF–HAF and HAF.
102
M. Simeunović, I. Djurović / Signal Processing 101 (2014) 99–103
of signal x1 ðnÞ. As it can be seen, the NU-CPF–HAF outperforms the CPF–HAF and HAF in terms of both the SNR threshold and the MSE. The SNR threshold of the NU-CPF–HAF is around 8 dB, that is by 1 dB and 8 dB lower with respect to the CPF–HAF and HAF, respectively. The MSEs of the NU-CPF–HAF are by about 3 dB (parameter a5) and 1 dB (parameter a4), and by about 5 dB (parameter a5 ) and 3 dB (parameter a4) lower with respect to the MSEs of CPF–HAF and HAF based estimates, respectively. For P¼3, the NU-CPF and CPF have the same SNR thresholds, around 3 dB, whereas the SNR threshold of the HAF is around 3 dB. The MSE is the smallest for the NU-CPF, while it is the largest for the HAF. Note that, above the SNR threshold, the MSE of the NU-CPF approaches the CRLB.
Due to noise, peak of (7) is dislocated from the true value Ω0 ¼ 2ð3a3 nΔ þa2 Þ for [3] 2 1 ∂δf N ðΩ0 Þ ^ ∂ f N ðΩ0 Þ ^ ¼ Ω0 Ω ; δΩ 2 ∂Ω ∂Ω
ð11Þ
where ∂2 g nN ðΩ0 Þ ∂g N ðΩ0 Þ ∂g nN ðΩ0 Þ ∂2 f N ðΩ0 Þ ¼ 2Re g N ðΩ0 Þ þ ; 2 2 ∂Ω ∂Ω ∂Ω ∂Ω ð12Þ ∂δg nN ðΩ0 Þ ∂g N ðΩ0 Þ n ∂δf N ðΩ0 Þ 2Re g N ðΩ0 Þ δg N ðΩ0 Þ : þ ∂Ω ∂Ω ∂Ω ð13Þ Using the following intermediate results:
5. Conclusion
g N ðΩ0 Þ N1 A2 expðjβ0 Þ;
In order to evaluate the CPF by the means of the FFT, in this paper we propose the non-uniformly signal sampling approach. By implementing the FFT algorithm in the CPF, three goals have been accomplished: the calculation complexity of the CPF is reduced, the accuracy is improved and the application of efficient fine estimation technique [10,11] is possible. The performance of the proposed algorithm is evaluated statistically and confirmed by numerical examples.
∂g N ðΩ0 Þ N2 C jA2 exp jβ0 Δ2 C 1 ; ∂Ω 2
∂2 g N ðΩ0 Þ N3 C A2 exp jβ0 Δ4 C 2 1 ; 2 3 ∂Ω N1
δg N ðΩ0 Þ ¼ A2 expðjβ0 Þ ∑ ηðmÞ; m¼0
N1
∂δg N ðΩ0 Þ ¼ jA2 exp jβ0 Δ2 C ∑ mηðmÞ; ∂Ω m¼0
Acknowledgment This work has been supported in part by the Ministry of Science of Montenegro and Foremont project. Appendix A This appendix provides the statistical analysis of the NU-CPF method. The derivation uses the first-order perturbation analysis, so the results hold asymptotically as the number of signal samples approaches infinity. To that end, we use the notation X Y to denote the equality of the dominant (first-order) terms of X and Y. Furthermore, to simplify the MSE derivation, we assume that SNR b 0 dB, i.e., A2 bs2 [3]. In the absence of noise, the NU-CPF is of the following form: N1
g N ðΩÞ ¼ NUCPFx ðn; ΩÞ ¼ A2 expðjβ0 Þ ∑ expðjðΩ0 ΩÞΔ2 CmÞ; m¼0
3
2
where β0 ¼ 2a3 ðnΔÞ þ 2a2 ðnΔÞ þ 2a1 nΔ; Ω0 ¼ 6a3 nΔ þ 2a2 and N 1 ¼ N2 =ð4CÞ. The noise νðnÞ causes perturbations in the NU-CPF that can be expressed by δg N ðΩÞ ¼ NUCPFx ðn; ΩÞ NUCPFs ðn; ΩÞ N1
¼ A2 expðjβ0 Þ ∑ ηðmÞexpðjðΩ0 ΩÞΔ2 CmÞ; m¼0
where pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi νðn þ mC Þ νðn mC Þ νðn þ mC Þνðn mC Þ pffiffiffiffiffiffiffiffi þ pffiffiffiffiffiffiffiffi þ pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi : ηðmÞ ¼ sðn þ mC Þ sðn mC Þ sðn þ mC Þsðn mC Þ
where X C Y means X is approximately equal to Y, we obtain ∂2 f N ðΩ0 Þ 1 C A4 Δ4 C 2 N 41 ; 6 ∂Ω2
ð14Þ
( ) N1 ∂δf N ðΩ0 Þ N1 n 4 2 2A N1 Δ C Im ∑ m η ð mÞ : ∂Ω 2 m¼0
ð15Þ
^ becomes Substituting (14) and (15) into (11), δΩ ( ) N1 N1 ηn ðmÞ 12Im ∑ m 2 m¼0 ^ : δΩ Δ2 CN 31
ð16Þ
From (16), we can conclude that the proposed estimator is ^ ¼ 0. asymptotically unbiased, i.e., EfδΩg Using the identity EfImfXgImfYgg ¼ 12 RefEfXY n g ^ can be expressed as EfXYgg; the MSE of δΩ ( n o N1 N1 72 N1 ^ 2 Re ∑ ∑ p E ðδΩÞ 2 C 2 Δ4 N61 p¼0q¼0 N1 n q E η ðpÞηðqÞ E ηn ðpÞηn ðqÞ : 2 ð17Þ Knowing that the noise is zero-mean white Gaussian, Eq. (17) reduces to 1 18 2 þ N1 SNR ^ 2g þ 1 : EfðδΩÞ 3 C 2 Δ4 N 41 SNR
ð18Þ
M. Simeunović, I. Djurović / Signal Processing 101 (2014) 99–103
Since a2 ¼ Ωð0Þ=2, from (18) it follows: 1 1 96C 2 þ 9 2þ N1 SNR SNR þ1 ¼ Efðδa^ 2 Þ2 g : Δ4 N 6 SNR 2C 2 Δ4 N41 SNR 3
References
ð19Þ
The parameter a3 can be obtained from the estimates of ΩðnÞ and Ωð0Þ as a3 ¼ ðΩðn1 Þ Ωð0ÞÞ=ð6n1 ΔÞ. Thus, we have Efðδa^ 3 Þ2 g ¼
2 ^ ^ 1 ÞÞ2 g 2EfδΩðn ^ 1 ÞδΩð0Þg ^ EfðδΩð0ÞÞ g þ EfðδΩðn : 2 2 36n1 Δ
In order to calculate Efðδa^ 3 Þ2 g; we have ^ 1 ÞδΩð0Þg. ^ EfδΩðn It can be shown that " # 18s2 ð2 þ s2 =A2 Þ N22 ^ ^ þ N1 ; EfδΩ ðn1 ÞδΩ ð0Þg 2 A C 1 C 2 Δ4 N31 N22 3
to
103
find
where N2 ¼ ðN=2 jn1 jÞ2 =C 2 . Here, parameters C2 and N2 correspond to the NU-CPF calculated at n ¼ n1 , while C1 and N1 correspond to the NU-CPF calculated at n ¼0. The final expression for Efðδa^ 3 Þ2 g then reads 1 " # 2þ N1 þ 3 N2 þ 3 N 22 þ3N 1 SNR 2 ^ Efðδa 3 Þ g ¼ þ : 3n21 Δ6 SNR 2N 41 C 21 2N42 C 22 N 31 N 22 C 1 C 2 ð20Þ The optimal n1 can be calculated by minimizing (20) with respect to n1. Assuming that C 1 ¼ N=2 and C 2 ¼ N=2 jn1 j; the minimization yields the following optimal value nopt 1 ⌊0:0955Nc. Substituting C 1 ¼ N=2, C 2 ¼ N=2 jn1 j and nopt into (19) and (20), we obtain the 1 final expressions for Efðδa^ 2 Þ2 g and Efðδa^ 3 Þ2 g given by (9).
[1] P. O'Shea, A fast algorithm for estimating the parameters of a quadratic FM signal, IEEE Trans. Signal Process. 52 (2) (2004) 385–393. [2] P. O'Shea, R. Wiltshire, A new class of multilinear functions for polynomial phase signal analysis, IEEE Trans. Signal Process. 57 (6) (2009) 2096–2109. [3] I. Djurović, M. Simeunović, S. Djukanović, P. Wang, A hybrid CPF–HAF estimation of polynomial-phase signals: detailed statistical analysis, IEEE Trans. Signal Process. 60 (10) (2012) 5010–5023. [4] S. Peleg, B. Porat, Linear FM signal parameter estimation from discrete-time observations, IEEE Trans. Aerosp. Electron. Syst. 27 (4) (1991) 607–616. [5] S. Barbarossa, A. Scaglione, G. Giannakis, Product high-order ambiguity function for multicomponent polynomial-phase signal modeling, IEEE Trans. Signal Process. 46 (3) (1998) 691–708. [6] P. O'Shea, Improving polynomial phase parameter estimation by using nonuniformly spaced signal sample methods, IEEE Trans. Signal Process. 60 (7) (2012) 3405–3414. [7] J.A. Legg, D.A. Gray, Performance bounds for polynomial phase parameter estimation with nonuniform and random sampling schemes, IEEE Trans. Signal Process. 48 (2) (2000) 331–337. [8] I. Djurović, M. Simeunović, B. Lutovac, Are genetic algorithms useful for the parameter estimation of FM signals? Digit. Signal Process. 22 (6) (2012) 1137–1144. [9] P. Wang, I. Djurović, J. Yang, Generalized high-order phase function for parameter estimation of polynomial phase signal, IEEE Trans. Signal Process. 56 (7) (2008) 3023–3028. [10] E. Aboutanios, B. Mulgrew, Iterative frequency estimation by interpolation on Fourier coefficients, IEEE Trans. Signal Process. 53 (4) (2005) 1237–1242. [11] I. Djurović, S. Djukanović, V.V. Lukin, An algorithm for the fine estimation of polynomial-phase signals, IEEE Trans. Aerosp. Electron. Syst. 48 (4) (2012) 3687–3693.
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Fast communication
Non-uniform sampled cubic phase function Marko Simeunović n, Igor Djurović University of Montenegro, Electrical Engineering Department, Cetinjski put bb, 81 000 Podgorica, Montenegro
a r t i c l e i n f o
abstract
Article history: Received 14 December 2013 Received in revised form 24 January 2014 Accepted 5 February 2014 Available online 13 February 2014
Parameter estimation of polynomial phase signals (PPSs) based on the cubic phase function (CPF) and its extensions cannot be performed by using the fast Fourier transform (FT) algorithm. Therefore, in order to express the CPF by means of the FT, in this paper we propose a scheme for the CPF evaluation based on non-uniform sampling. Calculation complexity of the estimation procedure is significantly reduced, whereas the accuracy is the same or better compared to the original algorithm. & 2014 Elsevier B.V. All rights reserved.
Keywords: Polynomial-phase signal Cubic phase function Fourier transform
1. Introduction The cubic phase function (CPF) has been proposed for the estimation of third order polynomial phase signals (PPSs) [1]. Several generalizations of the CPF to PPSs of higher orders exist in the literature, including the recently proposed hybrid CPF–high-order ambiguity function (CPF–HAF) [2,3]. However, since the CPF-based techniques cannot be evaluated using the fast Fourier transform (FFT), they are characterized by larger computational complexity than the FT-based techniques such as the HAF and product HAF (PHAF) [4,5]. With the requirement of OðN2 Þ complex multiplications and additions, where N is the number of signal samples, the CPF can be unfeasible in real-time applications, for example, in radar and sonar applications [3]. A procedure for parameter estimation of higher order PPSs based on the non-uniform signal sampling has been proposed in [6,7]. By sampling the signal at nonequidistant time instants, the procedure lowers the nonlinearity of the estimator function, which, in turn, improves its performance, signal to noise ratio (SNR) threshold and mean squared error (MSE). In this paper, we use a similar non-uniform sampling scheme to
n
Corresponding author. E-mail address: [email protected] (M. Simeunović).
http://dx.doi.org/10.1016/j.sigpro.2014.02.005 0165-1684 & 2014 Elsevier B.V. All rights reserved.
transform the CPF in a form suitable for the FFT algorithm. The modified CPF requires OðN log NÞ operations and has lower MSE than the standard CPF. It can be generalized for other CPF-based techniques. Here, we implement the modified CPF in the evaluation of the CPF–HAF. The rest of paper is organized as follows. In Section 2, we present the signal model and overview some of the most popular PPS estimators. The proposed method is presented in Section 3. In this section, the performance study is given as well. Numerical examples supporting the theoretical analysis are given in Section 4 followed by the conclusions in Section 5.
2. CPF and related estimators Consider the following signal model: P
!
xðnÞ ¼ sðnÞ þ νðnÞ ¼ A exp j ∑ ai ðnΔÞi þ νðnÞ; i¼0
n A ½ N=2; N=2;
ð1Þ where νðnÞ is complex zero-mean white Gaussian noise with variance s2 and s(n) is a P-th order PPS with the amplitude A and phase parameters ai ; i ¼ 0; …; P. The number of signal samples is N þ1 and Δ is the sampling rate. Here, N is even positive integer and Δ satisfies the
100
M. Simeunović, I. Djurović / Signal Processing 101 (2014) 99–103
Nyquist–Shannon sampling theorem. Our interest is to estimate fA; a0 ; …; aP g from x(n). The maximum likelihood (ML) estimation of parameters fA; a0 ; …; aP g of higher order PPSs is computationally complex since it requires maximization of a P-dimensional function [8]. Therefore, the phase differentiation (PD)-based techniques are used to reduce the search space [4,5,9]. The PD is performed recursively by the auto-correlation function
The main problem in the CPF evaluation is transforming the underlying signal by quadratic phase matching sequence ∑m ½expð jΩðmΔÞ2 Þ that cannot be evaluated using the FFT. Therefore, the complexity of the CPF-based estimation procedure performed over a search grid of size O(N) is OðN 2 Þ, whereas the HAF-based procedure evaluated using the FFT requires OðN log NÞ operations. This is the reason for proposing the non-uniform sampled CPF in the next section.
1 PDKxðnÞ ½n; τ1 ; τ2 ; …; τK ¼ PDKxðnÞ ½n þτK ; τ1 ; τ2 ; …; τK 1 n
1 fPDKxðnÞ ½n τK ; τ1 ; τ2 ; …; τK 1 g ;
PD0xðnÞ ½n ¼ xðnÞ;
ð2Þ
where K is the number of PDs, τ1 ; τ2 ; …; τK are the lag parameters and PDKxðnÞ ½n; τ1 ; τ2 ; …; τK is the PD operator applied on x(n). The highest order phase parameter aP can be estimated from the HAF by performing ðP 1Þ PDs and periodogram maximization: a^ P ¼
arg maxjHAFðωÞj ω
2P 1 P!ΔP 1 ∏Pi ¼ 11 τi
;
1 ½n; τ1 ; …; τP 1 expð jωΔnÞ: HAFðωÞ ¼ ∑PDPxðnÞ
ð3Þ
n
Once aP is estimated, lower order parameters and amplitude can be obtained from the dechirped signal xd ðnÞ ¼ xðnÞexpð ja^ P ðnΔÞP Þ by repeating the procedure [4]. The PD-based estimation becomes less accurate as the PPS order increases. Each PD generates additional interference terms caused by noise and additional cross-terms in case of several signal components. Moreover, the dechirping procedure causes the error propagation from higher to lower order phase parameters. This effect becomes more emphasized with larger P. The crossterms could be reduced using the product form of the HAF (PHAF) [5], but the main problems associated with the HAF still remain. Therefore, by reducing the number of PDs, issues associated with the PD implementation are also reduced. However, lower number of PDs increases the computational complexity of estimation due to increased dimensionality of the search space. The CPF is introduced for the estimation of cubic phase signals (P ¼3) as [1] CPFðn; ΩÞ ¼ ∑xðn þ mÞxðn mÞexpð jΩðmΔÞ2 Þ:
ð4Þ
m
In the absence of noise, the CPF peaks at the second order phase derivative ΩðnÞ ¼ 2ð3a3 nΔ þ a2 Þ and two highest order phase parameters can be estimated from (4) calculated at n ¼0 and n ¼ n1 . The CPF requires lower number of the PDs with respect to the HAF for cubic phase signals and is correspondingly more accurate than the HAF. The CPF is extended for the estimation of higher order PPSs as [3] CPFHAFðn; ΩÞ ¼ ∑PDP 3 ½n þ m; τ1 ; …; τP 3 m
PDP 3 ½n m; τ1 ; …; τP 3 expð jΩðmΔÞ2 Þ:
3. Non-uniform sampled CPF The parameter estimation of higher order PPSs using a non-uniform sampling scheme has been recently proposed in [6]. Here, we apply the same scheme in the CPF definition in order to enable its evaluation by the FFT algorithms. Maximization of the CPF function is usually performed in two steps. In the first step, a coarse estimate is obtained and that estimate is refined in the second step by searching over the predefined grid around the coarse estimate. The second step can be very computationally demanding. Therefore, several refine search strategies have been proposed for that purpose [10,11]. These strategies are able to obtain very precise estimate by calculating the objective function at several points only. However, these strategies require the evaluation of the objective function by the FFT. Therefore, the CPF evaluation using the FFT would additionally reduce the complexity of the CPF-based estimation procedure. In this paper, we propose pffiffiffiffiffiffi to modify the CPF (4) by substituting m with m ¼ Ck. In that case, the autocorrelation contained in (4) has the following form: pffiffiffiffiffiffi pffiffiffiffiffiffi x1 ðkÞ ¼ xðn þ CkÞxðn CkÞ ¼ A2 expfjð6a3 nΔ3 Ck þ2a2 Δ2 Ckþ 2a3 ðnΔÞ3 þ2a2 ðnΔÞ2 þ a1 nΔ þ 2a0 Þg þνx ðnÞ
ð6Þ
and the corresponding CPF can be written as NUCPFx ðn; ΩÞ ¼ ∑x1 ðkÞexpð jΩCΔ2 kÞ ¼ FTfx1 ðkÞg:
ð7Þ
k
Relationship (7) gives us the CPF representation by means of the FT. Function (7) will be referred to as the nonuniform sampled CPF (NU-CPF). Note that the NU-CPF peaks at ΩðnÞ ¼ 2ð3a3 nΔ þ a2 Þ, same as the CPF. In order to maintain the same nature of noise after resampling and to achieve high resampling factor, C has to be chosen as C N=2 jnj [6]. The NU-CPF requires the evaluation p offfiffiffiffiffi the atffi ffi signal x(n) pffiffiffiffiffi non-integer time instants, i.e., n þ Ck and n Ck. Unknown signal values can be obtained using the interpolation procedure which can be described by two steps [6, Section III C.]:
ð5Þ
This approach is referred as the CPF–HAF, and it uses ðP 3Þ PDs to transform a P-th order PPS to a cubic phase signal, whose parameters are in turn estimated by the CPF.
(a) interpolate x(n) by a factor 2 or 4 using some standard interpolation technique to obtain xi(n); (b) calculate unknown signal value at arbitrary time
M. Simeunović, I. Djurović / Signal Processing 101 (2014) 99–103
instant n0 using the following formula: x i ðn0 Þ ¼
n0 ⌊n0 c ⌈n0 ⌉ n0 x ð⌈n0 ⌉Þ þ x ð⌊n0 cÞ; ⌈n0 ⌉ ⌊n0 c i ⌈n0 ⌉ ⌊n0 c i
ð8Þ
where ⌈n0 ⌉ and ⌊n0 c are two sample times closest to n0. ⌈n0 ⌉ is the upper time value, while ⌊n0 c is the lower one.
Assuming that the interpolation in step (a) is performed by zero padding in the frequency domain, the evaluation of xi(n) for all values of n requires OðNlog NÞ operations. On the other hand, the number of evaluations of (8) in step (b) depends on n used in the NU-CPF and it can be shown that it equals O(N). Therefore, the overall complexity of the interpolation procedure is OðN log NÞ operations. Taking this into account and knowing that maximization of (7) can be performed by the FFT algorithm, the NU-CPF-based estimation procedure requires O(N log N) operations. Similar procedure can be applied for all generalizations of the CPF to higher order PPSs. Here, we consider the CPF–HAF approach. If we, instead of PDP 3 ½n 7 m; pffiffiffiffiffiffi τ1 ; …; τP 3 , evaluate PDP 3 ½n 7 Ck; τ1 ; …; τP 3 in (5), the corresponding non-uniform sampled CPF–HAF (NU-CPF–HAF) has the following form: h i pffiffiffiffiffiffi NUCPFHAFx ðn; ΩÞ ¼ ∑PDP 3 n þ Ck; τ1 ; …; τP 3 k
h i pffiffiffiffiffiffi PDP 3 n Ck; τ1 ; …; τP 3 expð jΩCΔ2 kÞ: The remaining procedure is the same as for (5). Statistical properties of the NU-CPF have been determined in Appendix by the first-order perturbation method [1,3]. The analysis have shown that the NU-CPF is asymptotically unbiased estimator, i.e. Efδa^ p g 0; p ¼ 2; 3; with
MSE [dB]
In this section, we evaluate the performance of the proposed method on PPSs of the third and fifth orders ! Pk
xk ðnÞ ¼ exp j ∑ ai ðnΔÞi þ νðnÞ;
k ¼ 1; 2;
where Δ ¼ 1, N ¼514, P 1 ¼ 3, P 2 ¼ 5 and ai is the i-th element of vector G ¼ f0:9; 0:44; 2:15 10 4 ; 6:42 6 ; 4:12 10 8 ; 1:51 10 10 g. The CPF and NU-CPF are used in the parameter estimation of x1 ðnÞ, while the CPF–HAF and NU-CPF–HAF are used to estimate the parameters of x2 ðnÞ. The HAF estimator is also considered. The HAF-, CPFand CPF–HAF-based estimations are performed following the instructions from [1,3,4]. Monte-Carlo simulations with 300 trials are used. MSEs of the two highest order phase parameters are shown in Fig. 1 where they are compared with the corresponding Cramér–Rao lower bounds (CRLBs). Fig. 1 (a) and (b) corresponds to the parameters of signal x2 ðnÞ, whereas Fig. 1(c) and 1(d) corresponds to the parameters
PPS P=5, estimation of a
4
NU-CPF-HAF CPF-HAF HAF CRLB
-160 -180 -200
-260
SNR [dB] 0
5
MSE [dB]
10
15
20
25
30
5
10
15
20
25
30
PPS P=3, estimation of a2 NU-CPF CPF HAF CRLB
-60 -80
-140
-100
-160 -180 -10
0
MSE [dB] NU-CPF CPF HAF CRLB
-120
SNR [dB]
-220
PPS P=3, estimation of a3
-100
SNR [dB] -5
0
5
10
15
ð10Þ
i¼0
-140
-240
-80
4. Numerical example
MSE [dB]
NU-CPF-HAF CPF-HAF HAF CRLB
1 826:61 2 þ SNR Efðδa^ 3 Þ2 g : N 7 Δ6 SNR ð9Þ
In (9), the MSE E fðδa^ 3 Þ2 g is calculated for n1 ⌊0:0955Nc that minimizes the MSE. By comparing the obtained results with that of the CPF [1], it can be concluded that, for high SNR values, the NU-CPF and the CPF have approximately the same accuracy for a2, whereas the MSE of the NU-CPF for a3 is lower than that of the CPF by about 1 dB. To conclude, with non-uniform sampling we accomplish two aims: reduced calculation complexity and improved accuracy.
5
-220
-280
the MSEs given by 1 48 2 þ SNR Efðδa^ 2 Þ2 g ; Δ4 N5 SNR
PPS P=5, estimation of a
-200
101
SNR [dB]
-120 -10
-5
0
5
10
15
Fig. 1. MSEs of the two highest order phase parameters estimated by the CPF, NU-CPF, CPF–HAF, NU-CPF–HAF and HAF.
102
M. Simeunović, I. Djurović / Signal Processing 101 (2014) 99–103
of signal x1 ðnÞ. As it can be seen, the NU-CPF–HAF outperforms the CPF–HAF and HAF in terms of both the SNR threshold and the MSE. The SNR threshold of the NU-CPF–HAF is around 8 dB, that is by 1 dB and 8 dB lower with respect to the CPF–HAF and HAF, respectively. The MSEs of the NU-CPF–HAF are by about 3 dB (parameter a5) and 1 dB (parameter a4), and by about 5 dB (parameter a5 ) and 3 dB (parameter a4) lower with respect to the MSEs of CPF–HAF and HAF based estimates, respectively. For P¼3, the NU-CPF and CPF have the same SNR thresholds, around 3 dB, whereas the SNR threshold of the HAF is around 3 dB. The MSE is the smallest for the NU-CPF, while it is the largest for the HAF. Note that, above the SNR threshold, the MSE of the NU-CPF approaches the CRLB.
Due to noise, peak of (7) is dislocated from the true value Ω0 ¼ 2ð3a3 nΔ þa2 Þ for [3] 2 1 ∂δf N ðΩ0 Þ ^ ∂ f N ðΩ0 Þ ^ ¼ Ω0 Ω ; δΩ 2 ∂Ω ∂Ω
ð11Þ
where ∂2 g nN ðΩ0 Þ ∂g N ðΩ0 Þ ∂g nN ðΩ0 Þ ∂2 f N ðΩ0 Þ ¼ 2Re g N ðΩ0 Þ þ ; 2 2 ∂Ω ∂Ω ∂Ω ∂Ω ð12Þ ∂δg nN ðΩ0 Þ ∂g N ðΩ0 Þ n ∂δf N ðΩ0 Þ 2Re g N ðΩ0 Þ δg N ðΩ0 Þ : þ ∂Ω ∂Ω ∂Ω ð13Þ Using the following intermediate results:
5. Conclusion
g N ðΩ0 Þ N1 A2 expðjβ0 Þ;
In order to evaluate the CPF by the means of the FFT, in this paper we propose the non-uniformly signal sampling approach. By implementing the FFT algorithm in the CPF, three goals have been accomplished: the calculation complexity of the CPF is reduced, the accuracy is improved and the application of efficient fine estimation technique [10,11] is possible. The performance of the proposed algorithm is evaluated statistically and confirmed by numerical examples.
∂g N ðΩ0 Þ N2 C jA2 exp jβ0 Δ2 C 1 ; ∂Ω 2
∂2 g N ðΩ0 Þ N3 C A2 exp jβ0 Δ4 C 2 1 ; 2 3 ∂Ω N1
δg N ðΩ0 Þ ¼ A2 expðjβ0 Þ ∑ ηðmÞ; m¼0
N1
∂δg N ðΩ0 Þ ¼ jA2 exp jβ0 Δ2 C ∑ mηðmÞ; ∂Ω m¼0
Acknowledgment This work has been supported in part by the Ministry of Science of Montenegro and Foremont project. Appendix A This appendix provides the statistical analysis of the NU-CPF method. The derivation uses the first-order perturbation analysis, so the results hold asymptotically as the number of signal samples approaches infinity. To that end, we use the notation X Y to denote the equality of the dominant (first-order) terms of X and Y. Furthermore, to simplify the MSE derivation, we assume that SNR b 0 dB, i.e., A2 bs2 [3]. In the absence of noise, the NU-CPF is of the following form: N1
g N ðΩÞ ¼ NUCPFx ðn; ΩÞ ¼ A2 expðjβ0 Þ ∑ expðjðΩ0 ΩÞΔ2 CmÞ; m¼0
3
2
where β0 ¼ 2a3 ðnΔÞ þ 2a2 ðnΔÞ þ 2a1 nΔ; Ω0 ¼ 6a3 nΔ þ 2a2 and N 1 ¼ N2 =ð4CÞ. The noise νðnÞ causes perturbations in the NU-CPF that can be expressed by δg N ðΩÞ ¼ NUCPFx ðn; ΩÞ NUCPFs ðn; ΩÞ N1
¼ A2 expðjβ0 Þ ∑ ηðmÞexpðjðΩ0 ΩÞΔ2 CmÞ; m¼0
where pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi νðn þ mC Þ νðn mC Þ νðn þ mC Þνðn mC Þ pffiffiffiffiffiffiffiffi þ pffiffiffiffiffiffiffiffi þ pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi : ηðmÞ ¼ sðn þ mC Þ sðn mC Þ sðn þ mC Þsðn mC Þ
where X C Y means X is approximately equal to Y, we obtain ∂2 f N ðΩ0 Þ 1 C A4 Δ4 C 2 N 41 ; 6 ∂Ω2
ð14Þ
( ) N1 ∂δf N ðΩ0 Þ N1 n 4 2 2A N1 Δ C Im ∑ m η ð mÞ : ∂Ω 2 m¼0
ð15Þ
^ becomes Substituting (14) and (15) into (11), δΩ ( ) N1 N1 ηn ðmÞ 12Im ∑ m 2 m¼0 ^ : δΩ Δ2 CN 31
ð16Þ
From (16), we can conclude that the proposed estimator is ^ ¼ 0. asymptotically unbiased, i.e., EfδΩg Using the identity EfImfXgImfYgg ¼ 12 RefEfXY n g ^ can be expressed as EfXYgg; the MSE of δΩ ( n o N1 N1 72 N1 ^ 2 Re ∑ ∑ p E ðδΩÞ 2 C 2 Δ4 N61 p¼0q¼0 N1 n q E η ðpÞηðqÞ E ηn ðpÞηn ðqÞ : 2 ð17Þ Knowing that the noise is zero-mean white Gaussian, Eq. (17) reduces to 1 18 2 þ N1 SNR ^ 2g þ 1 : EfðδΩÞ 3 C 2 Δ4 N 41 SNR
ð18Þ
M. Simeunović, I. Djurović / Signal Processing 101 (2014) 99–103
Since a2 ¼ Ωð0Þ=2, from (18) it follows: 1 1 96C 2 þ 9 2þ N1 SNR SNR þ1 ¼ Efðδa^ 2 Þ2 g : Δ4 N 6 SNR 2C 2 Δ4 N41 SNR 3
References
ð19Þ
The parameter a3 can be obtained from the estimates of ΩðnÞ and Ωð0Þ as a3 ¼ ðΩðn1 Þ Ωð0ÞÞ=ð6n1 ΔÞ. Thus, we have Efðδa^ 3 Þ2 g ¼
2 ^ ^ 1 ÞÞ2 g 2EfδΩðn ^ 1 ÞδΩð0Þg ^ EfðδΩð0ÞÞ g þ EfðδΩðn : 2 2 36n1 Δ
In order to calculate Efðδa^ 3 Þ2 g; we have ^ 1 ÞδΩð0Þg. ^ EfδΩðn It can be shown that " # 18s2 ð2 þ s2 =A2 Þ N22 ^ ^ þ N1 ; EfδΩ ðn1 ÞδΩ ð0Þg 2 A C 1 C 2 Δ4 N31 N22 3
to
103
find
where N2 ¼ ðN=2 jn1 jÞ2 =C 2 . Here, parameters C2 and N2 correspond to the NU-CPF calculated at n ¼ n1 , while C1 and N1 correspond to the NU-CPF calculated at n ¼0. The final expression for Efðδa^ 3 Þ2 g then reads 1 " # 2þ N1 þ 3 N2 þ 3 N 22 þ3N 1 SNR 2 ^ Efðδa 3 Þ g ¼ þ : 3n21 Δ6 SNR 2N 41 C 21 2N42 C 22 N 31 N 22 C 1 C 2 ð20Þ The optimal n1 can be calculated by minimizing (20) with respect to n1. Assuming that C 1 ¼ N=2 and C 2 ¼ N=2 jn1 j; the minimization yields the following optimal value nopt 1 ⌊0:0955Nc. Substituting C 1 ¼ N=2, C 2 ¼ N=2 jn1 j and nopt into (19) and (20), we obtain the 1 final expressions for Efðδa^ 2 Þ2 g and Efðδa^ 3 Þ2 g given by (9).
[1] P. O'Shea, A fast algorithm for estimating the parameters of a quadratic FM signal, IEEE Trans. Signal Process. 52 (2) (2004) 385–393. [2] P. O'Shea, R. Wiltshire, A new class of multilinear functions for polynomial phase signal analysis, IEEE Trans. Signal Process. 57 (6) (2009) 2096–2109. [3] I. Djurović, M. Simeunović, S. Djukanović, P. Wang, A hybrid CPF–HAF estimation of polynomial-phase signals: detailed statistical analysis, IEEE Trans. Signal Process. 60 (10) (2012) 5010–5023. [4] S. Peleg, B. Porat, Linear FM signal parameter estimation from discrete-time observations, IEEE Trans. Aerosp. Electron. Syst. 27 (4) (1991) 607–616. [5] S. Barbarossa, A. Scaglione, G. Giannakis, Product high-order ambiguity function for multicomponent polynomial-phase signal modeling, IEEE Trans. Signal Process. 46 (3) (1998) 691–708. [6] P. O'Shea, Improving polynomial phase parameter estimation by using nonuniformly spaced signal sample methods, IEEE Trans. Signal Process. 60 (7) (2012) 3405–3414. [7] J.A. Legg, D.A. Gray, Performance bounds for polynomial phase parameter estimation with nonuniform and random sampling schemes, IEEE Trans. Signal Process. 48 (2) (2000) 331–337. [8] I. Djurović, M. Simeunović, B. Lutovac, Are genetic algorithms useful for the parameter estimation of FM signals? Digit. Signal Process. 22 (6) (2012) 1137–1144. [9] P. Wang, I. Djurović, J. Yang, Generalized high-order phase function for parameter estimation of polynomial phase signal, IEEE Trans. Signal Process. 56 (7) (2008) 3023–3028. [10] E. Aboutanios, B. Mulgrew, Iterative frequency estimation by interpolation on Fourier coefficients, IEEE Trans. Signal Process. 53 (4) (2005) 1237–1242. [11] I. Djurović, S. Djukanović, V.V. Lukin, An algorithm for the fine estimation of polynomial-phase signals, IEEE Trans. Aerosp. Electron. Syst. 48 (4) (2012) 3687–3693.