Quadratic Frequency Modulation Signals Parameter Estimation Based on Two-Dimensional Product Modified Parameterized Chirp Rate-Quadratic Chirp Rate Distribution

In an inverse synthetic aperture radar (ISAR) imaging system for targets with complex motion, the azimuth echo signals of the target are always modeled as multicomponent quadratic frequency modulation (QFM) signals. The chirp rate (CR) and quadratic chirp rate (QCR) estimation of QFM signals is very important to solve the ISAR image defocus problem. For multicomponent QFM (multi-QFM) signals, the conventional QR and QCR estimation algorithms suffer from the cross-term and poor anti-noise ability. This paper proposes a novel estimation algorithm called a two-dimensional product modified parameterized chirp rate-quadratic chirp rate distribution (2D-PMPCRD) for QFM signals parameter estimation. The 2D-PMPCRD employs a multi-scale parametric symmetric self-correlation function and modified nonuniform fast Fourier transform-Fast Fourier transform to transform the signals into the chirp rate-quadratic chirp rate (CR-QCR) domains. It can greatly suppress the cross-terms while strengthening the auto-terms by multiplying different CR-QCR domains with different scale factors. Compared with high order ambiguity function-integrated cubic phase function and modified Lv’s distribution, the simulation results verify that the 2D-PMPCRD acquires higher anti-noise performance and obtains better cross-terms suppression performance for multi-QFM signals with reasonable computation cost.


Introduction
The high-resolution inverse synthetic aperture radar (ISAR) imaging has been widely used in the field of civil and military in the past few decades. The ISAR image processing and relative algorithms are described in [1,2]. The primary procedures for ISAR imaging algorithm are the range alignment [3][4][5][6] and the phase adjustment after range compression. For targets with smooth motions, conventional range-Doppler (RD) method is effective in reconstructing ISAR images. For targets with the slow maneuvering, the azimuth echo signals are always modeled as the linear frequency modulation (LFM) signal. The RD method can not acquire a high quality image due to the Doppler frequency shift caused by the chirp rate. Therefore, many algorithms for LFM signal parameters estimation have been proposed, such as Wigner-Ville distribution (WVD) [7], fractional Fourier transform (FrFT) [8], cubic phase function (CPF) [9], Lv's distribution (LVD) [10] and parameterized centroid frequency-chirp rate distribution (PCFCRD) [11]. However, the azimuth echoes of targets with complex motion should be modeled as the multicomponent quadratic frequency

The 2D-PMPCRD Method
The geometry for ISAR imaging of targets with complex motions used here is based on the model in [22]. The algorithms of the range alignment and the phase adjustment are described in [6][7][8][9]23,24]. In this paper, we mainly study the processing of the Doppler frequency shift. The ISAR usually transmits the LFM signals. According to [25,26], we can find that the azimuth echoes take the multi-QFM signals model: where t (− T s * N 2 ≤ t ≤ T s * N 2 ) denotes the slow time. N and T s denote the signal length and sampling interval. P is the number of QFM signal components. A i , φ 1,i , φ 2,i , and φ 3,i denote the amplitude, the centroid frequency, the CR and the QCR of the ith QFM signal, respectively. z(t) denotes the white Gaussian noise, which has a mean of zero and a variance of δ 2 .

The Principle of 2D-PMPCRD
In order to make the algorithm easily understandable, a noise-free QFM signal is presented as where A, φ 1 , and φ 2 , φ 3 denote the amplitude, the centroid frequency, the CR and QCR of the mono-QFM signal, respectively. Based on the analysis in [11,21], a novel PSSAF with multi-scaled factors for QFM signal parameters estimation is proposed. It can be defined as where * denotes the complex conjugation, and τ is the lag-time variable. D = {d i }, 1 ≤ i ≤ L denotes the set of multi-scale factors. According to [21], the scaled factor d i is introduced to complete the order reduction, which can benefit the anti-noise performance. h denotes the constant delay. It is introduced to improve the algorithm performance in low SNR [10,11]. Thus, the selection of d i and h is important for the proposed algorithm. The selection criterion of d i and h will be discussed in Section 3.2.
A new time-quadratic chirp rate analysis technique can be obtained by performing Fourier transform on R s (t, τ, D) with respect to τ: By substituting Equation (2) into Equation (4), we can obtain where C(D) denotes Equation (5) is similar to the time-chirp rate (TCR) analysis in [11]. Thus, it has a similar characteristic to the TCR technique and it concentrates along ω = 2φ 3 D. However, TQC will not concentrate along the t-axis and it will cause poor cross-terms suppression and poor anti-noise performance. In order to make Equation (5) concentrate along the t-axis, the matrix H(t, ω) = exp(−j2πωt 2 ) is firstly introduced to dechirp TQC. After the dechirp processing, we perform the Fourier transform for the dechirped TQC with respect to the t-axis. The modified parameterized chirp rate-quadratic chirp rate distribution (MPCRD) is obtained: where G p , ω and r denote the amplitudes, quadratic chirp rate domain with respect to τ and chirp rate domain with respect to t, respectively. From the analysis of Equation (7), it is easy to see that it coherently accumulates at (r, ω) = (4φ 2 D, 2φ 3 D) on chirp rate and quadratic chirp rate (CR-QCR) domains. QFM signal will be transformed into multi-peaks on CR-QCR domains. The peaks are related to the scale factors. We set d 1 as the reference scale factor, and the related peak as the reference peak. The relationships between peaks can be presented as where P 1,i is the zoom factors, and ω 1 and r 1 are the coordinates of reference peak. ω i and r i are the coordinates of peaks, which are related to scale factors d i . L denotes the number of the scale factors. From Equation (8), we can see that the peaks are related to the scale factors, which are selected by selection criteria. The selection criteria will be discussed in Section 3.2. According to the relationships between peaks and scale factors, the 2D-PMPCRD is proposed, which can be presented as

Implementation of the 2D-PMPCRD
The key steps of the proposed algorithm are the Fourier transform with respect to τ and t. In this section, we propose mNUFFT-FFT to transform QFM signals into CR-QCR domains. The mNUFFT-FFT perform modified NUFFT along the τ-axis and modified FFT along the t-axis to speed up the implementation.
The discrete form of Equation (3) can be expressed as Equation (10) based on the sampling scheme of Classen and Mecklenbrauker [27] (CM): where T s denotes the sampling interval. n = − N 2 , − N 2 + 1, ..., N−1 2 , m = − N 2 , − N 2 + 1, ..., N−1 2 . By performing the discrete Fourier transform for Equation (10) with respect to m, a new time-quadratic chirp rate technique with the optimal scale factor, which is presented as d 1 , can be represented as In Equation (11), the order of the integration variable is two. Thus, FFT can not be used to speed up the implement. However, through the analysis of [11,28,29], we can use NUFFT to speed up the Fourier transform along the m-axis without accuracy loss. The processing of the NUFFT can be described as where NUFFT m [ * ] denotes the NUFFT process with respect to the m-axis.
Then, the matrix H(n, ω) = exp(−j2πωn 2 ) is used to dechirp TQC. After that, FFT is used to speed up the Fourier transform implement with respect to n-axis. The processing can be described as From Equation (9), we can see that 2D-PMPCRD needs to transform the QFM signal into the same position. However, the auto-terms will be transformed into different positions with different scale factors. Thus, this paper proposes mNUFFT-FFT to speed up the implement and transform the auto-terms into the same position.
In the mNUFFT-FFT, in order to transform the QFM signal into the same position along the m axis, the NUFFT with respect to m axis can be rewritten as follows: where (14) is called mNUFFT, the auto-terms with different scale factors will be transformed into the same position along the m-axis after the mNUFFT. We can use mFFT with respect to the n-axis to transform the auto-terms with different scale factors into the same position along the n-axis after the dechirp processing. The processing of mFFT can be presented as where Finally, the discrete 2D-PMPCRD can be obtained where L is the number of scale factors.

Cross-Term Suppression Performance Analysis
In this section, the cross-term suppression of 2D-PMPCRD and MLVD are analyzed. The parameter selection criterion is also discussed in this section.

Cross-Term Suppression Comparison
In this section, the cross-term suppression for 2D-PMPCRD and MLVD is discussed under multi-QFM signals. In order to explain the cross-term problem, two QFM signals are considered firstly. It can be presented as In order to make the proposed algorithm more understandable, Equation (3) that denotes PSSAF of 2D-PMPCRD can be rewritten as where Submitting Equations (17) into (19), we can obtain where R auto (t, D) are the auto-terms, R cross (t, D) are the cross-terms. The proposed algorithm with one scale factor can be described as Equation (18) and the PSSAF in [11] have similar forms. Thus, Equation (18) follows the characteristic of PCFCRD [11]. If Equation (23) takes the LFM signal model, the cross-terms will form spurious peaks, which will influence the detection of the auto-terms. Equation (23) takes the LFM form when Through the analysis of Equation (23), two cases will be proposed to explain the cross-terms influences. They are listed as follows: In these two cases, the spurious peaks will appear on CR-QCR domains for the MPCRD. However, 2D-PMPCRD can greatly suppress the spurious peaks by introducing the multi-scale factor PASSF and mNUFFT-FFT. We will explain how the proposed algorithm suppresses the spurious peaks as follows.

Case One
When ∆φ 1 = 0, ∆φ 2 = 0, ∆φ 3 = 0, the R cross (t, τ, d i ) can be presented as follows: From Equation (24), we can see that the R cross takes two LFM signals. Therefore, the cross-terms will form two spurious peaks. The spurious peaks coordinates are presented as From Equation (22), it is easily seen that the auto-terms take two LFM signals form, which means that the auto-terms will form two true peaks. The coordinates of true peaks are presented as From the analysis of the coordinates of the true peaks and the spurious peaks, we can find that the coordinates of true peaks are proportional to the scale factors. It can be presented as However, the spurious peaks coordinates do not satisfy Equation (29). From Section 2.2, we can see that mNUFFT-FFT can be used to transform the auto-terms into the same position and the cross-terms into different positions on the CR-QCR domains. After performing the mNUFFT-FFT, we can enhance the true peaks and suppress the spurious peaks by multiplying all the CR-QCR domains. In order to prove the cross-term suppression performance, we consider the noiseless and noisy condition. Example 1. Two QFM signals P 1 , P 2 are considered to show the good cross-term suppression performance of the 2D-PMPCRD compared with MLVD in this example. The sampling frequency f s is set to 256 Hz. The signal length N and effective signal length N e are set as 624 and 512. The effective signal length N e is the length of the signal after the autocorrelation. Signal parameters are set as follows: A 1 = 1, φ 1,1 = 10, φ 2,1 = 10, The constant delay is set as 2 s. In this example, two scale factors are set as d 1 = 56, d 2 = 64 for 2D-PMPCRD, which are selected by parameter selection criterion and will be discussed in Section 3.2. Figure 1 gives simulation results of 2D-PMPCRD. Figure 2 gives the simulation results of MLVD with the scale factor equal to d 1 = 56.
In Figures 1 and 2 , Cr1, Cr2 and Au1, Au2 denote the spurious peaks and the true peaks. The CR and QCR axis denote the chirp rate and quadratic chirp rate of QFM signals. From Figure 2, we can see that there are four peaks on the CR-QCR domain by using MLVD. It is not hard to see that the spurious peaks are high enough to influence the detection of the true peaks. It is a great challenge for MLVD. Thus, the MLVD can not be used to estimate parameters of the QFM signals when multi-QFM signals have the same QCR. The 2D-PMPCRD overcomes the identifiability problem. From Figure 1, we can see that the true peaks are strengthened and the spurious peaks are greatly weakened after the 2D-PMPCRD operator. The simulation results have proven the theoretical analysis.

Example 2.
In this example, we analyze the cross-terms suppression under noisy multi-QFM signals. The parameters of QFM signals are set the same as Example 1. The input SNR is set as 0 dB. Figure 3 gives simulation results of 2D-PMPCRD. Figure 4 gives the simulation results of MLVD. In Figures 3 and 4, we can see that the true peaks are much higher than spurious peaks for the 2D-PMPCRD under noisy condition. However, the spurious peaks for the MLVD are high enough to influence the detection of the true peaks. Thus, the 2D-PMPCRD has better cross-term suppression ability than MLVD.

Case Two
When ∆φ 1 = 0, ∆φ 2 = 0, ∆φ 3 = 0, the R cross (t, τ, d i ) can be presented as follows: It is not hard to see that the coordinates of spurious peaks can be described as follows: Because φ 3,1 is equal to φ 3,2 and φ 2,1 is equal to φ 2,2 , we see that the auto-term will form a higher peak, which is different from case one. The position of the true peak is (4φ 2: d i , 2φ 3: d i ), where φ 2: = φ 2,1 = φ 2,2 and φ 3: = φ 3,1 = φ 3,2 , and it meets Equation (29). However, the coordinates of the spurious peaks do not meet Equation (29). Thus, the process to enhance the true peak and suppress the spurious peaks are the same as case one. In this case, we also consider the noiseless and noisy condition. Example 3. Two QFM signals P 3 , P 4 are considered to show the good cross-term suppression performance of the 2D-PMPCRD compared with MLVD in this example. The sampling frequency f s is set as 256 Hz. The signal length N and effective signal length N e are set as 624 and 512. Signal parameters are set as follows: A 3 = 1, φ 1,3 = 10, φ 2,3 = 10, φ 3,3 = 4 for P3; A 4 = 1, φ 1,4 = 40, φ 2,4 = 10, φ 3,4 = 4 for P4. The constant delay is set as 2 s. In this example, two scale factors, which are d 1 = 56 and d 2 = 64, are selected for 2D-PMPCRD. Figure 5 gives simulation results of 2D-PMPCRD. Figure 6 gives the simulation results of MLVD with the scale factor equal to d 1 = 56. In Figures 5 and 6, Cr3, Cr4 and Au3, Au4 denote the spurious peaks and the true peaks. From Figure 6, we can see that two spurious peaks appear on the CR-QCR domain, which will influence the detection of the true peak. Thus, MLVD can not be used under this condition. From the analysis of Equation (22), we can see that the true peaks will appear on the same position on the CR-QCR domain when the multi-QFM signals have the same QCR and CR. From Figure 6, we can see that the simulation results have proven the theoretical analysis.

Example 4.
In this example, we analyze the cross-terms suppression under noisy multi-QFM signals. The parameters of QFM signals are set the same as Example 3. The input SNR is set as 0 dB. Figure 7 gives simulation results of 2D-PMPCRD. Figure 8 gives the simulation results of MLVD.
In Figures 7 and 8, we can see that the simulation results are the same as Example 3. From Figures 5 and 7, we can see that the difference between Figures 5 and 7 is not obvious. In fact, there are differences between Figures 5 and 7 if we calculate the peak to mean background ratio (PMR), which is the ratio of the max value of the modulus of Q 2D-PMPCRD (r, w, D) to the mean value of the modulus of Q 2D-PMPCRD (r, w, D). The PMRs of Figures 5 and 7 are about 42.07 dB and 33.24 dB, respectively. We can see that the 2D-PMPCRD also has good cross-term suppression performance under noisy environment. From the above simulation analysis, it is easy to see that the 2D-PMPCRD has better cross-term suppression ability than MLVD and can solve the identifiability problems of MLVD.

Parameter Selection Criterion
In this section, the selection criterion of the constant delay h, the scale factor d i and the number of scale factors L are analyzed.

Selection of Constant Delay
According to the analysis in [10,11], we know that the constant delay h is related to signal length and anti-noise performance. It can benefit the anti-noise ability if we set a suitable value. Thus, in this section, we discuss the select of the constant delay h.
With the analysis in [11], we can find that the anti-noise ability achieves its best when h ≥ N e * T s . Since the effective length N e of the signal is finite in realistic application, h should not be increased without limit. Here, considering the anti-noise performance and the effective length of the signal, the constant delay h is set as In realistic applications, we have much data that are stored in the storage. Therefore, we can use them to satisfy the constant delay [10].

Selection of the Scales Factors and Number of Scale Factors
Through the analysis in [21], we can see that the scale factors can not only reduce the signal order, but also benefit the anti-noise performance. The selection criteria of scale factor d i is to obtain the smallest MSE of QCR and CR according to [15]. Through the analysis in [15], the optimal d i , which is denoted as d opt , can be presented as where * denotes the round up function. N denotes the signal length. According to [17], we can see that the optimal choice of L lags are d 1 = d 2 = · · · = d L−1 = d opt . Thus, the optimal scale factors of 2D-PMPCRD are d 1 = d 2 = · · · = d L−1 = 0.089N . However, in Section 3.1, we can see that the cross-terms form spurious peaks if the scale factors are set as the same. The spurious peaks influence the detection of the true peaks. To solve this problem, the simple method is to use different scale factors. Thus, the scale factors are usually selected as different ones in realistic application. Considering the accuracy of CR and QCR estimation, the sub-optimal scale factors should be selected around the optimal one [17].
From the above analysis, we can see that the selection of the number of scale factors L is also very important for the performance of 2D-PMPCRD. According to [17], estimation performance will be improved as the number of scale factors L increases. However, we can see that the anti-noise performance is improved slightly when L ≥ 2. In this paper, from Section 3.1, we can see that the proposed algorithm can acquire good cross-term suppression performance when L is set to 2. Furthermore, from Section 2, we can see that the computation cost will increase rapidly as the number of scale factors L increases. Thus, considering the computation cost and performance, L is usually selected as 2.

Anti-Noise Performance and Computation Cost Analysis
In this section, anti-noise performance and computation will be analyzed for 2D-PMPCRD. Both the mono-component QFM signals and multi-QFM signals are considered to prove the good anti-noise performance of 2D-PMPCRD.

Anti-Noise Performance Analysis
In this section, we propose Examples 5 and 6 to illustrate the anti-noise performance of 2D-PMPCRD.
Example 5. One QFM signal that denotes P 5 is considered. The sampling frequency f s is set as 256 Hz. The signal length N and effective signal length N e are set as 624 and 512. Parameters for P 5 are set as follows: A 5 = 1, φ 1,5 = 10, φ 2,5 = 10, φ 3,5 = 4. The set of scale factors are set to d 1 = 56, d 2 = 64. The constant delay h is set to 2 s. The test input SNRs in figures are SNR = [−8:1:−2]. In this example, 100 simulations are run for each SNR value. The input-output SNR [21] is used as a measurement and the result is shown in Figure 5a. The observed mean square error (MSE) and the Cramer-Rao Lower Bound (CRLB) of quadratic chirp rate and chirp rate in different SNRs are shown in Figure 5b,c, respectively.
From Figure 9, we can see that the anti-noise performance of HAF-ICPF is the worst. The anti-noise performance of MLVD is better than the HAF-ICPF. However, the threshold SNR of MLVD is still high. The anti-noise performance of 2D-PMPCRD is the best.
The HAF-ICPF uses HAF to reduce the coefficient order of QFM signals. However, HAF is sensitive to noise and suffers from the error propagation. Thus, from Figure 9a, we can see that the threshold SNR of HAF-ICPF is only −2 dB. The evidence of the simulation results can also be found in [15,20]. The MLVD introduces a novel PSSAF and the keystone transform to transform the QFM signals into CR-QCR domain. The MLVD can not only improve the anti-noise performance, but also overcome the error propagation by estimating the CR and QCR together on the CR-QCR domain. However, the MLVD only uses one scale factor, which is not suitable in low SNR. From Figure 9a, we can see that the threshold SNR of MLVD is −4 dB. The evidence of the simulation result can also be found in [21]. The 2D-PMPCRD uses a constant delay h and multi-scale factors D to improve the performance in low SNR. In addition, the multiplying operation of different CR-QCR domains can also improve the anti-noise performance. Therefore, the anti-noise performance of 2D-PMPCRD is the best. The threshold SNR of 2D-PMPCRD is −6 dB.
In Figure 9b,c, the MSEs of CR and QCR estimations are plotted as functions of SNR. From the Figure 9b,c, we can see that, for the 2D-PMPCRD, the MSEs of CR and QCR estimations increase abruptly at −6 dB. This is because, if the SNR is lower than the threshold SNR, there is high probability that the maximum peak on the CR-QCR domain are some random peaks, unrelated to (φ 2,5 , φ 3,5 ) of the QFM signal. This will cause the parameter estimation algorithm to fail. Therefore, the MSEs of CR and QCR increase abruptly at the threshold SNR. Example 6. Two QFM signals P 6 , P 7 are considered to show the good anti-noise performance of the 2D-PMPCRD compared with MLVD under noisy multi-QFM signals. The sampling frequency f s is set as 256 Hz. The signal length N and the effective signal length N e are set as 624 and 512. Signal parameters are set as follows: A 6 = 1, φ 1,6 = 10, φ 2,6 = 10, φ 3,6 = 4 for P 6 ; A 6 = 1, φ 1,7 = 40, φ 2,7 = 50, φ 3,7 = 12 for P 7 . The constant delay is set as 2 s. In this example, two scale factors, which are d 1 = 56 and d 2 = 64, are selected for 2D-PMPCRD. SNR is set as −3 dB. Figure 10a gives simulation results of 2D-PMPCRD. Figure 10b gives the simulation results of MLVD with the scale factor equal to d 1 = 56.
From Figure 10, we can see that 2D-PMPCRD can distinguish these two QFM signals. The true peaks of MLVD are immersed in the noise. Therefore, 2D-PMPCRD has better anti-noise performance than MLVD when noisy multi-QFM signals are considered.

Computation Cost Analysis
Computation cost is an important factor that should be considered in realistic application. In this section, the computation cost of 2D-PMPCRD will be discussed and the computation cost of the comparison algorithms HAF-ICPF and MLVD will also be analyzed.
For HAF-ICPF, the main procedures are HAF and ICPF, and their computation costs are (O(N)) [14] and (O(N 3 )) [30,31], respectively. Thus, the computation cost of HAF-ICPF is (O(N 3 )). For MLVD, the main procedures are PSSAF and keystone transform. The computation of PSSAF is (O(N 2 )) [21]. The keystone transform can be divided into FFT and IFFT [21]. Thus, the computation cost of keystone transform is similar to FFT. As we all know, the computation of FFT is (O (N 2 log 2 N)). Thus, the computation cost of MLVD is (O(N 2 log 2 N)) [21]. The main procedures of 2D-PMPCRD are PSSAF, dechirp operation, mNUFFT-FFT and L − 1 products of MPCRD. The computation cost of the PSSAF of 2D-PMPCRD is similar to MLVD, thus its computation cost is (O(N 2 )). The computation cost of dechirp operation is (O(N 2 )) [11]. From the analysis of Section 3.1, we know that the mNUFFT-FFT can be divided into NUFFT and FFT. Thus, the computation cost of the mNUFFT-FFT is the same as NUFFT and FFT, whose computation costs are (O(N 2 log 2 N)) [11,28,29]. In realistic application, L is much smaller than the signal effective length (e.g., 2). Thus, the computation of 2D-PMPCRD is (O(N 2 log 2 N)). The computation costs of these three algorithms are listed in Table 1. Table 1. Computation cost of algorithms.

Estimation Algorithm Computation Cost
From Table 1, we can see that the computational cost of HAF-ICPF is the worst. The computational cost of 2D-PMPCRD is equivalent to MLVD. Therefore, 2D-PMPCRD can obtain better cross-terms suppression performance and anti-noise performance with reasonable computational cost. Compared with HAF-ICPF and MLVD, the proposed algorithm is more suitable for the real-time ISAR imaging system.

Conclusions
In this paper, we propose the 2D-PMPCRD to estimate CR and QCR of the QFM signals by introducing a novel multi-scale PSSAF and the mNUFFT-FFT. The PSSAF utilizes the multi-scale factors and the constant delay to improve the anti-noise performance. The mNUFFT-FFT transforms