Analysis of Maneuvering Targets with Complex Motions by Two-Dimensional Product Modified Lv’s Distribution for Quadratic Frequency Modulation Signals

For targets with complex motion, such as ships fluctuating with oceanic waves and high maneuvering airplanes, azimuth echo signals can be modeled as multicomponent quadratic frequency modulation (QFM) signals after migration compensation and phase adjustment. For the QFM signal model, the chirp rate (CR) and the quadratic chirp rate (QCR) are two important physical quantities, which need to be estimated. For multicomponent QFM signals, the cross terms create a challenge for detection, which needs to be addressed. In this paper, by employing a novel multi-scale parametric symmetric self-correlation function (PSSF) and modified scaled Fourier transform (mSFT), an effective parameter estimation algorithm is proposed—referred to as the Two-Dimensional product modified Lv’s distribution (2D-PMLVD)—for QFM signals. The 2D-PMLVD is simple and can be easily implemented by using fast Fourier transform (FFT) and complex multiplication. These measures are analyzed in the paper, including the principle, the cross term, anti-noise performance, and computational complexity. Compared to the other three representative methods, the 2D-PMLVD can achieve better anti-noise performance. The 2D-PMLVD, which is free of searching and has no identifiability problems, is more suitable for multicomponent situations. Through several simulations and analyses, the effectiveness of the proposed estimation algorithm is verified.


Introduction
In the past decades, radar signal processing has been a very important research field [1][2][3], especially estimating radar parameters [4][5][6]. The parameters of radar signals mainly include frequency [7], phase [8], and direction-of-arrival [9][10][11]. In this paper, we only focus on research into the frequency of the quadratic frequency modulation (QFM) signal [2,4]. Recently, many researchers have paid much attention to presenting effective inverse synthetic aperture radar (ISAR) imaging algorithms, which are important for high-resolution ISAR imaging and recognition of moving targets [12,13]. In general, the primary steps for ISAR imaging are range alignment [14,15] and phase adjustment [16,17], and translational-induced phase error correction. Then, the conventional range-Doppler (RD) method is utilized to reconstruct a well-focused ISAR image. However, in practical application, as the target's motions are not smooth, the conventional RD method is not valid. For slow maneuvering targets, the azimuth echoes can be modeled as a linear frequency modulated (LFM) signal [18][19][20][21], and as the Doppler frequency shifts, we cannot obtain high-resolution imaging using the conventional RD. The chirp rate is identified as the cause of the target image defocus. Recently,

2D-PMLVD with QFM Signal
In this section, the model of the QFM signals and the principle of the 2D-PMLVD are introduced, respectively.
For ISAR imaging of targets with complex motions, the geometry used here is based on the model in [31,32]. In this paper, range alignment and translational induced phase error correction, which can be referred to in [14][15][16][17], will not be discussed in detail. We will only focus on the processing of the Doppler frequency shift. Suppose that the radar transmits the LFM signal. According to the analyses in [31,32], after the range alignment and the translational induced phase error correction, azimuth echoes of the lth range cell take the form of multicomponent QFM signals where P is the phase of the pth component, t denotes the slow time and its sampling frequency is the same as the pulse repetition frequency (PRF), z(t) denotes the additive complex white Gaussian noise with a variance of δ 2 , φ 1,p , φ 2,p , φ 3,p denote the centroid frequency, the chirp rate and the quadratic chirp rate, respectively.

2D-PMLVD with Mono-QFM Signal
In order to facilitate the understanding of the algorithm, we first consider a noise-free QFM signal. It is represented as where φ(t) is the phase function φ 1 , φ 2 and φ 3 denote the centroid frequency, the chirp rate and the quadratic chirp rate, respectively. According to analyses in [21,33], a novel multi-scale PSSF is defined as where a denotes a constant time-delay, τ denotes the lag time variable, and * denotes the complex conjugation.
is the set of scale factors, and the selection criteria will be presented in Section 2.3.1. According to [33,34], the scale factor d i and the corresponding lag-time τ are utilized to complete the order reduction. Besides that, d i can reduce the number of self-correlations, which benefits the anti-noise performance.
Substituting (2) into (3) yields It is easily seen from Equation (4) that the time variable t and lag variable τ couple with each other in exponential phase terms. For every scale factor, Equation (4) suffers the same influence of the coupling; we only need to analyze Equation (4) with one scale factor to explain the influence. Thus, Equation (4) is represented as The form of Equation (5) is the same as the self-correlation function defined for the MLVD [33] and LVD [21]. According to [21,33], we have known that the coupling influence energy accumulation.
In order to analyze the influence of the coupling for the energy accumulation of the proposed algorithm, we use the energy accumulation method presented in [21,33], referred to as EA. The EA method is also the basic energy accumulation method used for the proposed algorithm. The EA method is divided into two steps. Firstly, perform FFT on Equation (5) with respect to t. Secondly, perform FFT on the first step with respect to τ. After the first step, the signal energy will peak along the inclined line f t = 2φ 3 d i (τ + a). It is easily seen that the signal energy accumulation cannot be accomplished by the second step. Thus, in order to effectively accumulate the energy of the signal, the coupling needs to be removed. Here, the idea of keystone transformation is borrowed to rescale the time axis, which is defined as , τ where a is a constant time-delay and h is a scaling factor different from d i . The choice of parameters a and h will be presented in Section 2.3.1. Φ(t, τ) is the phase function with respect to (t, τ). K is the scaling operator. t n is denotes scaled time.
Applying the scaling operator K to (5), we have It is easily seen from (7) that the coupling has been removed by introducing the new time variable t n . Then, by applying the first step of EA to Equation (7), the signal energy becomes a beeline f t n = 2φ 3 d i /h. Thus, the signal energy can be accumulated by the second step. The process can be expressed as where δ(·) denotes the Dirac delta function, FFT(·) is the Fast Fourier transform operator. Equation (8) is similar to MLVD [33]. It can be easily seen from Equation (8) that a sole peak appears on the 2-D frequency domain f τ − f t n . The coordinate of the peak is f τ,i = 2φ 2 d i , f t n ,i = 2φ 3 d i /h . According to the coordinate of the peak, the parameters φ 1 , φ 2 , and φ 3 can be estimated. However when the EA method is applied to the multi-scale PSSF, the process is represented as It is easily seen from Equation (9) that the mono-component signal will be translated into multi-peaks on f τ − f t n . The number of peaks is the same as the number of scale factors. Suppose that d 1 is the reference factor, and the peak which is correspond to d 1 is the reference peak. The relationships between reference peak and other peaks are represented as where P 1,i is zoom factor. f τ,1 and f t n ,1 are the coordinates of the reference peak on the 2-D frequency domain. f τ,i and f t n ,i are the coordinates of the peak which corresponds to d i . It can be easily seen from Equation (10) that scale factor influences the location of the peak. Due to the values of scale factors have been selected by the selection criteria, which will be presented in Section 2.3.2, the relationships in Equation (10) can be easily obtained. According to these relationships, an improved energy accumulation method is presented, which is called mSFT t,i − FFT τ,i . The mSFT t,i − FFT τ,i not only inherits advantages of decoupling but also improves the performance of energy accumulation. By using the mSFT t,i − FFT τ,i , the peaks of the QFM signal are transformed into different 2-D frequency domains, but at the same position under different scale factors. After multiplying all the 2-D frequency domains, a higher peak will appear. The mSFT t,i − FFT τ,i will be presented in detail in Section 3. The processing is named 2D-PMLVD, which is represented as 2D-PMLVD not only suppresses cross terms, but also improves the anti-noise ability, which will be proved in Sections 2.2 and 4.1.
2D-PMLVD is divided into four steps. Firstly, select one of the scale factors as the reference one and calculate the relationships between reference factor and other factors by Equation (10). Secondly, transform the QFM signal into different 2-D frequency domains based on mSFT t,i − FFT τ,i . Thirdly, multiply all the 2-D frequency domain to obtain the higher peak. Fourthly, detect the peak value to estimate the signal parameters.
For example, suppose that the set of scales contains two factors which are represented as D = {d 1 , d 2 }. According to Equation (10), we get where d 1 is the reference scale. After applying mSFT t,i − FFT τ,i to the QFM signal, two peaks will appear on different 2-D frequency domains. The reference peak will appear on f τ,1 − f t n ,1 , and the other peak will appear on f τ,2 − f t n ,2 . However, due to the mSFT t,i − FFT τ,i operator, the relation between the different frequency domain is represented as In order words, the new frequency interval becomes [0, 1/P 1,2 ]. Thus, because of the transformation, the two peaks appear at the same position. Then, a higher peak is obtained by multiplying the different frequency domains. According to the analysis in [25], the proposed algorithm can be improved when the number of scale factors is increased. Then, the computational cost will increase simultaneously. Hence, the number of scale factors will not be too large. In this paper, the number is selected as 2. By this selection, the cross terms can be effectively suppressed which are shown in Figures 1 and 2.

2D-PMLVD with Multi-QFM Signals
In this part, 2D-PMLVD with multi-QFM signals is analyzed. The proposed algorithm not only improves the ability to suppress cross terms, but also improves the performance of the EA method. The EA method, which is presented in MLVD [33] and LVD [21], has been proved to be effective for suppressing cross terms. Now, model the noise-free multi-QFM signals as To formulate the cross term problem arising from multi-QFM signals, we first consider two signals (P = 2) in Equation (14), which is represented as For ease of understanding, the multi-scale PSSF can be written as Submitting (14) into (16b), we obtain where R auto (t, D) and R cross (t, D) denote auto-terms and cross terms with different scale factors, respectively. In order to describe the shortcoming of the EA method and explain the principle of the energy accumulation method in the proposed algorithm, we firstly consider one scale factor in Equation (17). It is represented as According to [21], Equation (16a) is similar to the form of the self-correlation function defined for the LVD. Thus, under multi-QFM signals, if Equation (18c) takes the form of the LFM signal, the cross terms will accumulate as an auto-term. Let LFM signals will exist in the cross terms as in the following two cases [32,35]: ∆φ 1 = 0, ∆φ 2 = 0, and ∆φ 3 = 0. 2. ∆φ 1 = 0, ∆φ 2 = 0, and ∆φ 3 = 0. In these two cases, the spurious peaks will appear on f τ − f t n , which greatly interferes with the performance of the detection. However, 2D-PMLVD overcomes the problem by introducing the multi-scale PSSF and mSFT t,i − FFT τ,i . This will be analyzed in detail as follows.

Case One
When ∆φ 1 = 0, ∆φ 2 = 0, and ∆φ 3 = 0, the R cross (t, d i ) is represented as where φ 3,: = φ 3,1 = φ 3,2 . It can be easily seen from Equation (19) that R cross (t, d i ) takes the form of two LFM signals. This means that the cross terms will accumulate two spurious peaks. The coordinates of the spurious peaks are It can be easily seen from Equation (18b) that R auto (t, d i ) takes the form of two LFM signals. The coordinates of the true peaks are For different scale factors, the coordinates of the true peaks meet Equation (10), and can be represented as However, through the analysis of Equation (20), the coordinates of the spurious peaks do not meet Equation (10). In other words, the coordinates are not proportional to the scale factor. Thus, the spurious peaks cannot be transformed into the same position by the operator mSFT t,i − FFT τ,i . After multiplying all the frequency domain, the spurious peaks will not be enhanced. Compared with the true peaks, the spurious peaks are "suppressed". The result is based on the premise that the scale factors are different. This is because if the scale factors are equivalent, the spurious peaks are in the same position, which means that the cross terms cannot be suppressed. This is proved by Example 1. The selection criterion for the scale factor is introduced in Section 2.3.2. Example 1. We consider two QFM signals P1 and P2. Signal parameters are set as follows: A 1 = 1, φ 1,1 = 89, φ 2,1 = 42, φ 3,1 = 4 for P1; and A 2 = 1, φ 1,2 = 10, φ 2,2 = 20, φ 3,2 = 4 for P2. The sampling frequency F s is set to 256 Hz. The effective signal length N is equal to 256. Both a and h are set to one. The quadratic chirp rate of P1 and P2 are equal to each other. Figure 1a shows the results of MLVD with a scale factor equal to 56. Au1 and Au2 denote the true peaks of the auto-terms. Because φ 3,2 = φ 3,1 , two spurious peaks, denoted by Cr1 and Cr2, exist in the 2-D frequency domain alongside the true peaks. It can be easily seen that the spurious peaks are even higher than the true ones in Figure 1a. This creates a serious challenge for signal detection. Figure 1b shows the results of 2D-PMLVD with two scale factors equal to 56 and 64. Because the multi-scale PSSF and the operator mSFT t n ,i − FFT τ,i were applied, the spurious peaks of the cross terms are suppressed, and the true peaks of the auto-terms are greatly strengthened. Figure 1c,d gives the contour of Figure 1a,b, respectively.

Case Two
When ∆φ 1 = 0, ∆φ 2 = 0, and ∆φ 3 = 0, R auto (t, d i ) is represented as where φ 2,: = φ 2,1 = φ 2,2 . And R cross (t, d i ) is represented as Unlike case one, R auto (t, d i ) only takes the form of one LFM signal in case two. That is because the two LFM signals are the same. It means that, by the proposed algorithm, only a higher peak of the auto-terms is obtained. The coordinates of the true peak are (2φ 2,: d i , 2φ 3,: d i ). It can be easily seen that the coordinates are proportional to the scale factor.
It can be easily seen from Equation (24) that R cross (t, d i ) also takes the form of two LFM signals. The coordinates of the two peaks are (φ 1,2 − φ 1,1 + 2φ 2,: d i , 2φ 3,: d i ) and (φ 1,1 − φ 1,2 + 2φ 2,: d i , 2φ 3,: d i ), respectively. Thus, the coordinates are not proportional to the scale factor. This is the same as in case one.
The process for suppressing the cross terms is the same as in case one. After applying 2D-PMLVD to case two, the spurious peaks can effectively be suppressed and the identifiability is improved. This is shown in Example 2.

Example 2.
We consider two QFM signals P3 and P4. Signals parameters are set as follows: A 3 = 1, φ 1,3 = 10, φ 2,3 = 20, φ 3,3 = 10 for P3; and A 4 = 1, φ 1,4 = 30, φ 2,4 = 20, φ 3,4 = 10 for P4. The sampling frequency F s is set to 256 Hz. The effective signal length N is equal to 256. Both a and h are set to one. The quadratic chirp rate and chirp rate of P3 and P4 are the same. Figure 2a shows the results of MLVD with a scale factor equal to 56. Because the quadratic chirp rate and chirp rate are the same, the auto-term is accumulated into one peak, which is represented as Au3(Au4) . Cr3 and Cr4 denote the spurious peaks. Figure 2b shows the results of 2D-PMLVD with two scale factors. By means of 2D-PMLVD, the spurious peaks of the cross terms are suppressed, and the true peak is greatly strengthened. Figure 2c,d gives the contour of Figure 2a,b, respectively.

Parameter Selection and Resolution
In this part, the selection criteria of the parameters a, h and scale factors are presented. In practical applications, the QFM signal is of finite-length and discrete. This means that the selection criteria only need to be discussed in a finite-length discrete setting.

Selection of a and h
We know that the optimal value of ah is equal to one [21,36]. Define the discrete parameter pair (a, h) as (m/ f s , h). In order to directly read the quadratic chirp rate from the 2-D frequency domain, the optimal value of h is set to one. According to the length of redundancy, the selection of the parameters follows two criteria.
Criterion 1: According to the analysis in [21], as long as the signal is long enough to guarantee at least 1 s redundancy, that is to say N all / f s − N 2D−PMLVD / f s ≥ 1, the parameters are N all and N 2D−PMLVD denote the length of the signal and the effective length of the signal, respectively. In this case, the parameters of a and h are equal to one, and the previously processed data can be used as redundancy information to perform the proposed algorithm. Because 2D-PMLVD includes a scale factor set that is different from the LVD, the effective length of 2D-PMLVD is shorter than that of LVD. It means that 2D-PMLVD is much more suited to fulfilling the condition of redundancy.
Criterion 2: However, in practical applications, the total sampling time may be very short, even less than 1 s, or the frequencies of the chirp signal may vary quickly with time even if the length of signal is long enough. The parameter selection criterion is not satisfied by the above conditions. The new condition is N all / f s − N 2D−PMLVD / f s < 1. To guarantee the condition that a is equal to one, the parameters are represented as The derivation of Equation (26) is in accordance with [21], which aims at making h as close as possible to one under the condition of ah = 1.

Selection of d i
According to the analysis in [33], we know that the scale factors not only reduce order, but also benefit anti-noise performance. Thus, their selection is very important for 2D-PMLVD. According to the analysis in [37], we know that the optimal d i value is obtained by minimizing the values for the MSE of chirp rate and quadratic chirp rate. This is represented as where • denotes the round up operator. We have seen in Section 2.2 that, in order to suppress the cross terms, the scale factors in D will be different. When we move away from the optimal choice, we will certainly obtain a sub-optimal set. As far as accuracy is concerned, it would be better to have the effective signal length be as close to the optimal length as possible, which suggests that the selection of the other scale factors in the set D should be close to d opt i . For a mono-QFM signal, the values of scale factors may be equivalent, as there is no cross term. However, in actual processing, we do not know the QFM signal is mono-component or multi-component. Thus, the scale factors in D are different.

Resolution
In the proposed algorithm, which is similar to the frequency resolution of Fourier transform, the resolution for representing the parameters of CR and QCR relate to the signal length, which makes sense, since the 1-D QFM signals are converted into the 2-D single frequency signals in the CR-QCR domain. In this paper, due to the sampling scheme of Claasen and Mecklenbräuker (CM) and the finite-length effect, the resolution for representing the parameters of CR and QCR are different. Moreover, the resolutions also relate to the scale factors. We have seen that the multi-scale PSSF can be divided into Equation (16a,b), and that Equation (16a) is the same as PSIAF, which is presented in [21]. Thus, some theorems in [21] can also be applied to this paper. Applying Equation (16b) to a QFM signal, the auto-terms can be represented as different LFM signals with different scale factors and the centroid frequency and the chirp rate of the ith LFM signal can be expressed as Suppose ∆ denotes sampling interval. According to the Appendix B in reference [21], we have known that the resolutions of centroid frequency and the chirp rate are expressed as where N i denotes the length of the ith LFM signal; and N i is represented as N i = N all − 2d i , where N all denotes the length of the QFM signal. According to Equations (28) and (29), the corresponding resolution for ϕ 2,p and ϕ 3,p is represented as It can be easily seen from Equation (30) that different scale factors correspond to different resolutions; and the resolutions are influenced by the worst resolutions in the proposed algorithm. Thus, the resolutions are influenced not only by the signal length but also by the scale factors.

Implementation
In this section, the discrete implementation of 2D-PMLVD will be discussed by using the sampling scheme of Claasen and Mecklenbräuker (CM) [38]. The discrete form of (4) expressed as The key step of the implementation is the scaling transform. Similar to the keystone transform, the operator K can be fulfilled by using DFT-IFFT or SFT-IFFT, which are analyzed in detail in [19]. According to the analysis in [19], the SFT-IFFT transformation only uses complex multiplications and FFT based on the scaling principle, and is more suitable for implementation than DFT-IFFT. If the Fourier transform is represented as then the representation of SFT is expressed as S(λ, f ) = F (r(t)) = r(t) exp(−j2πλ f t)dt (33) where λ = (2mT s + a)h/N. According to the analysis in Sections 2.1 and 2.2, the energy of the auto-terms is expected to be accumulated into the same position under different scale factors. However, the EA method, which is presented in [21,33], is only valid for mono-scale factors, and is invalid for multi-scale factors. According to the analysis of SFT, an improved energy accumulation method is proposed, called modified SFT-FFT (mSFT t,i − FFT τ,i ). The operator is divided into two parts: a performing mSFT t,i [•] operator along the t axis and a performing FFT τ,i [•] operator along the τ axis.
The mSFT t,i [•] is different from Equation (33), which is represented as where R(t, τ, d i ) denotes the PSSF when scale factor is where f 1 denotes the reference frequency domain and f ∆ i denotes the ith frequency domain, which corresponds to the ith scale factor. Compared with the interval of f 1 , the interval of f ∆ i becomes [0, 1/P 1,i ] by introducing P 1,i . The variable λ in mSFT t,i [•] is the same as in Equation (33). Thus, the coupling between t and τ can also be removed by the mSFT t,i [•] operator. Based on Equation (10), we know that P 1,i denotes the relationships between reference energy peak position and the ith energy peak position. By using the zoom factor P 1,i , the mSFT t,i [•] operator changes the frequency interval. Then, under different scale factors, the signal energy of the auto-term is accumulated at the same position on the f t n axis. In other words, under different scale factors, the signal energy of the auto-term is accumulated into the same beeline which is represented as f t n = X(τ), where X(τ) denotes the function of τ. Then, applying the FFT τ,i [•] operator along the τ axis, the signal energy of the auto-term is accumulated as a sharp peak in the 2-D frequency domain. For the ith scale factor, the number of FFT point is Num 1 /P 1,i , where Num 1 denotes the number of the FFT point for the reference scale factor. Thus, mSFT t,i − FFT τ,i is represented as Through analyses of the implementation, mSFT t,i [•] not only inherits advantages of the SFT with regard to the transformation of the time-domain, but also improves the performance in the frequency-domain. Moreover, the computation cost does not increase too much. The flowchart of mSFT t,i − FFT τ,i is shown in Figure 3.   Step1): Substitute QFM signals and D into Equation (3), and R s (t, τ, D) is obtained.
Step2): for i = 1 : L do Apply the mSFT t,i operator to R s (t, τ, d i ) along the t axis by Equation (34). Apply the FFT τ,i [•] operator to Equation (34) along the τ axis, and we can obtain the data, represented as β i . end for Step3): Multiply all β i , we can obtain different peaks, which correspond to different QFM components.
Step4): Detect the peak value to estimate the signal parameter chirp rate and quadratic chirp rate. Output: the QFM signal parameter chirp rate and quadratic chirp rate.

Analyses of Anti-Noise Performance and Computational Cost
In this section, the anti-noise performance and computational cost, which play important roles in parameters estimation [32][33][34], will be analyzed for 2D-PMLVD. HAF-ICPF [30], MLVD [33], and the algorithm presented in [34] are chosen as comparisons.

Computational Cost
The anti-noise performance of PMVLD has been analyzed in the Section 4.1. In this section, the computational cost will be analyzed for PMVLD.
For HAF-ICPF, the main implementation procedures include HAF   Table 1 lists the computational costs of these four estimation algorithms.  The input-output SNR comparison of the four algorithms is shown in Figure 4a, and the solid line denotes the result of the matched filter for P5. HAF-ICPF contains HAF and ICPF. HAF is used to reduce the order of the QFM signal, and then ICPF is used to accumulate the energy and estimate the parameter. We know that the threshold SNR of HAF is low [2]. However, according to [40], ICPF will improve the threshold SNR. Therefore, the threshold SNR of HAF-ICPF is −2dB. The evidence for this result can be obtained as described in [33,37]. The algorithm presented in [34], which employs a fourth-order autocorrelation function, is based on coherent accumulation [34]. The algorithm presented in [34] defines a novel parametric autocorrelation function to complete the order reduction and energy accumulation, which benefits anti-noise performance, although redundancy information is not utilized. Thus, in Figure 4a, the threshold SNR of the algorithm presented in [34] is −3 dB. In contrast to the algorithm presented in [34], MLVD utilizes redundancy information, and therefore the anti-noise performance is based on a better anti-noise algorithm. However, MLVD only utilizes one scale factor, which is not favorable for anti-noise performance. Thus, its threshold SNR is −4 dB, which conforms to the result in [33]. In Figure 4a, the threshold SNR of 2D-PMLVD is −5 dB, which is better than the other three algorithms. The advantages include: (1) the energy accumulation of 2D-PMLVD is coherent; (2) in the defined multi-scale PSSF, the redundancy information and multi-scale factors not only reduce the number of self-correlations, but also benefit anti-noise performance [33,37]; and (3) by the mSFT t n ,i − FFT τ,i and production, the cross terms are effectively suppressed. In addition, it is worthwhile noting that the error propagations exist in HAF-ICPF as a result of the sequential estimation of chirp rate and the quadratic chirp rate. However, this problem does not exist in the other three algorithms.
For 2D-PMLVD, the observed MSEs for the chirp rate and the quadratic chirp rate estimation are plotted in Figure 4b,c as functions of SNR. The corresponding Cramer-Rao bounds (CRBs), which can be obtained as described in [41,42], are also shown by the solid line in the Figure 4. Obviously, in Figure 4b,c, both the MSEs of the chirp rate and quadratic chirp rate estimations are close to CRB when SNR ≥ −5dB. Therefore, the threshold SNR of the 2D-PMVLD is −5 dB for MSEs, and these results conform to those in Figure 4a.

Computational Cost
The anti-noise performance of PMVLD has been analyzed in the Section 4.1. In this section, the computational cost will be analyzed for PMVLD.
For HAF-ICPF, the main implementation procedures include HAF (O(N)) [2] and ICPF O N 3 [40], where N denotes signal length. Thus, the HAF-ICPF computational cost is O N 3 [30,33,34]. For the algorithm presented in [34], its main implementation procedures include a parametric autocorrelation function O N 2 , the GSCFT O N 2 log 2 N [34] and an FFT operation O N 2 log 2 N . Thus, its computational cost is O N 2 log 2 N [34]. For MLVD, the main implementation procedures include the parametric symmetric self-correlation function O N 2 , the keystone transform O N 2 log 2 N and the FFT operation O N 2 log 2 N . Thus, the MLVD computational cost is also O N 2 log 2 N [33]. For 2D-PMLVD, the main implementation procedures include the defined parametric symmetric self-correlation function O N 2 , the mSFT t,i − FFT τ,i O N 2 log 2 N , and (L − 1) products of frequency domain. However, L is usually a small number (e.g., 2), meaning that the additional cost is not excessive. Thus, the computational cost of 2D-PMLVD is O N 2 log 2 N . Table 1 lists the computational costs of these four estimation algorithms.

HAF-ICPF
O N 3 The algorithm presented in [34] O N 2 log 2 N MLVD O N 2 log 2 N 2D-PMLVD O N 2 log 2 N

Conclusions
In this paper, we propose 2D-PMLVD for QFM signals by employing multi-scale PSSF and mSFT t,i − FFT τ,i . Due to the characteristics of multi-scale PSSF, 2D-PMLVD can utilize mSFT t,i − FFT τ,i to remove coupling and accumulate energy to the same position under different scale factors. Based on the products of different 2-D frequency domains, cross terms are suppressed and auto-terms are enhanced. The main contributions of the proposed method include the following: (1) it eliminates the procedure of parameter searching; (2) it improves identifiability for multi-QFM signals; (3) it improves anti-noise performance; and (4) it achieves a balance between low computational cost and high performance with low SNR.