Multicomponent chirp signals analysis using product cubic phase function

Digital Signal Processing 16 (2006) 654–669 www.elsevier.com/locate/dsp

Multicomponent chirp signals analysis using product cubic phase function ✩ Pu Wang ∗ , Jianyu Yang School of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu 610054, PR China Available online 6 October 2006

Abstract This paper presents an algorithm for estimating the parameters of multicomponent chirp signals. The estimator is based on the cubic phase function (CPF), which is efficient to estimate the parameters of monocomponent first-, second-, and third-order polynomial phase signal. When the CPF is dealing with multicomponent chirp signals, the spurious peaks arise and hence the identifiability problem exists. A new approach based on the product cubic phase function (PCPF) is proposed to remove the identifiability problem. The occurrence probability of spurious peak, and the effect of noise to the estimator are statistically studied. The simulation examples are provided for validating the theoretical analysis. © 2006 Elsevier Inc. All rights reserved. Keywords: Parameter estimation; Chirp signal; Cubic phase function

1. Introduction Linear frequency-modulated (LFM) signals, also known as chirp signals, are often found in signal processing and communication applications such as radar, sonar, biomedicine, seismic analysis, and mobile communications. These signals are nonstationary and usually embedded in noise. There are various approaches proposed for estimating this kind of signal in this literature. The time-frequency distributions, such as the Wigner–Ville distribution (WVD) and its related bilinear class [1,2], are efficient to reveal the instantaneous frequency (IF) over the time-frequency plane. Among the parametric methods, the maximum likelihood estimation (MLE) [4] plays an important role. The discrete form of the MLE, known as discrete chirps Fourier transform (DCFT), is presented in [5]. Due to the exhausting two-dimensional (2-D) maximization of the MLE, the suboptimal methods are preferred to fast implementation. The phase unwrapping method [3] is proposed based on the finite difference operator. It however suffers from the incapability to analyze the multicomponent signals. The polynomial phase transform (PPT) [6,7] with order 2, Radon–Ambiguity transform (RAT) [8] and fractional autocorrelation [9] convert the signal into one-dimensional (1-D) parameter space in which the chirp rate is only interested. Then, the demodulation is applied to estimate other parameters. The fractional Fourier transform (FrFT) [10] has received considerable attention as well for its improved noise rejection and ability to resolve closely-spaced chirp signals. ✩ This paper was in part presented at the 31st International Conference on Acoustics, Speech, and Signal Processing, Toulouse, France, May 14–19, 2006. * Corresponding author. E-mail address: [email protected] (P. Wang).

1051-2004/$ – see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.dsp.2006.09.002

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

655

In this paper, we consider the recently proposed transform, called cubic phase function (CPF) [11,12], which first estimates the instantaneous frequency rate (IFR) and then uses the IFR as initial step to estimate other parameters. When the CPF is applied to multicomponent chirp signals, the spurious peak can be identified and thus the identifiability problem occurs. Theoretical analysis shows that the spurious peak may appear within the observation time 2 with probability 1 − (1/2)(K −K)/2 , where K is the number of components, provided that the phase parameters are independently and uniformly distributed. To solve this problem, the product CPF (PCPF) is presented. By exploring different time dependencies of the auto-terms and cross-terms (including spurious peaks), the PCPF aligns the auto-terms and disaligns the cross-terms, and the later multiplication enhances the auto-terms while weakening the cross-terms. This is motivated by the idea of the product higher-order ambiguity function (PHAF) proposed in [13], which discerns the auto-terms and cross-terms by exploiting different lag dependencies. Herein, we utilize different time dependencies. 2. Brief review of cubic phase function We consider a chirp signal defined as N −1 N −1 ,..., , (1) 2 2 where N is assumed odd and φ(n) is the instantaneous phase of the chirp signal. The CPF for a signal x(n) is defined by x(n) = Aeφ(n) = Aej (a0 +a1 n+a2 n ) , 2

∞ CPF(n, Ω) =

n=−

x(n + τ )x(n − τ )e−j Ωτ dτ, 2

(2)

0

where Ω represents the IFR that is defined in (3) d 2 φ(n) . dn2 Substituting x(n) in (2) with (1) yields IFR(n) =

2 j 2(a0 +a1 n+a2 n2 )

∞

CPF(n, Ω) = A e

(3)

ej (2a2 −Ω)τ dτ. 2

(4)

0

It is easy to obtain that, when Ω = 2a2 , the CPF achieves the maximum along the IFR. Once the IFR is estimated, the a2 can be determined by a2 = Ω/2. The CPF-based estimator for a2 is statistical efficient, since its mean-square error (MSE) approaches the Cramer–Rao lower bound (CRLB) asymptotically at high SNR at n = 0 [11]. When the CPF is applied to multicomponent chirp signals, however, the identifiability problem arises. This claim can be verified by the following section. 3. Problem formulation 3.1. Appearance of spurious peak The multicomponent chirp signals can be modeled in discrete time as x(n) =

K 

Ak ej (ak,0 +ak,1 n+ak,2 n ) , 2

k=1

n=−

(N − 1) (N − 1) ,..., , 2 2

(5)

where K is the number of components, {ak,i }2i=0 is the ith-order phase coefficient for the kth component, and Ak is the amplitude which is assumed to be real and positive constant. To formulate the identifiability problem arising from multicomponent signals, we first consider two chirp signals (K = 2): x(n) = A1 ej (a1,0 +a1,1 n+a1,2 n ) + A2 ej (a2,0 +a2,1 n+a2,2 n ) . 2

2

(6)

656

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

Substituting (6) to (2) yields 2 CPF(n, Ω) = A21 ej 2(a1,0 +a1,1 n+a1,2 n )

∞ e

j (2a1,2 −Ω)τ 2

dτ

2 + A22 ej 2(a2,0 +a2,1 n+a2,2 n )

0

∞ + A1 A2 z(n)

∞

ej (2a2,2 −Ω)τ dτ 2

0

ej (a1,2 +a2,2 −Ω)τ ej {(a1,1 −a2,1 )+2(a1,2 −a2,2 )n}τ dτ 2

0

∞ + A1 A2 z(n)

ej (a1,2 +a2,2 −Ω)τ ej {(a2,1 −a1,1 )+2(a2,2 −a1,2 )n}τ dτ, 2

(7)

0

where z(n) = ej {(a1,0 +a2,0 )+(a1,1 +a2,1 )n+(a1,2 +a2,2 )n } . From (7), the CPF presents two peaks at Ω = 2a1,2 and Ω = 2a2,2 corresponding to the first two terms in the right of (7), whereas two cross-terms disperse along the IFR domain since the positions of cross-terms vary with time. It is interesting to see what happens in (7) if a1,1 = a2,1 . Assuming a1,1 = a2,1 , the last two terms reduce to 2

∞ A1 A2 z(n)

e

j (a1,2 +a2,2 −Ω)τ 2 j 2(a1,2 −a2,2 )nτ

e

∞ dτ,

A1 A2 z(n)

0

ej (a1,2 +a2,2 −Ω)τ ej 2(a2,2 −a1,2 )nτ dτ. 2

(8)

0

In this case, (7) at n = 0 has the form as ∞ CPF(0, Ω) = A21 ej 2a1,0

e

j (2a1,2 −Ω)τ 2

∞ dτ

+ A22 ej 2a2,0

0

ej (2a2,2 −Ω)τ dτ 2

0

+ 2A1 A2 ej (a1,0 +a2,0 )

∞

ej (a1,2 +a2,2 −Ω)τ dτ. 2

(9)

0

It is clear that another peak (spurious peak), which maximizes at Ω = a1,2 + a2,2 , will be present at n = 0. Generally speaking, the cross-terms merge into spurious peaks at the time point that satisfies the following equation: (a2,1 − a1,1 ) + 2(a2,2 − a1,2 )nc = 0,

(10)

where nc represents the time position of the presentation of spurious peak. The CPF slice at n = nc therefore is given by 2 CPF(nc , Ω) = A21 ej 2(a1,0 +a1,1 nc +a1,2 nc )

∞ e

j (2a1,2 −Ω)τ 2

dτ

2 + A22 ej 2(a2,0 +a2,1 nc +a2,2 nc )

0

∞ + 2A1 A2 z(nc )

∞

ej (2a2,2 −Ω)τ dτ 2

0

ej (a1,2 +a2,2 −Ω)τ dτ. 2

(11)

0

As a result, the CPF presents a spurious peak locating at Ω = a1,2 + a2,2 again. The identifiability problem for multicomponent case with K > 2 can be analyzed in similar way of K = 2. Thus, we can conclude that, for K > 2, • The more components exist, the more cross-terms arise. The number of cross-terms is K 2 − K; • The cross-terms merge into (K 2 − K)/2 spurious peaks at n = nl,c , where nl,c satisfy the following equations set: (ak,1 − am,1 ) + 2(ak,2 − am,2 )nl,c = 0,

1  k = m  K;

l = 1, . . . ,

K2 − K . 2

(12)

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

657

3.2. Occurrence probability of spurious peak To identify the occurrence probability of spurious peak over the observation time [−(N − 1)/2:(N − 1)/2] (abbreviate it as [−N/2:N/2]), according to (12), we need to first compute the the probability density functions (PDF) of b1 = [am,1 − ak,1 ] and b2 = 2[ak,2 − am,2 ], and then the PDF of nl,c = b1 /b2 over [−N/2:N/2]. Prior to the analysis, some preliminaries are provided. 1) To avoid ambiguities arising from the cyclic nature of spectral transforms of sampled signals [6,12], it is assumed that |a1 |  π and that |a2 |  π/N ; 2) ak,1 obeys uniform distribution within [−π:π]  1 , −π  ak,1  π, p(ak,1 ) = 2π (13) 0, else, and ak,2 obeys uniform distribution within [−π/N :π/N ] N , −π/N  ak,2  π/N, p(ak,2 ) = 2π 0, else,

(14)

where ak,1 is independent of ak,2 ; 3) ak,i is independently and identically distributed (i.i.d.) with am,i , where 1  k = m  K. Under Preliminary 3), the PDF of b1 and b2 can be obtained, respectively (see Appendix A), ⎧ 2π+b 1 ⎪ ⎨ 4π 2 , −2π  b1 < 0, p(b1 ) = 2π−b1 , 0  b1  2π, ⎪ 4π 2 ⎩ 0, else, and p(b2 ) =

(15)

⎧ 4πN +b2 N 2 ⎪ ⎪ ⎨ 16π 2 , −4π/N  b2 < 0, ⎪ ⎪ ⎩

4πN −b2 N 2 , 16π 2

0,

0  b2  4π/N, else.

(16)

Since ak,1 is independent of ak,2 , b1 is independent of b2 . Based on the independent property, the PDF of nl,c = b1 /b2 located in [−N/2:N/2] is (see Appendix B) ⎧ ) ⎨ 2(nl,c +N , −N/2  nl,c < 0, 3N 2 p(nl,c ) = 2(n −N ) (17) ⎩ l,c , 0  nl,c  N/2. 3N 2 Therefore, the occurrence probability of spurious peak in [−N/2:N/2] can be calculated by integrating up p(nl,c ) over [−N/2:N/2] yields,

N 1 N Pr −  nl,c  = . (18) 2 2 2 It means that two cross-terms between any two components may merge into a spurious peak within [−N/2:N/2] with probability 50%, provided that the phase parameters are independently and uniformly distributed. For K chirp signals, the occurrence probability of spurious peak within [−N/2:N/2] hence is

(K 2 −K)/2 N 1 N Pr −  nl,c  =1− . 2 2 2 For example, K = 4, the occurrence probability of spurious peak is about 98.44%.

(19)

658

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

3.3. Summary Based on the analysis above, the algorithm in [11,12] is not suitable to estimate multicomponent chirp signals due to the existing of cross-terms and its resulting spurious peaks. Furthermore, there is not a principle which appropriately selects the time positions to avoid the identifiability problem. Thus, we need to improve the CPF-based algorithm for multicomponent chirp signals. 4. Product cubic phase function The above section establishes the existing of cross-terms and spurious peaks. We also find that the cross-terms disperse across the IFR domain and turn into spurious peaks when (12) is satisfied, whereas the auto-terms concentrate on several straight lines that locate at Ω = 2ak,2 . Moreover, the auto-terms are independent on time while the crossterms are related to the time positions. Hence, it provides the basis for discerning the auto-terms from the cross-terms, even for the cross-terms give rise to spurious peaks. In this paper, we combine several CPF slices at different time positions by multiplying the corresponding CPF slices, which results in the product CPF (PCPF). 4.1. Definition Given the set of L different time positions nl , where nl ∈ [−(N − 1)/2:(N − 1)/2], the CPF(nl , Ω) corresponding to the set of time are computed by (2) and then the PCPF as the product of the CPFs at different time positions is defined as PCPF(Ω) =

L 

CPF(nl , Ω).

(20)

l=1

From (20), the PCPF will peak at Ω = 2ak,2 . It provides a number of advantages: • The multiply operation amplifies the auto-terms due to the fact that the auto-terms align, while weakening the cross-terms due to its dispersion along the IFR domain; • The spurious peaks are suppressed and even eliminated; • The noise rejection is improved with respect to that of CPF, when L is sufficiently large. 4.2. Estimation algorithm First, assume that the number of components K is known and the amplitudes of the components are equal or little different. The detailed procedure for multicomponent chirp signal is as follows: Step 1. Initialize k = 1, x1 (n) = x(n). Step 2. Estimate aˆ k,2 by searching the highest peak in the PCPF. Step 3. Estimate aˆ k,1 by dechirping and finding the Fourier transform peak: (N −1)/2  2 aˆ k,1 = arg max xk (n)e−j (ak,1 n+aˆ k,2 n ) . ak,1 n=−(N −1)/2

Step 4. Estimate aˆ k,0 and Aˆ k by evaluating −1)/2 1 (N 2) −j ( a ˆ n+ a ˆ n k,1 k,2 xk (n)e Aˆ k = , N n=−(N −1)/2

(N −1)/2   −j (aˆ k,1 n+aˆ k,2 n2 ) aˆ k,0 = angle xk (n)e . n=−(N −1)/2

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

659

Step 5. Peel off the estimated kth strong component by xk+1 = xk − Aˆ k ej (aˆ k,0 +aˆ k,1 n+aˆ k,2 n ) . 2

And set k = k + 1. Step 6. Repeat Steps 2–5 until k = K. For single chirp signal, just simply set K = 1. For multicomponent signals with unbalanced amplitudes, the algorithm extension can be performed similarly using the procedure in [14]. 4.3. Detection of a peak As explained in Section 4.1, the magnitude of PCPF forms a peak corresponding to a2 . However, due to the existing of noise and multicomponent signals, one has to decide whether the maximum can be determined as a peak. This is related to the choice of a suitable criterion and a threshold. In this paper, the decision criterion proposed in [14,15] is employed. In the selected criterion, a ratio P , which is defined as the ratio of the peak and its immediate neighborhood floor, is used and defined as P (τ ) =

maxΩ∈[−π/N,π/N ) |PCPF(Ω)| , medianΩ∈[Ω0 −ζ,Ω0 +ζ ] |PCPF(Ω)|

(21)

where Ω0 is the position of the maximum and ζ describes the neighborhood range around the maximum. Then, we compare it with a threshold γ and declare it is indeed a peak if P  γ . In [14,15], the authors recommended that γ ∈ [4, 5]. Obviously, a higher threshold increases the probability of detection, while decreasing the probability of false alarms. We applied this criterion with γ = 5 and ζ = 7π/N 2 to Examples 2 and 5, and the ratios P are 2.4624 and 25.1672, respectively. Compared the ratio with the threshold γ , it can be said that, under this criterion, the PCPF can easily identify a distinct peak, while the CPF fails to detect the chirp signal. Remark I. If the decision criterion in [14,15] is applied to the noise-free signal in Example 3, the ratios of CPF and PCPF are 2.9925 and 27.4268, respectively. This means that, under this criterion, the CPF is not able to detect any peaks in its spectrum even when the signal observations are noise-free. Remark II. There is a relationship between the ratio γ and signal length N . Larger observations lead to higher peak, but heavier computation burden. The PCPF can thus be treated as an alternative to improve the ratio γ without dramatically increasing computations and with much less observations. 4.4. Statistical analysis of the a2 estimate Since the estimation algorithm is iterative, it inevitably suffers from the error propagation effect. Error in the a2 estimate will propagate to the latter estimates. Hence, the statistical analysis of the a2 estimate is the most critical step and is detailed in Appendix C. For other estimates, the statistical properties can be deduced in a similar way in [16]. Table 1 gives the asymptotic MSE of the a2 estimate using alternative methods such as the HAF-based and CPF-based methods. In Table 1, the PCPF is specified as the one with two time positions, n1 = 0.4N and n2 = 0. Compared with other methods, the PCPF has good trade-off between the noise rejection and ability for analyzing multicomponent chirp signals, i.e., the HAF-based method can adapt to multicomponent chirp signal, but has higher asymptotic MSE than Table 1 Asymptotic MSEs (AMSE) of different methods and CRLB for the a2 estimate AMSE (PCPF)

AMSE (CPF)

AMSE (HAF)

CRLB

90(1.0183+1/1.9905SNR) SNR×N 5

90(1+1/2SNR) SNR×N 5

96(1+1/2SNR) SNR×N 5

90 SNR×N 5

660

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

PCPF, while the CPF approaches the CRLB, but is disable to analyze the multicomponent signals. Another advantage of the PCPF is the lower of the threshold SNR as increasing of the L, i.e., −5 dB of L = 3 and −6 dB of L = 5. Remark III. The statistical analysis also indicts that the the choice of n1 = 0 and n2 = 0 has minimum variance of the a2 estimate, i.e., the asymptotic MSE approaches the CRLB at high SNR, however, this choice leads the PCPF with L = 2 to be the square of CPF, which means that the time selection cannot improve the ability for analyzing multicomponent signals. To deal with multicomponent signals, the selected time positions should be sufficiently spaced to misalign the cross-terms and spurious peaks. Therefore, there is a trade-off between the estimate variance and the ability for multicomponent signals. In this paper, we select n1 = 0.4N and n2 = 0. For the PCPF with L > 2, the time positions should be sufficiently spaced as well. 4.5. Computational cost To compare the PCPF with the CPF in terms of computational cost, we have that, at each iteration, the PCPF requires • The computation of L PCPF’s; • L − 1 vector product. Hence, the PCPF requires L times more computations than the CPF. Fortunately, the computation of CPF can be reduced to the order of N log2 N with the subband decomposition technique [11,12]. Moreover, L is always a small number, which means that the additional cost is not excessive. 5. Numerical analysis In this section, comprehensive numerical examples are provided to illustrate the proposed algorithm. Example 1. Consider two chirp signals with the same A and a1 . The selected parameters are listed below: 1st component: A1 = 1, 2nd component: A2 = 1,

a1,0 = 0, a2,0 = 0,

a1,1 = 2π/5, a2,1 = 2π/5,

a1,2 = π/5N, a2,2 = −2π/5N.

The signal length N = 515. Figure 1 gives the CPF slice at n = 0, |CPF(0, Ω)|. As shown in Fig. 1, two cross-terms merge to a spurious peak at Ω = a1,2 + a2,2 = −π/5N , while the auto-terms locate at Ω = 2a1,2 = 2π/5N and Ω = 2a2,2 = −4π/5N , respectively. Example 2. Consider two chirp signals with the same A but different a1 . The selected parameters are listed below: 1st component: A1 = 1, 2nd component: A2 = 1,

a1,0 = 0, a2,0 = 0,

a1,1 = π/5, a2,1 = 2π/5,

a1,2 = π/5N, a2,2 = −2π/5N.

In this example, the time point that satisfies (12) is n ≈ 86. Then we plot the CPF slice at n = 86, |CPF(86, Ω)|. As shown in Fig. 2, the cross-terms merge to a spurious peak at Ω = a1,2 + a2,2 = −π/5N again. Example 3 (To testify the improvement of noise rejection). Consider single chirp signal that is exactly the first component in Example 1 embedded in −5 dB additive white Gaussian noise. The selected time positions used for the PCPF are [−100, −50, 0, 50, 100]. We plot the |CPF(0, Ω)| and |PCPF(Ω)| in Figs. 3 and 4, respectively. While the peak is overwhelmed by the noise in Fig. 3, the PCPF in Fig. 4 shows a distinct peak at Ω = 2a2 . And the enhancement of SNR from CPF to PCPF is obvious. Example 4 (To demonstrate suppression of cross terms and spurious peaks). Consider the same chirp signals in Example 2.

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

661

Fig. 1. The CPF at n = 0 of two chirp signals with same a1 .

Fig. 2. The CPF at n = 86 of two chirp signals with different a1 .

In this example, the PCPF has been computed by using five different time positions. The time positions are the followings: [−100, −50, 0, 86, 100]. It is worthy noting that the time sets include the position at which spurious peak arises (n = 86). Figure 5 shows the PCPF of the same signals in Fig. 2. Two peaks corresponding to the auto-terms in PCPF are evident, whereas the spurious peak is almost eliminated. Example 5 (To evaluate the performance of the a2 estimate for single chirp signal using PCPF). The generated signal is the same signal in Example 3 and the SNR varies from −8 to 8 dB in steps of 1 dB. For each SNR, 200 independent experiments are realized. The time sets are [0] for L = 1, [0, 0.4N ] for L = 2, [−0.4N , 0, 0.4N ] for L = 3, and [−0.4N , −0.2N , 0, 0.2N , 0.4N ] for L = 5, respectively. The MSEs of the a2 estimate corresponding to different sets of time positions are evaluated and presented in Fig. 6. We can observe the measured results agree with the theoretical MSE for PCPF with L = 2 at high SNR. At low SNR, the threshold effect arises at about SNR = −4 dB due to the nonlinearity.

662

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

Fig. 3. The CPF at n = 0 of single chirp signal in −5 dB noise.

Fig. 4. The PCPF of the same signal as in Fig. 3.

One advantages of the PCPF is the decrease of the SNR threshold as increasing the sets of time positions. It can be observed that in Fig. 7, the PCPF with L = 5 shows a threshold SNR at −6 dB, which is about 2 dB lower than that of the CPF-based estimator. Example 6 (To evaluate the performance of the a2 estimate for multicomponent chirp signals using PCPF). Due to the PCPF is proposed to solve the identifiability problem when the CPF deals with multicomponent chirp signals, we present the performance of the PCPF when multicomponent signals are encountered. Two chirp signals are considered, which are the same signals in Example 1. And the SNR varies from −8 to 8 dB in steps of 1 dB. For each SNR, 200 independent experiments are realized. The MSEs of the a2 estimate versus the SNR are obtained for a number of L sets of time positions such as L = 1, 3, 5. The time sets are the same as that in Example 5. Note that the spurious peak will present at n = 0, which is the choice of time position for the PCPF with L = 1.

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

663

Fig. 5. The PCPF of the same signal as in Fig. 2.

Fig. 6. Comparison between theoretical results, measured results and CRLB for the a2 estimate with L = 2.

From Fig. 8, we can observe that the CPF at n = 0 cannot deal with multicomponent signals as the SNR increases due to the occurrence of spurious peak. However, the PCPF with L > 1 removes the identifiability problem and derives small variance at low SNR. 6. Discussions For multicomponent higher-order FM signals, the CPF presents spurious peaks. The PCPF, however, cannot discern the spurious peaks from the auto-terms, due to the auto-terms are also related to the time positions. Thus, some furtherdeveloped methods are required and will be considered in the future. 7. Conclusion In this paper, the spurious peaks are first identified when the CPF deals with multicomponent chirp signals. By exploring different time dependence of the auto-terms and cross-terms, a PCPF-based algorithm is proposed for re-

664

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

Fig. 7. Measured MSEs for the a2 estimate (monocomponent case).

Fig. 8. Measured MSEs for the a2 estimate (multicomponent case).

moving the identifiability problem arising from the spurious peaks. The first-order permutation analysis is performed for deriving the asymptotic MSE of the a2 estimate. Compared with the CPF, the PCPF has a slight higher asymptotic MSE at high SNR but has a lower threshold SNR and the ability for analyzing multicomponent chirp signals. The computation cost and time selection are also discussed. A number of numerical examples are provided to identify the theoretical results. Appendix A. Proof of (15) and (16) Lemma 1. If random variables X and Y are independent, identically and uniformly distributed over [−a, a], the PDF of Z = X − Y is ⎧ 2a+z ⎪ ⎨ 4a 2 , −2a  z < 0, pZ (z) = 2a−z , 0  z  2π, (A.1) ⎪ ⎩ 4a 2 0, else.

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

665

Fig. 9. The illustration of G.

Proof. To compute the PDF of Z, the following equation exists: ∞ pZ (z) =

pXY (x, x − z) dx,

(A.2)

−∞

where pXY (x, x − z) is the joint PDF for X and Y . Under the independent assumption, the equation above can be rewrote as ∞ pZ (z) =

pX (x)pY (x − z) dx.

(A.3)

−∞

Since X and Y are uniformly distributed, the pXY is  1 , (x, z) ∈ G, pXY (x, x − z) = 4a 2 0, else,

(A.4)

where G = {(x, z): |x|  a, |x − z|  a} shown in Fig. 9. When 0  z  2a, 1 pZ (z) = 2 4a

a−z dx = −a

2a − z . 4a 2

(A.5)

When −2a  z < 0, 1 pZ (z) = 2 4a

a dx = −a−z

When |z| > 2a, pZ (z) = 0.

2a + z . 4a 2

(A.6)

2

Appendix B. Proof of (17) The PDF of Z = X/Y can be calculated by [17] ∞ pZ (z) = −∞

|y|pXY (zy, y) dy.

(B.1)

666

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

Fig. 10. The illustration of G1 , G2 , G3 , and G4 .

Herein, let X and Y denote b1 and b2 , respectively. And the |y|pXY (x, y) can be obtained as ⎧ 4πN −N 2 Y ⎪ y 2π−zy , (y, z) ∈ G1 , ⎪ 2 4π 16π 2 ⎪ ⎪ ⎪ 2 ⎪ 2π+zy 4πN +N Y ⎨ −y , (y, z) ∈ G2 , 4π 2 16π 2 |y|pXY (x, y) = 2π−zy 4πN +N 2 Y ⎪ ⎪ , (y, z) ∈ G3 , −y 4π 2 ⎪ 16π 2 ⎪ ⎪ ⎪ ⎩ 2π+zy 4πN −N 2 Y , (y, z) ∈ G4 , y 4π 2 16π 2

(B.2)

where G1 , G2 , G3 and G4 are shown in Fig. 10. Thus, if 0  z  N/2, 0 pZ (z) = −4π/N

2π + zy 4πN + N 2 Y −y dy + 4π 2 16π 2

4π/N 

y

2π − zy 4πN − N 2 Y 2(z − N ) dy = 2 2 4π 16π 3N 2

(B.3)

2π + zy 4πN − N 2 Y 2(z + N ) dy = . 4π 2 16π 2 3N 2

(B.4)

0

and if −N/2  z < 0, 0 pZ (z) = −4π/N

2π − zy 4πN + N 2 Y −y dy + 4π 2 16π 2

4π/N 

y 0

Substituting z of nl,c , (17) exists. Appendix C. Asymptotic MSE of the a2 estimate In this appendix, we perform the first-order permutation analysis of the a2 estimate using the PCPF with L = 2. Initially, assume that the selected time positions for the PCPF with L = 2 are n1  0 and n2 = 0. Note that the choice of n2 = 0 arises from the asymptotic optimal for the a2 estimate using the CPF at n = 0 [11,12]. For the noise-free chirp signal, the PCPF with L = 2 is gN (Ω) =

N1  m=0

s(n1 + m)s(n1 − m)e

−j Ωm2

N2 

s(n2 + m)s(n2 − m)e−j Ωm , 2

(C.1)

m=0

where N1 = (N − 1)/2 − n1 and N2 = (N − 1)/2 − n2 . The PCPF will peak at a global maximum Ω0 = 2a2 . Then additive noise v(n) has the perturbation, δgN (Ω), on the gN (Ω) and moves the peak to Ω = 2a2 + δΩ. In case of high SNR, the expression of the permutation can be approximated by keeping only the cross-terms containing not more than one noise factor and the noise terms containing two noise factors. Therefore, the perturbation can be approximated as

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

δgN (Ω) =

N1 

s(n1 + m)s(n1 − m)e−j Ωm

m=0

+

2

N2 

xvs (n2 )e−j Ωm

667

2

m=0

N2 

s(n2 + m)s(n2 − m)e−j Ωm

m=0

2

N1 

xvs (n1 )e−j Ωm , 2

(C.2)

m=0

where xvs (n1 ) = s(n1 + m)v(n1 − m) + s(n1 − m)v(n1 + m) + v(n1 − m)v(n1 + m) and xvs (n2 ) = s(n2 + m)v(n2 − m) + s(n2 − m)v(n2 + m) + v(n2 − m)v(n2 + m). From the derivations in [16], the asymptotic MSE of the peak fluctuations is   E{α 2 } E (δΩ)2 ≈ , β2

(C.3)

where E{·} denotes expected value and where  ∗ (Ω ) ∗ (Ω )  ∂ 2 gN ∂gN (Ω0 ) ∂gN 0 0 , α = 2Re gN (Ω0 ) + ∂Ω ∂Ω ∂Ω 2   ∗ (Ω ) ∂δgN ∂gN (Ω0 ) ∗ 0 β = 2Re gN (Ω0 ) + δgN (Ω0 ) . ∂Ω ∂Ω

(C.4) (C.5)

Using the above equations, the intermediate results are deduced: gN (Ω0 ) ≈ A4 K 2 N1 N2 ,  1 ∂gN (Ω0 ) ≈ −j A4 K 2 N1 N23 + N13 N2 , ∂Ω 3  ∂ 2 gN (Ω0 ) 1 9N1 N25 + 10N13 N23 + 9N15 N2 , ≈ −A4 K 2 45 ∂Ω 2 N1 N2   2 2 δgN (Ω0 ) ≈ A2 KN2 xvs (n1 )e−j Ω0 m + A2 KN1 xvs (n2 )e−j Ω0 m , m=0

(C.6) (C.7) (C.8) (C.9)

m=0



N1 N2   N22 N12 ∂δgN (Ω0 ) 2 2 2 −j Ω0 m2 2 2 m + m + ≈ −j A KN2 xvs (n1 )e xvs (n2 )e−j Ω0 m , − j A KN1 ∂Ω 3 3 m=0

m=0

(C.10) where K = ej [2a0 +a1 (n1 +n2 )+a2 (n1 +n2 ) ] . Substituting of the above results in (C.4) and (C.5) yields 2

−8A8 (N12 N26 + N16 N22 ) , 45   β ≈ −2A6 N1 N2 Im{Γ } , α≈

(C.11) (C.12)

where Γ = KN2

N1  m=0



N2  N12 ∗ N22 ∗ 2 j Ω0 m 2 2 m − m − xvs (n1 )e xvs (n2 )ej Ω0 m . + KN1 3 3 2

(C.13)

m=0

Using (C.11) and (C.12), we derive the first-order approximation for δΩ, δΩ = −

45[Im{Γ }] β ≈− . α 4A2 (N1 N25 + N15 N2 )

(C.14)

Its expected value, the bias of the a2 estimate, can be shown to be (to first-order approximation) zero. To compute the mean-square of (C.14), it is necessary to compute the values of E{Γ Γ ∗ } and E{Γ Γ }. For brevity, only the key formulae are listed:

668

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

Fig. 11. The theoretical MSE of the a2 estimate as function of n1 /N (n2 = 0).



23 5 2 4 2 5 7 6 4 2 5 4 4 5 2 N N + N N − N N2 + σ N N + N N E{Γ Γ } ≈ 2A σ 90 1 2 45 1 2 90 1 45 1 2 45 1 2

8 6 1 7 N N2 − N12 N25 + N1 N26 δ(N − 4n1 ), + 2A2 σ 2 45 1 6 90 E{Γ Γ } = 0, ∗



2 2

where u(·) denotes the unit step function. Hence we get



4 4 7 7 8 3 4 4 7 2 14 2 23 7 4 12 4 4 7 4 N N + N N − N N + 2A σ N N + N N E{β } ≈ 4A σ 90 1 2 45 1 2 90 1 2 45 1 2 45 1 2

8 8 3 1 4 7 7 N1 N2 − N1 N2 + N13 N28 δ(N − 4n1 ). + 4A14 σ 2 45 6 90

(C.15) (C.16)

(C.17)

Using (C.3), (C.11), (C.17), and δa2 = δΩ/2, we get  45[(46N15 N22 + 16N12 N25 − 14N16 N2 ) + (8N15 N22 + 8N12 N25 )/SNR]  E (δa2 )2 ≈ 256SNR(N1 N25 + N15 N2 )2 +

45(32N16 N2 − 30N12 N25 + 14N1 N26 ) 256SNR(N1 N25 + N15 N2 )2

δ(N − 4n1 ),

(C.18)

where SNR = A2 /σ 2 . It is shown that the variance of aˆ 2 depends on the values of N1 and SNR (n2 = 0 results in N2 = N/2). For any given N , SNR, E{(δa2 )2 } can be minimized by approximately choosing n1 . Figure 11 shows the theoretical MSE of the a2 estimate at SNR = 10 and −2 dB. From this plot, n1 = 0 gives rise to asymptotic minimum variance of a2 given as  90(1 + 1/2SNR)  , (C.19) E (δa2 )2 = N 5 SNR however, this choice leads the PCPF to be the square of CPF. Recall that the selected time positions should be sufficiently spaced to misalign the cross-terms and spurious peaks. Therefore, there is a trade-off between the estimated variance and the performance to analyze multicomponent signal. In this paper, we select n1 = 0.4N . Under this selection, the asymptotic MSE of the a2 estimate at high SNR is   91.6505 + 45.2152/SNR E (δa2 )2 ≈ . (C.20) SNRN 5

P. Wang, J. Yang / Digital Signal Processing 16 (2006) 654–669

669

By comparison, the CRLB for estimating the a2 is approximately   90 . (C.21) E (δa2 )2 ≈ SNRN 5 Thus, the asymptotic MSE of the a2 estimate at n1 = 0.4N and n2 = 0 is 1.83% above the CRLB at high SNR. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

L. Cohen, Time-Frequency Analysis, Prentice Hall, Englewood Cliffs, NJ, 1995. F. Hlawatsch, G.F. Boudreaux-Bartels, Linear and quadratic time-frequency signal representations, IEEE SPIE Mag. 9 (1992) 21–67. P.M. Djuric, S.M. Kay, Parameter estimation of chirp signals, IEEE Trans. Signal Process. 38 (1990) 2118–2126. T. Abotzoglou, Fast maximum likelihood joint estimation of frequency and frequency rate, IEEE Trans. Acoust. Speech Signal Process. 22 (1986) 708–715. X.-G. Xia, Discrete Chirp–Fourier transform and its application to chirp rate estimation, IEEE Trans. Signal Process. 48 (2000) 3122–3133. S. Peleg, B. Friedlander, The discrete polynomial phase transform, IEEE Trans. Signal Process. 43 (1995) 1901–1914. S. Peleg, B. Friedlander, Multicomponent signal analysis using the polynomial-phase transform, IEEE Trans. Aerospace Electron. Syst. 32 (1996) 378–387. M. Wang, A.K. Chan, C.K. Chui, Linear frequency-modulated signal detection using radon-ambiguity transform, IEEE Trans. Signal Process. 46 (1998) 571–586. O. Akay, G.F. Boudreaux-Bartels, Fractional convolution and correlation via operator methods and application to detection of linear FM signals, IEEE Trans. Signal Process. 49 (2001) 979–993. L.B. Almeida, The fractional Fourier transform and time-frequency representations, IEEE Trans. Signal Process. 42 (1994) 3084–3091. P. O’Shea, A new technique for instantaneous frequency rate estimation, IEEE Signal Process. Lett. 8 (2002) 251–252. P. O’Shea, A fast algorithm for estimating the parameters of a quadratic FM signal, IEEE Trans. Signal Process. 52 (2004) 385–393. S. Barbarossa, A. Scaglione, G.B. Giannakis, Product high-order ambiguity function for multicomponent polynomial-phase signal modeling, IEEE Trans. Signal Process. 46 (1998) 691–708. M.Z. Ikram, G.T. Zhou, Estimation of multicomponent polynomial phase signals of mixed orders, Signal Process. 81 (2001) 2293–2308. G.T. Zhou, M.Z. Ikram, Unsupervised detection and parameter estimation of multicomponent sinusoidal signals in noise, in: Proceedings of the 34th Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, CA, November 2000, pp. 842–846. S. Peleg, B. Porat, Linear FM signal parameter estimation from discrete-time observations, IEEE Trans. Aerospace Electron. Syst. 27 (1991) 607–614. R.A. Johnson, Miller & Freund’s Probability and Statistics for Engineers, Prentice Hall, Englewood Cliffs, NJ, 2004.

Pu Wang was born in Sichuan, China, in 1981. He received the B.S. degree in electronic engineering and the M.S. degree in communication and information systems, both from the University of Electronic Science and Technology of China (UESTC), Chengdu, China, in 2003 and 2006, respectively. He is currently a teaching assistant at the School of Electronic Engineering, UESTC, where he is pursuing his Ph.D. degree in signal and information processing. His current research interests include nonstationary signal processing, time-frequency analysis, detection and parameter estimation, robust estimation, and independent component analysis. He is a student member of IEEE. Jianyu Yang was born in Sichuan, China, in 1963. He received the B.S. degree from the National University of Defense Technology, Changsha, China, in 1984, and the M.S. and Ph.D. degrees from the University of Electronic Science and Technology of China (UESTC), Chengdu, China, in 1987 and 1991, respectively. He was an Associate Professor from 1993 to 1999, and has been a Professor since 1999 at UESTC. In 2005, he worked as a Visiting Professor in RLE of Massachusetts Institute of Technology (MIT). His current research interests include radar systems and imaging, passive millimeter wave imaging, signal detection and estimation, and digital signal processing. He is a senior member of CIE, a member of IEEE and IEE, a senior editor of the Chinese Journal of Radio Science, and a member of several academic committees in China.