Parameter estimate of multi-component LFM signals based on GAPCK

JID:YDSPR AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.1 (1-12)

Digital Signal Processing ••• (••••) ••••••

1

Contents lists available at ScienceDirect

67 68

2 3

Digital Signal Processing

4

69 70 71

5

72

6

www.elsevier.com/locate/dsp

7

73

8

74

9

75

10

76

11 12 13 14 15 16

Parameter estimate of multi-component LFM signals based on GAPCK ✩

77

Tong Gu, Guisheng Liao, Yachao Li ∗ , Yifan Guo, Yan Huang

79

82 83

a r t i c l e

i n f o

84

a b s t r a c t

85

19 20 21

Article history: Available online xxxx

22 23 24 25 26 27

80 81

National Laboratory of Radar Signal Processing, Xidian University, Xi’an 710071, China

17 18

78

Keywords: Linear frequency modulated signal The lower-order non-linearity Parameter estimates Generalized adjustable parameter correlation kernel

28 29 30 31

In this paper, we propose a fast and robust parameter estimate method for multi-component linear frequency modulated (LFM) signals from a finite number of noisy discrete time observations. First, a new kernel called generalized adjustable parameter correlation kernel (GAPCK) is introduced to avoid the coupling terms between time and lag. Then, Matched Fourier transform (MFT) and Robust Energy accumulation (REA) are utilized to obtain the chirp rate estimate of the received LFM signals. This estimate is used to compensate the time-quadratic phase term of the GAPCK, and then Fast Fourier transform (FFT) along time-axis is performed to estimate the constant coefficient. Moreover, the asymptotic statistical properties of parameter estimates are derived. The proposed method has low computational complexity and favorite performance under low signal-to-noise ratio (SNR) due to the loworder non-linearity of the GAPCK. Finally, simulated and real data are provided to verify the robustness and effectiveness of the proposed method. © 2020 Elsevier Inc. All rights reserved.

86 87 88 89 90 91 92 93 94 95 96 97

32

98

33

99

34

100

35

1. Introduction

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54

In recent years, parameter estimate of linear frequency modulated (LFM) signals have received considerable attention due to its extensive application in various fields, such as radar, sonar, communications, and seismic analysis [1–4]. In Inverse synthetic aperture radar (ISAR) applications, due to the relative motion between radar and target, the received signals can be modeled as LFM signals whose parameters reveal the motion characteristics about the target, i.e., velocity and acceleration [5–7]. To estimate parameters of LFM signals, various approaches have been proposed, such as the maximum likelihood estimation (MLE) [8–11], the time-frequency (TF) analysis and their extensions [12–28]. The MLE is statistically optimal among proposed methods, but a two-dimensional (2-D) search is required, which increases the computational complexity. To avoid 2-D search, several suboptimal methods [9,10] based on MLE are proposed, which can be fast implemented. However, these methods are only suitable for parameter estimate of mono-component LFM signal. To estimate parameters of multi-component LFM signals, the Fractional Fourier

55 56 57 58 59 60 61 62 63 64 65 66

✩

This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grants 61001211, 61303035, 614701283, the Key R&D Program of Shaanxi Province under Grant 2017ZDXM-SF-094, the National Key R&D Program of China under Grant 2018YFB2202500. Corresponding author. E-mail addresses: [email protected] (T. Gu), [email protected] (G. Liao), [email protected] (Y. Li), [email protected] (Y. Guo), [email protected] (Y. Huang).

*

https://doi.org/10.1016/j.dsp.2020.102683 1051-2004/© 2020 Elsevier Inc. All rights reserved.

transform (FRFT)-based method is proposed in [11]. This method converts LFM signal into one-dimensional (1-D) fractional order space for achieving dimensionality reduction, but it still requires a 1-D search. For the TF analysis, it can be generally divided into two categories: Linear TF transform and Non-linear TF transform. The Linear TF transform, such as the Short-time Fourier transform (STFT), Gabor transform [19], and Wavelet transform, can not simultaneously obtain a good resolution in both time and frequency domains. To achieve a better time-frequency resolution, Non-linear TF transform is investigated, such as Radon/Hough-Wigner transform (RWT/HWT) [16,17], Polynomial Wigner-Ville distribution (PWVD) [20–22], Second-order Wigner-Ville distribution (2D-WVD) [23], Polynomial-Phase transform (PPT) [24], and Product Cubic Phase function (CPF/PCPF) [25,26]. In RWT/HWT, Cartesian coordinate is converted into Polar coordinate, which is indispensable and has great computational burden. The 2D-WVD has low computational complexity, but it involves four-order non-linearity, increasing the number of the cross-terms. Therefore, the parameter estimate performance degrades severely when the signal-to-noise ratio (SNR) is under 0 dB. The CPF-based methods are asymptotically efficient for parameter estimate of mono-component LFM signal, but it presents spurious peaks whose energy is the twice of true peaks for parameter estimate of multi-component LFM signals. Moreover, the energy distribution of CPF is non-homogeneous, which will influence the robustness of the CPF-based methods. To solve the aforementioned problems, a fast and robust parameter estimate method based on generalized adjustable param-

101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:YDSPR

AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.2 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

eter correlation kernel (GAPCK) is proposed in this paper. It only involves third-order non-linearity and can accurately estimate parameters of multi-component LFM signals under low SNR. The proposed method can be summarized into three steps: First, the operation of GAPCK is implemented which not only eliminates the coupling terms between time and lag, but also has low-order nonlinearity, introducing fewer cross-terms. Second, two-step accumulation processing based on Matched Fourier transform (MFT) and Robust Energy accumulation (REA) is adopted to estimate chirp rates of multi-component LFM signals and suppress the impact of the cross-terms and noise. Third, using the chirp rate estimates, the compensate function is constructed to eliminate the time-quadratic phase term of the GAPCK, then Fast Fourier transform (FFT) along time-axis is implemented to estimate the constant coefficient. The remainder of this paper is organized as follows. Section 2 discusses Non-linear TF analysis of LFM signal. Section 3 presents the fast and robust parameter estimate method for multicomponent LFM signals based on GAPCK. Section 4 gives the computational complexity and statistical property of the proposed method. Section 5 validates the robustness and effectiveness of the proposed method by the simulated and real data. Finally, the conclusion is drawn in Section 6.

23 24

2. TFD of LFM signal

25 26 27

Assume x (t ) is made up of K-component noiseless LFM signals with constant amplitude, i.e.,

28 29 30

x (t ) =

33 34 35 36



A i · exp j2π ai ,1 t + ai ,2 t 2

41 42 43 44

(q)

K x (t , τ ) =

47 48 49

52 53 54 55 56 57



Si =

60 61 62 63 64 65 66

(2)

1

(3)

−1 ⇔ complex conjugate

Note that the WVD kernel is a special case of the MLICK function. The auto-term of WVD kernel can be represented as

WVD2auto (t , τ ) = x (t − τ ) · x∗ (t + τ )

=

K 

(4)

i =1

The WVD kernel has two-order non-linearity, which means the number of the cross-terms is less. However, the coupling term between t and τ exists in WVD kernel. Performing FFT with respect to τ , we can obtain

Y (t , f τ ) =

K 

2





| A i | δ f τ − 2ai ,1 + 4ai ,2t



2D − WVD4auto (t , τ ) = x t +

 

× x t− =

K 

τ 2





(5)

i =1

where f τ denotes the linear coefficient variable, δ (•) denotes the Dirichlet function. From (5), the LFM signal is distributed as multiple straight lines in the TF domain. Thus, it can only obtain the instantaneous frequency (IF) of the LFM signal and is hard to directly estimate parameters of LFM signal due to the existence of



− τ0 x∗ t +

τ 2

τ 2



x∗ t −

+ τ0



2

71

70

(6)

72 73 74

| A i | exp (− j4π ai1 τ0 − j8π ai2 τ0 t )

75 76

where τ0 is a fixed lag coefficient. From (6), one can see that the coupling terms between t and τ have been effectively removed. Therefore, FFT can be operated directly as

Y ( ft ) =

=

68 69

i =1

K 

67

τ

4

77 78 79 80 81

(4 )

82

K x_auto (t , τ ) exp (− j2π f t t ) dt



| A i |4 exp (− j4π ai1 τ0 ) δ f t − 4ai ,2 τ0

83



(7)

84 85

i =1

86

The peaks will produce at position f t = 4ai ,2 τ0 . However, the impact of the cross-terms become more serious due to the fourorder non-linearity of the 2D-WVD kernel. When x (t ) contains K-components, the number of the cross-terms is K 4 − K , which will increase SNR threshold for anti-noise performance according to [30]. To overcome these shortcomings, we propose a fast and robust parameter estimate method based on GAPCK in the following. 3. Parameter estimate of multi-component LFM signal based on GAPCK

87 88 89 90 91 92 93 94 95 96 97 98 99

3.1. Parameter estimate of multi-component LFM signals

100

From above discussion, we know that the parameters of LFM can not be directly estimated due to the existence of the coupling terms. To avoid the presence of coupling terms, the high-order correlation kernel is proposed, but it will bring more cross-terms. Thus, the coupling terms and high-order non-linearity is a mutual restrictive relationship, so how to optimize their relationship is the key to parameter estimate. In this paper, we propose a fast and robust parameter estimate method based on GAPCK, which not only avoids the coupling terms, but also has less cross-terms. The auto-terms of K-component LFM signals based on GAPCK can be expressed as

101 102 103 104 105 106 107 108 109 110 111 112 113

GAPCK 3auto (t , τ ) =



| A i |2 exp( j2π 2ai ,1 τ + 4ai ,2t τ )

58 59

x(si ) (t + u i · τ + p i )

where τ denotes lag, u i , p i are the coefficients to be determined, q is the order of the MLICK, and (si ) represents

50 51

 i =1

45 46

(1)

q

39 40



where t and A i denote time and constant amplitude, respectively, ai ,1 and ai ,2 denote the constant coefficient and chirp rate of the ith-component, respectively. The multi-linear instantaneous correlation kernel (MLICK) can be written as

37 38



i =1

31 32

K 

the coupling terms. To overcome this shortcoming, the high-order kernel function is proposed, such as the 2D-WVD kernel which is

⎧ ⎪ ⎨

K 

114

| A i |3

115

i =1  2 u 1 s1 + u 22 s2 + u 23 s3 ai2 τ 2 + (s1 + s2 + s3 ) ai2 t 2

⎞⎫ ⎪ ⎬ ⎜ +2 (u 1 p 1 s1 + u 2 p 2 s2 + u 3 p 3 s3 ) ai2 τ ⎟ · exp j2π ⎝ ⎠ + (u 1 s1 + u 2 s2 + u 3 s3 ) ai1 τ + 2 (u 1 s1 + u 2 s2 + u 3 s3 ) ai2 t τ ⎪ ⎪ ⎩ ⎭ ⎛

+2 ( p 1 s1 + p 2 s2 + p 3 s3 ) ai2 t + (s1 + s2 + s3 ) ai1 t

116 117 118 119 120

(8)

121

In (8), the coupling terms is a function of u i and p i . To avoid the coupling terms, these coefficients should satisfy the following conditions

123



u 1 s1 + u 2 s2 + u 3 s3 = 0 u 21 s1

+ u 22 s2

+ u 23 s3

or u 1 p 1 s1 + u 2 p 2 s2 + u 3 p 3 s3 = 0

122 124 125 126

(9)

Since the phase terms of τ do not contain the information of ai1 , estimating ai2 along lag-axis is convenient. For simplicity, one set of solutions to (9) is considered, which is

127 128 129 130 131 132

JID:YDSPR AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.3 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

1 2 3

⎧ s1 = 1 ⎨ u1 = 1 p1 = 1 u 2 = 1 p 2 = −1 s2 = 1 ⎩ u3 = 2 p3 = 0 s3 = −1

4

Substituting (10) into (8) yields

5

K 

6 7

GAPCK 3auto (t , τ ) =

8 9 10 11 12 13 14

⎧ ⎪ ⎪ ⎨

=

19 20 21

 T1

24 25 26 27 28 29 30 31 32 33 34 35 36 37

f τ 2 , t =

⎛

(11) x (nt ) =

" " fs "aˆ i1 " ≤ ,

| A i |3 i =1 ⎤⎞ 2

2

46



H com = exp − j2π aˆ j2 t 2

T2

=



f τ 2 , t = T 1

K  i =1,i = j









50 51

55 56 57 58 59 60 61 62 63 64 65 66

T3



Performing FFT on T 2

f









f τ 2 , f t = FFTt T 2

τ2

, t

2

(14)



f τ 2 , t

with respect to t, we can obtain

(17)

2N

⎧ ⎫ K ⎨  | A |3 δ  f − a + 2 a − a t ⎪ ⎬  ⎪ t i i1 i2 j2 = δ f τ 2 − (−2) · ai2 i =1,i = j ⎪ ⎪ ⎩ " "3  ⎭ + " A j " · δ f t − a j1

(15) where f t is the constant coefficient variable. Then we can estimate aˆ j1 and aˆ j2 in chirp rate-constant coefficient plane, simultaneously. Different from the fixed high-order kernel, the GAPCK adopts the optimal constraints between multi-parameters to determine its final model. The orders of GAPCK is lower than that of 2D-WVD kernel. That means GAPCK will have better performance under low SNR due to its low-order non-linearity. Therefore, we can estimate parameters more accurately under low SNR. Furthermore, although

83 84 85 86 87 89 90 91 93 94 95

= x (nt + nτ + 1) · x (nt + nτ − 1)

96

(18)

97 98 99

where nτ represents the discrete lag. s2 Performing MFT along lag-axis yields

100 101

 T 1 nt ,

f

τ2





−1)/2 ( N

=

102

GAPCK 3auto (nt , nτ )

103 104

i =−( N −1)/2

  · exp − j2π f τ 2 n2i

105

(19)

Note that this step can be computed by the fast algorithm [29] and the computational complexity of this step will be reduced from  



O N

2

to O N log2N .

106 107 108 109 110 111 112

s3 Extracting the peak location

113 114

For mono-component LFM signal, there are no external crossare accumulated along nτ and the peak of "terms.  The auto-terms "

" " " T 1 nt , f τ 2 " is located at −2ai2 as shown in Fig. 1 (a). Therefore,



82

88

" " "aˆ i2 " ≤ f s

· x (nt + 2nτ )

(13)



49

54

81

∗



f τ 2 , t · H com

79 80

GAPCK 3auto (nt , nτ )

  # | A i | exp j2π ai2 − a j2 t + ai1 t δ f τ 2 − (−2) · ai2   " "3 + " A j " exp ( j2π ai1 t ) δ f τ 2 − (−2) a j2 3

78

(16)

s1 Initialize K = 1 and transform x (nt ) with the GAPCK as



48

53

A i exp j2π ai1nt + ai2 nt2

92

Multiplying (12) by (13) yields



71

77



Steps of the proposed method can be summarized as follows:

accumulated along lag-axis and the peaks are located at −2ai2 without considering the impact of the cross-terms. Then, the peak extraction technology is applied to extract the chirp rate. Compensating function with respective to the jth estimated chirp rate aˆ j2 can be written as

47

52



where f τ 2 denotes the chirp rate variable. The auto-terms are

44 45



(12)

40

43

K 

where nt represents the discrete time, N represents the length of the signal, and the sampling rate is assumed to be f s . Hence, ai1  ai1 / f s and ai2  ai2 / f s2 . To avoid the ambiguities arising from the periodicity of the sampled signal, the following constraints should be satisfied [31]:

i =1

42

70

76

nt ∈ [0, N − 1]



−2ai2 τ     · exp ⎝ j2π ⎣ +ai2t 2 ⎦⎠ · exp − j2π f τ 2 τ 2 d τ 2 +ai1t K       | A i |3 exp j ai2t 2 + ai1t δ f τ 2 − (−2) · ai2 =

41

69

75

i =1

 K

⎡

38 39



| A i |3 · exp j2π −2ai2 τ 2 + ai2t 2 + ai1 t



68

74

The discrete form of LFM signals is

when x (t ) contains K-component, the number of the cross-terms is K 3 − K which is much smaller than that of 2D-WVD. Performing MFT on (11), we can get

22 23





i =1

17 18

K 

67

73

3.2. Realization

i =1

⎛

the orders of GAPCK are higher than that of WVD, the GAPCK avoids the existence of the coupling terms, which facilitates the execution of fast algorithm. Thus, the proposed method has advantages over these methods based on WVD kernel and 2D-WVD kernel.

72

| A i |3

⎞⎫ u 21 s1 + u 22 s2 + u 23 s3 ai2 τ 2 ⎪ ⎪ ⎜ ⎟⎬ ⎜ ⎟ · exp j2π ⎝ + (s1 + s2 + s3 ) ai1 t ⎠⎪ ⎪ ⎪ ⎪ ⎩ ⎭ + (s1 + s2 + s3 ) ai2 t 2

15 16

(10)

3

we can select arbitrary time cell to extract peak location. For multi-component LFM signals, there are two fatal problems. One is that the cross-terms is still inevitable, although the GAPCK avoids the coupling terms and has low-order non-linearity. Therefore, the impact of the cross-terms should be considered. For convenience, we consider a two-component LFM signals with simulation parameters in Table 1. Fig. 1 (b) shows the TF domain image of (19) and Fig. 1 (c) is the 677th cell of Fig. 1 (b) image. As shown in Fig. 1 (c), there are the spurious peaks in 677th cell, although the auto-terms of twocomponent LFM signals are accumulated along nτ . Thus, what we are most interested in is that why this phenomenon happens? Let us first consider the impact of the cross-terms with K-component LFM signals as

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:YDSPR

AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.4 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

4

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94 95

29 30 31 32 33 34 35

Fig. 1. (a) The chirp rate-time domain of auto term with mono-component LFM signal with SNR = 5 dB; (b) the chirp rate-time domain of two-component LFM signal with SNR = 5 dB; (c) the 677th cell of (b) image; (d) the chirp rate-time domain of two-component LFM signal with SNR = −5 dB; (e) the REA results of (d) image.

36 37 38 39 40 41 42

+

K K K   

44 45

Table 1 Simulation parameter for performance analysis.

∗

xi (t + τ + 1) x j (t + τ − 1) xl (t + 2τ )

i =1 j =1|i = j l=1 K K  

K K K   

⎛

49 50 51 52 53



⎡

ai1 + 2ai2 + a j1 − 2a j2 − al1 nt

57 58 59 60 61 62 63

R2 =

K 

K 

i =1 j =1|i = j

⎛

⎤⎞ ⎥⎟ ⎥⎟ ⎥⎟ ⎥⎟ ⎥⎟ ⎥⎟ ⎦⎠

A 2i A j

⎡







2ai1 − a j1 nt + 4ai2 − 4a j2 nt nτ

⎤⎞

⎥⎟  ⎜ ⎢  · exp ⎝ j2π ⎣ + 2ai1 − 2a j1 nτ + 2ai2 − 4a j2 n2τ ⎦⎠ 2  + 2ai2 − a j2 nt (22)

64 65 66

20 Hz 1034

102

1st constant coefficient ( a 11 )

0.484

104

1st chirp rate ( a 12 )

4.8356e-4

105

2nd constant coefficient ( a 21 )

0.0968

106

2nd chiro rate ( a 22 ) Signal-to-noise ratio

1.2e-3 5 dB/-5 dB

107





103



108 110





2ai2 + 2a j2 − 4al2 nt + ai1 + 2ai2 + a j1 − 2a j2 − 2al1 = NM

(21)

54 55

101

Sampling frequency Sampling numbers



⎜ ⎢ + 2a + 2a − 4a n n ⎜ ⎢ t τ i2 j2 l2 ⎜ ⎢  ⎜ ⎢ + a + 2a + a − 2a · exp ⎜ j2π ⎢ i1 i2 j1 j2 − 2al1 nτ  ⎜ ⎢ ⎝ ⎣ + ai2 + a j2 − 4al2 n2τ  + ai2 + a j2 − al2 nt2

48

Value

109

A i A j Al

i =1 j =1|i = j l=1

47

100

Parameter



The first component R1 and second component R2 can be expressed as

R1 =

99



i =1 j =1|i = j

46

56

(20)

xi (t + τ + 1) xi (t + τ − 1) x∗j (t + 2τ )

43

97 98

GAPCK 3cross (t , τ )

=

96

From (21) and (22), we can get that when the phase satisfies the relationship (M is any integer):









111 112

4ai2 − 4a j2 nt + 2ai1 − 2a j1 = NM

113

(23)

114

the spurious peaks will appear in the location of ai2 + a j2 − 4al2 and 2ai2 − 4a j2 . According to (23) and simulation parameters, Fig. 1(c) shows that the Peak 1 and Peak 4 are the spurious peaks caused by the cross-terms and the amplitude of them is similar to the true Peak 2 and Peak 3, which is consistent with the theoretical analysis. Thus, the robustness of the peak extraction of multi-component LFM signals should be considered due to the impact of the cross-terms. The other fatal problem is that the above analysis only considers the spurious peak caused by the cross-terms, while the high levels of background noise is ignored. From Fig. 1 (d) one can see that the result is worse with SNR = −5 dB after performing s1 step. Since the MFT operation is a 1-D energy accumulation, it is not enough to extract the true peaks. To address these problems, REA operation, which accumulates energy along nt is adopted as follows:

116

115 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:YDSPR AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.5 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

1 2

REA =

5 6 7 8 9

Y

f τ 2 , n i



or REA =

i =1

3 4

N  

N 

 Y1

f τ 2 , n i

 (24)

i =1

s4 Compensate (19) with the compensation function according to the extracted chirp rate in s3

H com (nt ) = exp



− jaˆ i2nt2





(25)

20





23

=



K 

&

" "3 " A j " · exp j2π

'

j =1|i = j

* + | A i |3 · exp ( j2π [ai1t])

26 27

FDQPT ICPF

 T3





f τ 2 , f t = FFTt T 2 3

= | Ai | δ



f





a j2 − aˆ i2 t +a j1 t

2

()

30

(26)

MLE

2D-WVD

2D-MF

ICPF

CPF

GAPCK

Runtime(s)

50.2

15.4

12.5

11.8

9.4

9.6

33 34 35 36 37

40 41 42 43 44 45



f τ 2 , nt





(27)

τ 2 − (−2) ai2 · δ ( f t − ai1 )

48 49

Parameter estimate

Asymptotic MSE

54 55 56 57 58 59 60 61 62 63 64 65 66

76 77 78

82 83 84

aˆ 1

6 2 ·SNR· N 3

aˆ 2

90 4 ·SNR· N 5

 

+

1 SNR2

1.5 1.5 + SNR +

0 .5 SNR2

1.5 +

2 SNR

 

6 2 ·SNR· N 3 90 2 ·SNR· N 5

88 89 90 91 92 94

location of the maximum value of T 3

95

f

τ 2 , ft .

96 97 98

s6 Estimate the amplitude of the ist component as

" " " " N −1 2 " "     " 2 ˆA i = 1 "" x (nt ) exp − j2π aˆ i1nt + aˆ i2 nτ " " N" N − 1 " "n=− 2

99 100 101

(28)

102 103 104

To suppress the impact of noise and weak objectives, the Relax [32] is adopted. Then, repeating S1-S6, the residual signal is





ˆ i exp j2π aˆ i1nt + aˆ i2nt2 xnew (nt ) = x (nt ) − A

 (29)

4. Performance analysis of the proposed method In this section, we mainly analyze the proposed method in terms of computational complexity and anti-noise performance in detail (i.e., asymptotic mean-square error, MSE). 4.1. Computational complexity

52 53

75

87

CRB

Then, aˆ i1 and re-estimate aˆ i2 can be extracted  according to the

50 51

74

93

46 47

73

86

38 39

72

85

Table 4 Theoretical MSE and CRB for estimated parameter.

31 32

71

81

Methods

28 29

70

80

Table 3 Comparison of runtime with different methods.

s5 Performing FFT with respect to nt , we can obtain

24 25

PPT

69

79



21 22

2D-WVD

GAPCK

T 2 nt , f τ 2 = T 1 nt , f τ 2 · H com (nt )

16

19

MLE

Then, multiplying (19) by (25) yields

14

18

 O L2  O K · ( N 2 + N 2 log2 N ) O (2K · LN log2 N )   O K · N 2 log2 N O ( K · N ( N − 1)) O ( K · ( N log2 N ))   O K · N2 O ( K · ( N log2 N ))

CPF

12

17

68

Computational complexity

2D-MF

11

15

67

Table 2 Comparison of computational complexity. Algorithm

From Fig. 1 (e), we can see that the true peaks are extracted accurately and the result is consistent with that in Fig. 1 (c).

10

13

5

Based on the solution of (10), the computational complexity of the proposed method contains four steps: the implementation of GAPCK (18), MFT (19), REA (24), and FFT (27). Owing to the fast algorithm, the computational complexity of performing MFT and  GAPCK will be reduced from O N 2 to O ( N log2 N ). The computational cost of REA and FFT is O ( N ) and O ( N log2 N ), respectively. Due to the influence of weak objectives, the Relax technique will be adopted. Thus, the total computational cost of the proposed method is O ( K · ( N log2 N )) (i.e. K is the iterative number or the number of LFM components). Compared with these existing methods (i.e., MLE [8], 2D-WVD, 2D-MF [33], CPF, PCPF, FDQPT, PPT, and ICPF [23–28]), the computational complexity and runtimes of different methods are listed in Tables 2 and 3, respectively. The simulation is performed with

Fig. 2. Estimate MSEs, (a) MSE of a2 , (b) MSE of a1 .

105 106

MATLAB R2016b environment, Intel Core i7-6700, 3.40 GHz processor, 16.0 GB RAM, and Microsoft Windows 7. The simulation parameters in Table 1 are adopted. In Table 2, L and K represent the searching time and the number of LFM components, respectively. The value of L is usually greater than N which is determined by the incremental step. In Tables 2 and 3, we can find that the computational complexity of the proposed method is much lower than the most methods. Although the computational complexity of CPF/PCPF is same as the proposed method, the estimation performance of CPF/PCPF will degrade severely under low SNR, which will be demonstrated in the subsequent statistical property analysis.

107 108 109 110 111 112 113 114 115 116 117 118 119

4.2. Statistical properties

120 121

The first-order statistical analysis [10] of parameter estimates is given in Appendix A and the results are shown in Table 4. The SNR is defined as SNR = A 2 /σ 2 . When SNR → ∞, the theoretical MSE of a2 and a1 is about 50 percent higher than the corresponding CRB [34] due to the three-order correlation property. In order to further evaluate the performance, the proposed method is compared with the CPF, PCPF, and ICPF. The simulation parameters are shown in Table 1 with N=256 and different SNRs (i.e., −8 dB ∼ 8 dB). For each SNR, 100 Monte Carlo simulations are performed. The MSEs of the estimated parameters are shown in Fig. 2.

122 123 124 125 126 127 128 129 130 131 132

JID:YDSPR

AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.6 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

6

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78 79

13 14

Fig. 3. The simulation with time-variant amplitude. (a), (b) and (c) by the proposed method with time-variant amplitude (initial SNR = 5 dB).

80 81

15 16 17 18 19 20 21 22 23 24

In Fig. 2, the estimation performance of the proposed method degrades severely until input SNR falls below −5 dB, while the estimation performance of CPF/PCPF begins to degrade severely when input SNR falls below 0 dB. That is to say the proposed method has a low SNR threshold. The SNR threshold of ICPF is almost the same as the proposed method. However, the computational complexity of ICPF is much greater than the proposed method as shown in Table 2. Therefore, the proposed method is a fast and robust parameter estimate method.

25 26

5. Simulation results

27 28 29 30 31 32 33 34 35 36

This section presents several numerical results obtained by the simulated data and real data to verify the effectiveness and robustness of the proposed method. The simulation parameters are shown in Table 5. The simulation verification is divided into three parts. The first part verifies the performance of multi-component LFM signals with time-variant/constant amplitude. Second, the anti-noise performance comparison with other methods is presented with different SNRs. Third, the effectiveness of the proposed method with real data is studied.

37 38 39

5.1. Performance analysis of multi-component LFM signals with time-variant/constant amplitude

40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

In this subsection, the capability of the proposed method to estimate parameters of multi-component LFM signals is verified. Assume that the amplitudes of multi-component LFM signals are constant with SNR = 5 dB. The simulation results and estimated parameters are shown in Fig. 4 (g), (h) and (i). Obviously, not only the peaks of multi-component LFM signals can be well focused, but also the value of parameters is estimated accurately. However, in many applications, each amplitude in multi-component LFM signals is usually different and varies with time. To verify the effectiveness of the proposed method, we assume that the amplitudes of multi-component LFM signals are A(n) = Aexp(−(n/N)2/2) (i.e. initial SNR = 5 dB with n = 0). Fig. 3 shows the simulation results with time-variant amplitude. We can find that the results are similar with those in Fig. 4 (g), (h) and (i). From Fig. 3, Fig. 4 (g), (h) and (i), we can see that the proposed method has a better parameter estimate performance for multi-component LFM signals whether the amplitude is constant or time-variant.

58 59 60

63 64 65 66

83

Parameter

Value

Sampling frequency Sampling numbers

200 Hz 257 

84 85 86

1st constant coefficient ( a 11 )

0.075

87

1st chirp rate ( a 12 )



2.5e-4

88

2st constant coefficient ( a 21 )

0.05

89

2st chirp rate ( a 22 )

1.25e-4

90

3st constant coefficient ( a 31 )

0.025

91

3st chirp rate ( a 32 ) Signal-to-noise ratio

3.75e-4 −5 dB/0 dB/5 dB

92









In this subsection, we verify the performance of the proposed method from the aspects of anti-noise. Moreover, the proposed method are compared with other methods under different SNRs (i.e., SNR = 5, 0, −5 dB). The comparison results are shown in Tables 6, 7 and 8.

93 94 95

Table 6 Estimated values in the simulation of SNR = 5 dB. Comp1

MLE 2D-WVD CPF/PCPF ICPF GAPCK

Comp2

96 97

Comp3

98

a11

a12

a21

a22

a31

a32

99

0.0745 0.0745 0.0747 0.0744 0.0744

2.519e-4 2.537e-4 2.522e-4 2.536e-4 2.517e-4

0.0501 0.0505 0.0510 0.0506 0.0511

1.253e-4 1.216e-4 1.234e-4 1.251e-4 1.261e-4

0.0248 0.0258 0.0255 0.0263 0.0267

3.769e-4 3.774e-4 3.762e-4 3.754e-4 3.724e-4

100 101 102 103 104 105 106

Table 7 Estimated values in the simulation of SNR = 0 dB. Comp1

MLE 2D-WVD CPF/PCPF ICPF GAPCK

Comp2

107 108

Comp3

a11

a12

a21

a22

a31

a32

0.0755 0.0766 0.0768 0.0748 0.0764

2.519e-4 2.544e-4 2.566e-4 2.544e-4 2.472e-4

0.0501 0.0512 0.0516 0.0508 0.0511

1.231e-4 1.215e-4 1.244e-4 1.257e-4 1.260e-4

0.0258 0.0267 0.0266 0.0253 0.0267

3.694e-4 3.724e-4 3.776e-4 3.759e-4 3.646e-4

109 110 111 112 113 114 115 116

Table 8 Estimated values in the simulation of SNR = −5 dB. Comp1

MLE 2D-WVD CPF/PCPF ICPF GAPCK

Comp2

117 118

Comp3

119

a11

a12

a21

a22

a31

a32

120

0.0755 0.0774 0.0944 0.0747 0.0764

2.481e-4 2.642e-4 6.536e-4 2.535e-4 2.475e-4

0.0521 0.0522 0.0426 0.0512 0.0520

1.193e-4 1.225e-4 3.251e-4 1.255e-4 1.183e-4

0.0239 0.0272 0.0463 0.0261 0.0267

3.769e-4 3.724e-4 7.754e-4 3.758e-4 3.684e-4

121

5.2. Performance analysis of multi-component LFM signals with different SNRs

61 62

82

Table 5 Simulation parameter for performance analysis.

122 123 124 125 126

In Tables 6, 7 and 8, we can see that the performance of MLE is relatively robust and the estimated results are closer to true values. However, it is time-consuming in practice due to the 2-D search. Compared with 2D-WVD and ICPF, the proposed method has lower computational complexity and can achieve similar estimate accuracy. The performance of CPF/PCPF degrades severely

127 128 129 130 131 132

JID:YDSPR AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.7 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

7

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

Fig. 4. The comparison with different SNRs. (a), (b) and (c) by the proposed method with SNR = −5 dB, (d), (e) and (f) by the proposed method with SNR = 0 dB, (g), (h) and (i) by the proposed method with SNR = 5 dB.

with SNR = −5 dB, which is consistent with the previous theoretical analysis. In contrast, the proposed method still has a better performance with SNR = −5 dB. Fig. 4 shows the distributions of multi-component LFM signals in chirp rate-constant coefficient domain with different SNRs based on GAPCK. In Fig. 4, with the decrease of SNRs, the amplitude of background noise increases rapidly. However, the estimate errors are slightly increased with the decrease of SNRs with the proposed method. Meanwhile, to avoid the produce of strong spurious peaks, the Relax technique is utilized in Fig. 4. To summarizes, the proposed method is suitable for parameter estimate with low SNR, and has great robust performance (i.e., it is robust against additive noise) and light computational burden, which is consistent with the aforementioned analysis.

58 59

5.3. Performance analysis of real data

60 61 62 63 64 65 66

In this subsection, the Ground/Airborne-based ISAR data is utilized. First, the phase parameter of Ground-based ISAR data can be considered to be approximately mono-component LFM signal due to the uniform rotation feature of airplane. Furthermore, to validate the robustness of the proposed method, we choose the 129th and 86th range samples with different energy level as shown in

Fig. 5. In the Ground-based ISAR data, the carrier frequency 5.5 GHz, bandwidth 400 MHz, pulse repetition frequency (PRF) is 100 Hz, and the azimuth sampling is 1000. Fig. 6 shows the simulation results of Ground-based ISAR data. In Fig. 6, these three images from left to right are GAPCK image, true TFD image, and the estimated TFD image, respectively. Fig. 6 (a), (d) shows the GAPCK images of two range samples. The energy of the two GAPCK images is consistent with that of images in Fig. 5. Compared with the other four images, the estimated TFD images shown in Fig. 6 (c), (f) are similar with Fig. 6 (b), (e), which is consistent with the property of uniform rotation. However, there exist some differences between the true TFD image and the estimated TFD image due to the effect of noise. Due to the complex three-dimensional motion (yaw, pitch and roll) of ship with the oceanic wave, the three-dimensional positions will be coupled in phase terms. The ISAR imaging model and processing procedure can be referred to [35]. Here, we only consider the parameter estimate problem. The Airborne-based data is utilized to validate the parameter estimate performance for multicomponent LFM signals. In Airborne-based data, the carrier frequency is 9.25 GHz, the bandwidth is 500 MHz and the PRF is 200 Hz, and the azimuth sampling is same as the Ground-based ISAR data.

106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:YDSPR

8

AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.8 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93 94

28 29 30

Fig. 6. The simulation with Ground-based ISAR data. (a), (d) are GAPCK images with 86th and 129th range sampling, (b), (e) are the true TFD images with 86th and 129th range sampling, (c), (f) are the estimated TFD images with 86th and 129th range sampling.

95 96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110 111

45 46 47

Fig. 7. The simulation with Ground-based ISAR data. (a) is the true TFD image, (b) is the estimated TFD image by MLE, and (c) is the estimated TFD image by the proposed method. Table 9 Estimated values with Ground-based data.

49 50

115 116

51

Estimate parameter

ai1 /PRF

ai2 /PRF2

52

MLE Proposed method

−1.657 −1.657

−1.6266e-6 −1.2513e-6

53 55 56 57 58 59

Fig. 5. Ground-based ISAR data from airplane with uniform rotation.

61 63 64 65 66

117 118 119 120

54

62

113 114

48

60

112

In Fig. 7, these three images from left to right are the true TFD image, estimated TFD image by the proposed method, and the image obtained by MLE, respectively. Due to the 2-D search, we can approximately believe that the estimated parameters of MLE are true. The results in Fig. 7 (b) and (c), are similar and are approxi-

mately consistent with that in Fig. 7(a). Thus, we can conclude that the proposed method has better parameter estimate performance for multi-component LFM signals and is applicable in practice. Noted that the detailed estimation results are listed in Tables 9 and 10 with the form of aˆ i /PRF i (i = 1, 2) [36].

121

6. Conclusion

127

122 123 124 125 126 128

In this paper, a fast and robust parameter estimate method for multi-component LFM signals is proposed. Firstly, the GAPCK is constructed, which is utilized to transform LFM signal into the time-lag domain. Secondly, the operation of MFT and REA are

129 130 131 132

JID:YDSPR AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.9 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

1 2

Table 10 Estimated values with Airborne-based data.

+

3

ai1 /PRF

4

MLE

Proposed

MLE

Proposed

−1.634 −1.775 −1.458

−0.1617 −0.1747 −0.1466

3.1362e-6 1.9264e-5 −1.2536e-5

3.1281e-6 1.9144e-5 −1.2638e-5

5 6 7

Comp1 Comp2 Comp3

ai2 /PRF2

8

9

⎧ S n+m+1 S n+m−1 W n∗+2m + S n+m+1 W n+m−1 S n∗+2m ⎪ ⎪ ⎪ ⎪ ⎪ ∗ ⎨ + S n+m+1 W n+m−1 W ∗ n+2m + W n+m+1 S n+m−1 S n+2m

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬

67 68 69

⎪ + W n+m+1 S n+m−1 W n∗+2m + W n+m+1 W n+m−1 S n∗+2m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ + W n+m+1 W n+m−1 W n∗+2m  

= b30 exp j −2a2 2m2 + a1 n + an2



70 71 72 73

+ β (n, m)

(34)

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

taken to estimate the chirp rate. Thirdly, after compensation operation, the FFT along time-axis is utilized to transform the signal into the chirp rate-constant coefficient domain. Finally, the constant coefficient is estimated according to the location of the maximum value. The advantages of the proposed method are summarized as follows: 1) The proposed method eliminates the coupling term and has no search operation. Hence, it avoids the operation of decoupling and has light computational burden; 2) Due to the lowerorder non-linearity of the GAPCK, the proposed method has fewer cross-terms and better anti-noise performance. In addition, when the original signal contains multi-components, each component of the multi-component LFM signals can be reconstructed combining with the Relax technique. The theoretical analysis, simulated and real data show the effectiveness and robustness of the proposed method.

Since 0 ≤ n ≤ N − 1 and

⎧ ⎧ ⎪ ⎪ ⎨ −n − 1 ≤ m ≤ N − n − 2 ⎨0≤n+m+1≤ N −1 0 ≤ n + m − 1 ≤ N − 1 ⇒ −n + 1 ≤ m ≤ N − n ⎪ ⎪ ⎩ ⎩ 0 ≤ n + 2m ≤ N − 1 −n/2 ≤ m ≤ N − n − 1/2

30 31 32 33 34 35 36 37 38 39

Appendix A. First order perturbation analysis This appendix presents the derivations of the first order statistical results of the phase parameters and some key formulae were derived in [34]. Assume g N () is a complex function with N samples, which depends on a real variable . The maximum value of g N () is at  = 0 . δ g N () is the perturbed function, which will make a small offset for 0 (i.e., 0 + δ ). The MSE of δ  can be approximately expressed as



42 43 44



E (δ )2 ≈

40 41

where

+

E D

2

, (30)

C2



∂ 2 g∗N (0 ) ∂ gN (0 ) ∂ g∗N (0 ) ∂ 2 f N (0 ) = 2 Re g  + ( ) N 0 ∂ ∂ ∂ 2 ∂ 2

#

45 46 47 48 49

(31)

*  ∂δ g∗N (0 ) ∂ gN (0 ) ∗ ∂δ f N (0 ) = 2 Re gN (0 ) + δ gN (0 ) ∂ ∂ ∂

(32)

50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

77 78

(35)

N −n−1/2



gN () =

m=−n/2



E {•} denotes expectation operator.

81 82 83

S n+m+1 S n+m−1 S n∗+2m exp − j 2m2

84

 (36)

89

N −n−1/2



δ gN () =

  β (n, m) exp − j 2m2

90

(37)

92

The maximum value of g N () is at 0 = −2a2 . The detailed derivations are given by

 

N −n−1/2



g N (0 ) =

b30 exp j a1 n + an2



99

(38)

  ∂ gN (0 ) = − j 2 b30 exp j a1 n + an2 ∂ & ) 3  3 # 2 b30 K N −n−1 n ≈ −j +   ∂ 2 gN (0 ) = −4 b30 exp j a1 n + an2 2 ∂ & )5  5 # −4 b30 K N −n−1 n ≈ +

 N −n−1/2

2



m=−n/2

102

(39)

103 104 105 106



107 108

m4

m=−n/2

109

(40)

110 111

2

112



N −n−1/2

100 101

m2

2

2



95

98

 N −n−1/2

5

94

97

≈ b30 K {( N + 1) /2}

3

93

96

m=−n/2

δ gN (0 ) =

91

β (n, m) exp − j 0 2m2

113

 (41)

m=−n/2

114 115 116



N −n−1/2



m2 β (n, m) exp − j 0 2 m2



117

(42)

118 119

m=−n/2

120

 

Assume the original signal is given by 2

Z n = S n + W n = b0 exp j a1 n + a2 n 

2



+ Wn

K = exp j a1 n + an2

121



(43)

(33)

After GAPCK operation, the signal can be expressed

Z GAPCK = Z n+m+1 Z n+m−1 Z n∗+2m

+ , = { S n+m+1 + W n+m+1 } { S n+m−1 + W n+m−1 } S n∗+2m + W n∗+2m = S n+m+1 S n+m−1 S n∗+2m

86 88

The perturbation in g N () is

where the constant K is defined by

 

85 87

∂δ gN (0 ) = − j 2 ∂

1) Asymptotic MSE of a2

79 80

m=−n/2

28 29

76

the intersection in (35) is −n/2 ≤ m ≤ N − n − 1/2. The function g N () is given by

26 27

74 75

9

123 124

Substituting (38)–(43) into (31) and (32) yields

⎧' ( -& ) 5  5 . ⎫ ⎪ ⎪ N +1 N −n−1 n ⎪ ⎪ ⎪ + /5 ⎪ ⎪ ⎪ ⎨ ⎬ 2 2 2 2 ∂ f N (0 ) 6 4 . ≈ 2b  2 ) & 0 3  n 3 ⎪ ⎪ ∂ 2 N −n−1 ⎪ ⎪ ⎪ ⎪ ⎪− ⎪ + /9 ⎩ ⎭ 2

122

2

125 126 127 128 129 130 131

(44)

132

JID:YDSPR

AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.10 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

10

1 2 3 4 5 6 7 8 9 10

∂δ f N (0 ) ∂ ⎧ ⎡ & ) ⎤⎫ N +1 ⎪ ⎪ 2 ⎪ ⎪ m ⎪ N − n − 1 / 2 ⎪ ⎥⎪ ⎪  ⎢ ⎪ ⎪ 2 ⎢ ⎥ ⎪ ⎪ . ⎨ ⎬ & )3  3 ⎢ ⎥ 3 2 N − n − 1 n ⎣ ⎦ = 2b0  Im m=−n/2 − + 3 / ⎪ ⎪ ⎪ ⎪ 2 2 ⎪ ⎪ ⎪ ⎪   ⎪ ⎪ ⎪ ⎪ ⎩ · K · β ∗ (n, m) · exp j 0 2m2 ⎭ =

2b30 2 Im {

13 14 15

18

⎛ E

+

ηη

∗

,

⎞

3b40 σ 2

⎜ ⎟ ≈ ⎝ +3b20 σ 4 ⎠ +σ 6

20 22 23 24 25 26



27

E

2



2

2

E

35 36 37 38 39 40

∂δ f N (0 ) ∂

(2 #

42



E (δa2 )2 ≈

≈

1 4

&

47

σ

56 57 58 59 60 61

66

79

720 ∂δ f N (0 ) ≈ 3 Im {η} ∂ b 0 2 N 6

H com = exp − jaˆ 2  n + j2aˆ 2  m 2 2

(51)

2

2

1 44 SNR

1.5 +

E 1.5 SNR

∂δ f N (0 ) ∂ +

0.5

84 85

= exp − j (a2 + δa2 ) 2n2 + j (2a2 + 2δa2 ) 2m2



(52)

88 89

(1 )



(1)

(1 )

Zn = Sn + W n





2 2

exp − j δa2  n



90

  exp − j2δa2 2m2

91 92

95

∂ 2 f N (0 ) ∂ 2

90

4 · SNR · N 5

1.5 +

  exp − j2δa2 2 m2 ≈ 1 − j2δa2 2m2

104

SNR

+

0.5 SNR2

(56)

(1)

2

(1) m2 S n

N 

δ gN () =

(50)

2) Asymptotic MSE of a1 The analysis of aˆ 1 is considered under the assumption that m = 0 or m = N. In this case, δa2 can be approximately given by

105 107

(1 )



2 2

+ W n exp − j δa2  n

108

 (57)

109 110 111 112

(1)

Sn exp {− j n}

(58)

  + , N (1) (1 )  W n exp − j δa2 2 n2 − j2δa2 2m2 S n exp {− j n}

N 

114 115 116 117

(59)

118 119

The maximum value of g N () is at 0 = a1 . The detailed derivations are given by

g N (0 ) =

103

113

n =1

*

102

106

n =1

,

99 101

gN () =

1.5

98

(55)

The functions g N () and δ gN () are given by

(49)

97

100

Z n = S n − j2δa2 

*− 2

SNR2

+



(54)

Since (15) or (27) is operated along n, the phase about m can be seen as constant. Hence, the influence of m can be ignored. Meanwhile, the following approximation holds for larger N

(1 )

2 For any given values of , N, + and ,SNR, E (δa2 ) depends on the choice of n. Fig. 8 shows E (δa2 )2 as a function of n/ N with SNR 1. It can be seen that the minimum value appears at m = 0 or m = N, and it can be written as



93 94

Hence, (53) can be rewritten as

⎫⎫ ⎧ ⎧ 2 ⎬⎬ ⎨& N + 1 ) ⎨  N −n−1 3 +  n 3 /9 2 2   /   5 N −n−1 5 ⎩ − N +1 ⎩ 2 + n2 /5 ⎭⎭ 2 2

E (δa2 )2 ≈

86 87

the output is

)



81 83

    (1) W n = β exp − ja2 2n2 exp ja2 2m2

2

(2 # 

80 82





120 121 122 123

b30 exp { ja1 n} exp {− j 0 n} ≈ b30 N

(60)

n =1

≈ −j

124 125 126

∂ gN (0 ) = − j ∂

64 65

78

(

(1 )

62 63

∂ 2 f N (0 ) δ a2 ≈ − ∂ 2

(− 1 '

S n = b30 exp { j [a1 n]}

2

54 55

'

96

⎞

2

⎜ ⎟ ≈ 24 b60 ⎝ +3b20 σ 4 ⎠

 * '



46

53

3b40

Finally, using (45), (48), and δa2 = δ / (2τ ), we can get

43

52

⎛

2

41

51

77

Fig. 8. MSE error of aˆ 2 as function of m/ N.

(47)

+σ 6 ⎫⎫ ⎧⎧ & ) & )3  3 .2 ⎪ ⎪ ⎪ N + 1 N − n − 1 n ⎪⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ + /9 ⎪ ⎪ ⎪ ⎪ ⎬⎪ ⎨⎨ ⎬ 2 2 2 × (48) . )2 & )5  5 & ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ N +1 N −n−1 n ⎪ ⎪ ⎪ ⎪ ⎪ + /5 ⎪ ⎭⎪ ⎩⎪ ⎭ ⎩−

34

50

76

(53)

η2 ≈ 0

'

32

49

74

Then, we can get

31

48

73

where

30

45

72



29

44

71

After the compensation operation as

⎧& )2 -& )3  3 .2 ⎫ ⎪ ⎪ N +1 N −n−1 n ⎪ ⎪ ⎪ ⎪ ⎪ + /9 ⎪ ⎪ ⎪ ⎬ ⎨ 2 2 2 (46) × ) -& )5  5 . ⎪ & ⎪ ⎪ ⎪ ⎪ ⎪ N + 1 N − n − 1 n ⎪− ⎪ ⎪ + /5 ⎪ ⎭ ⎩

21

33

70

75

where η is clearly defined in the above expression. To compute the mean square + , value of (45), we first compute the values of E {ηη∗ } and E η2 .

19

28

69

(45)

16 17

68

η}

11 12

67

N 

127

nb30 exp { ja1 n} exp {− j 0 n}

128 129

n =1

b30 N 2  2

130

(61)

131 132

JID:YDSPR AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.11 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

1 2 3 4 5

 ∂ 2 gN (0 ) = − j · − j 2 b30 n2 exp { ja1 n} exp {− j 0 n} 2 ∂ N

n =1

≈−

b30 N 3 2

6 7 8

3

δ gN (0 ) ≈

11

14 15 16

21 22 23

3

∂ f N (0 ) ≈− ∂ 2

27 28 29 31 32



'

E

∂δ f N (0 ) ∂

(1)

nW n exp − j δa2 2n2

(64)

2δ a 2 4

3 b30 N 4

-

(65)

 N & 

120b30 N2

n−

n =1

.

N

)

2

#

35 36

S∗Wn

(1 )

Im (η)

(66)

(2 #

≈

2 b60  12

3b40

σ

2

+ 4b20

4

σ + 2σ

6

 (67)

37

Due to δa1 = δ / (t ), the final result is given by





E (δa1 )2 ≈

6

2 · SNR · N 3



1.5 +

2 SNR

+

1 SNR2

*

(68)

38 39 40

Declaration of competing interest

41 42 43

The authors declared that they have no conflicts of interest to this work.

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

71

References [1] X.G. Xia, Discrete chirp-Fourier transform and its application to chirp rate estimation, IEEE Trans. Signal Process. 48 (11) (2000) 3122–3133. [2] P.-L. Shui, Z. Bao, H.-T. Su, Non-parametric detection of FM signals using timefrequency ridge energy, IEEE Trans. Signal Process. 56 (5) (2008) 1749–1760. [3] J. Wang, H.-T. Wang, Y. Zhao, Direction finding in frequency-modulated-based passive bistatic radar with a four-element adcock antenna array, IET Radar Sonar Navig. 5 (8) (2011) 807–813. [4] A. Sommer, J. Ostermann, Backprojection subimage autofocus of moving ships for synthetic aperture radar, IEEE Trans. Geosci. Remote Sens. (2019) 1–11, Early access. [5] Y. Gao, M.D. Xing, Isar imaging and cross-range scaling for maneuvering targets by using the NCS-NLS algorithm, IEEE Sens. J. 19 (13) (2019) 4889–4897. [6] Y. Wang, X.F. Chen, 3-D InISAR imaging of the ship target based on joint cross S-method algorithm, IEEE Trans. Geosci. Remote Sens. Lett. 16 (7) (2019) 1080–1084. [7] L. Wu, X.Z. Wei, ISAR imaging of targets with complex motion based on discrete chirp Fourier transform for cubic chirps, IEEE Trans. Geosci. Remote Sens. 50 (10) (2012) 4201–4213. [8] T. Abotzoglou, Fast maximum likelihood joint estimation of frequency and frequency rate, IEEE Trans. Aerosp. Electron. Syst. 22 (6) (1986) 708–715. [9] P.M. Djuric, S. Kay, Parameter estimation of chirp signals, IEEE Trans. Signal Process. 38 (12) (1990) 2118–2126. [10] S. Peleg, B. Porat, LFM signal parameter estimation from discrete-time observations, IEEE Trans. Aerosp. Electron. Syst. 27 (4) (1991) 607–614.

74 76 77 78

[17] Y. Sun, P. Willett, Hough transform for long Chirp detection, IEEE Trans. Aerosp. Electron. Syst. 38 (2) (2002) 553–569.

79

[18] V.C. Chen, L. Hao, Time-Frequency Transforms for Radar Imaging and Signal Analysis, AretchHouse, Inc., Boston, 2002.

81

80 82

[19] L. Cohen, Time-frequency distributions-a review, Proc. IEEE 77 (7) (1989) 941–981.

83

[20] B. Boashash, P. O’Shea, Polynomial Wigner–Ville distributions and their relationship to time-varying higher-order spectra, IEEE Trans. Signal Process. 42 (1) (1994) 216–220.

85

[21] M. Morelande, B. Senadji, B. Boashash, Complex-lag polynomial Wigner–Ville distribution, in: IEEE TENCON-Speech and Image Technologies for Computing and Telecommunications, vol. 8(5), 1997, pp. 43–46.

33 34

72

[16] J.C. Wood, D.T. Barry, Linear signal synthesis using the Radon-Wigner transform, IEEE Trans. Signal Process. 42 (8) (1994) 2105–2111.



where S = S n+m+1 S n+m−1 S n∗+2m . After some tedious but straightforward derivations, we can obtain

30

[13] M. Spigai, C. Tison, J. Souyris, Time-frequency analysis in high-resolution SAR imagery, IEEE Trans. Geosci. Remote Sens. 49 (7) (2011) 2699–2711.

75

(63)

2 b60 N 4

+

70

2 b30 N 3

24 26

69

[12] F. Hlawatsch, G.F. Boudreaux-Bartels, Linear and quadratic time-frequency signal representations, IEEE Trans. Signal Process. 9 (2) (1992) 21–67.

[15] E. Sejdic, K.A. Lowry, J. Roche, M.S. Redfern, J.S. Brach, A comprehensive assessment of gait accelerometry signals in time, frequency and time-frequency domains, IEEE Trans. Neural Syst. Rehabil. Eng. 22 (3) (2014) 603–612.

∂δ f N (0 ) ≈ −2 N Im ∂

25

68

73

n =1

6

67

[14] S. Peleg, B. Friedlander, The discrete polynomial phase transform, IEEE Trans. Signal Process. (1995) 1901–1914.

Substituting (60)–(64) into (31) and (32) yields 2

[11] O. Akay, G.F. Boudreaux-Bartels, Fractional convolution and correlation via operator methods and an application to detection of LFM signals, IEEE Trans. Signal Process. 49 (5) (2001) 979–993.

  exp − j δa2 2n2 exp {− j 0 n}

N 

× exp {− j 0 n} −

18 20

2δ a 2

∂δ gN (0 ) ≈ − j ∂ 0

17 19

(1 )

Wn

−j

10

13

N  n =1

9

12

(62)

11

84 86 87 88 89

[22] B. Boashash, P. O’Shea, Design of higher order polynomial Wigner–Ville distributions, IEEE Trans. Signal Process. 47 (9) (1999) 2608–2611.

90

[23] P.H. Huang, Fast SAR imaging method for ground moving target using so-WVD transform, IEEE Trans. Geosci. Remote Sens. 54 (4) (2016) 320–335.

92

91

[24] S. Peleg, B. Friedlander, Multicomponent signal analysis using Polynomial-Phase Transform, IEEE Trans. Aerosp. Electron. Syst. 32 (1) (Jan. 1996) 378–387.

93

[25] M.Z. Ikram, K. Abed-Meraim, Y. Hua, Fast quadratic phase transform for estimating the parameters of multicomponent chirp signals, Digit. Signal Process. (1997).

95

[26] P. Wang, J. Yang, Multicomponent chirp signals analysis using product cubic phase function, Digit. Signal Process. (2006).

94 96 97 98

[27] P. Wang, H. Li, I. Djurovic, B. Himed, Integrated cubic phase function for LFM signal analysis, IEEE Trans. Aerosp. Electron. Syst. 46 (3) (July 2010) 963–977.

99

[28] P. O’Shea, A new technique for instantaneous frequency rate estimation, IEEE Signal Process. Lett. 8 (2002) 251–252.

101

100 102

[29] Shengli Wang, Shiguo Li, jin lin Ni, guang yi Zhang, A new transform Matched Fourier transform, Acta Electron. Sin. 29 (3) (2001) 403–405.

103

[30] B. Porat, B. Friedlander, Asymptotic statistical analysis of the higher order ambiguity function for parameter estimation of polynomial-phase signals, IEEE Trans. Inf. Theory 42 (3) (1996) 995–1001.

105

[31] S. Peleg, B. Friedlander, The discrete polynomial-phase transform, IEEE Trans. Signal Process. 43 (8) (1995) 1901–1914.

104 106 107 108

[32] J. Li, P. Stoica, Efficient mixed spectrum estimation with applications to target feature extraction, IEEE Trans. Signal Process. 44 (2) (1996) 281–295.

109

[33] S.Q. Zhu, A new method for radar high-speed maneuvering weak target detection and imaging, IEEE Trans. Geosci. Remote Sens. Lett. 11 (7) (2014) 1175–1179.

111

[34] 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. [35] Y. Jiang, S. Sun, Three-dimensional aircraft isar imaging based on shipborne radar, IEEE Trans. Aerosp. Electron. Syst. 52 (5) (2016) 2504–2518. [36] X. Bai, R. Tao, Z. Wang, Y. Wang, ISAR imaging of a ship target based on parameter estimation of multicomponent quadratic frequency modulated signals, IEEE Trans. Geosci. Remote Sens. 52 (2) (2014) 1418–1429.

110 112 113 114 115 116 117 118 119 120 121 122

Tong Gu was born in Anhui, China, in 1992. He received a B.S. degree in information engineering from Tongling University, Tongling, China, in 2013 and an M.S. degree in electrical engineering from Xidian University, Xi’an, China, in 2016, where he is currently working toward a Ph.D. from the National Laboratory of Radar Signal Processing, Xidian University. His research interests include synthetic apeture radar (SAR)/ISAR imaging, multiple input multiple output synthetic apeture radar (MIMO-SAR), ground moving target indication (GMTI) and three dimensional synthetic apeture radar (3D-SAR).

123 124 125 126 127 128 129 130 131 132

JID:YDSPR

12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

AID:102683 /FLA

[m5G; v1.261; Prn:6/02/2020; 15:07] P.12 (1-12)

T. Gu et al. / Digital Signal Processing ••• (••••) ••••••

Guisheng Liao (M’96-SM’16) was born in Guangxi, China in 1963. He received the B.S. degree from Guangxi University, Guangxi, China, and the M.S. and Ph.D. degrees from Xidian University, Xian, China, in 1985, 1990, and 1992, respectively. He is currently a Yangtze River Scholars Distinguished Professor at the National Laboratory of Radar Signal Processing and serves as Dean of the School of Electronic Engineering, in Xidian University. Since 2009, he has been the evaluation expert for the international cooperation project of Ministry of Science and Technology in China. Since 2007, he has been the lead of Yangtze River Scholars Innovative Team and devoted in advanced techniques in signal and information processing. Since 2006, he has served as the panelists for the medium and long term development plan in high-resolution and remote sensing systems. From 1999 to 2000, he has been a Senior Visiting Scholar in the Chinese University of Hong Kong, Hong Kong. He is the author and coauthor of several books and of more than 200 publications. His research interests include array signal processing, spacetime adaptive processing, radar waveform design, and airborne/space surveillance and warning radar systems.

20 21 22 23 24 25 26 27 28 29 30

Yachao Li was born in Jiangxi, China, in May, 1981. He received M.S and Ph.D. degrees in electrical engineering from Xidian University, Xi’an, China, in 2005 and 2008, repectively. He is currently a Professor with Xidian University. His current research interests include synthetic apeture radar (SAR)/ISAR imaging, missile-borne SAR imaging, ground moving target indication (GMTI), matching and orientation of SAR image, real-time signal processing based on FPGA and DSP technology, and distributed radar.

Yifan Guo was born in Henan, China, in 1990. He received a B.S. degree in information engineering from North China University of Water Resources and Electric Power, Zhengzhou, China, in 2012 and an M.S. degree in electrical engineering from Xidian University, Xi’an, China, in 2015, where he is currently working toward a Ph.D. from the National Laboratory of Radar Signal Processing, Xidian University. His research interests include bistatic multiple input multiple output radar (Bi-MIMO), synthetic aperture radar (SAR), ground moving target indication (GMTI) and space time adaptive processing (STAP).

67 68 69 70 71 72 73 74 75 76 77 78

Yan Huang (M’19) received the B.S. degree in electrical engineering, and the Ph. D. degree in signal and information processing, both from Xidian University, Xi’an, China, in 2013 and 2018, respectively. He was studying as a visiting Ph.D. student in Electrical and Computer Engineering department at University of Florida from Sep. 2016 to July 2017, and in Electrical and Systems Engineering department at the Washington University in St. Louis from July 2017 to Aug. 2018. He is currently an assistant professor at the State Key Laboratory of Millimeter Waves, Southeast University. His research interests include machine learning, synthetic aperture radar, image processing, remote sensing

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65

131

66

132