An Efficient ISAR Imaging of Targets with Complex Motions Based on a Quasi-Time-Frequency Analysis Bilinear Coherent Algorithm

The inverse synthetic aperture radar (ISAR) imaging for targets with complex motions has always been a challenging task due to the time-varying Doppler parameter, especially at the low signal-to-noise ratio (SNR) condition. In this paper, an efficient ISAR imaging algorithm for maneuvering targets based on a noise-resistance bilinear coherent integration is developed without the parameter estimation. First, the received signals of the ISAR in a range bin are modelled as a multicomponent quadratic frequency-modulated (QFM) signal after the translational motion compensation. Second, a novel quasi-time-frequency representation noise-resistance bilinear Radon-cubic phase function (CPF)-Fourier transform (RCFT) is proposed, which is based on the coherent integration of the energy of auto-terms along the slope line trajectory. In doing so, the RCFT also effectively suppresses the cross-terms and spurious peaks interference at no expense of the time-frequency resolution loss. Third, the cross-range positions of target’s scatters in ISAR image are obtained via a simple maximization projection from the RCFT result to the Doppler centroid axis, and the final high-resolution ISAR image is thus produced by regrouping all the range-Doppler frequency centroids. Compared with the existing time-frequency analysis-based and parameter estimation-based ISAR imaging algorithms, the proposed method presents the following features: (1) Better cross-term interference suppression at no time-frequency resolution loss; (2) computationally efficient without estimating the parameters of each scatters; (3) higher signal processing gain because of 2-D coherent integration realization and its bilinear function feature. The simulation results are provided to demonstrate the performance of the proposed method.


Introduction
Thanks to the ability to produce high-resolution microwave imagery for the non-cooperative target nearly regardless of weather condition, the inverse synthetic aperture radar (ISAR) presents a range of applications in the field of national defense surveillance [1][2][3][4]. The range-Doppler (RD) ISAR The rest of the paper is organized as follows. In Section 2, the characteristic of the received signal for the maneuvering target is discussed. In Section 3, the RCFT derivation and the ISAR imaging algorithm of the maneuvering target based on RCFT are presented. The experimental results and computational complexity analysis are given in Section 4. Section 5 is the conclusion of this paper.

ISAR Imaging Model of Maneuvering Target
In this section, we do not consider the motion compensation problem, which means that the standard range alignment and phase adjustment are implemented beforehand. In Figure 1, the ISAR geometric configuration of a target with the complex motion is depicted, where XOY is a Cartesian coordinate and the origin O is the position of target rotating center. In Figure 1, the R 0 denotes the initial distance of origin O to radar platform, and v r , a r , and γ r respectively represent the radial velocity, the acceleration, and the acceleration rate. With those notations, the instantaneous rotation angle where t m indicates the slow time variable, ω(t m ) is the rotation angular velocity at time t m , and Ω, Ω and Ω respectively denote the initial angular velocity, angular acceleration, and angular acceleration rate.
From (1), R s (t m ) is the slant range of the scatter P with position x p , y p at time t m , given by In practice, the coherent integration time is usually short, and therefore the target rotation angle during that time is small, i.e., 3-5 • [8,9]. From that, sin θ(t m ) and cos θ(t m ) are approximated by θ(t m ) and 1. With also the assumption that motion compensation is completed, the received azimuth signal in a range bin is where f c , c, A k and K respectively represent the carrier frequency, the velocity of the wave propagation, the magnitude of the kth point scatterer, and the number of point scatterers in one range cell, n(t m ) is additive complex white Gaussian noise with a variance of δ 2 . From (3), it is seen that the received signal in a range cell is a multi-component PPS. In a generic form, (3) is rewritten as A k exp j2π b k,0 + b k,1 t m + b k,2 t 2 m + +b k,3 t 3 m + n(t m ) where φ k (t m ) is the phase term, and b k,0 , b k,1 , b k,2 and b k,3 respectively represent the initial phase, the centroid frequency, the chirp rate, and the change rate of chirp rate. From (4) again, it is clearer that the azimuth signal in a range bin is a multicomponent QFM signal. In this case, the conventional FT is not suitable to process QFM signal, and the second and high-order terms in (4) will deteriorate the image if they are not properly handled. Therefore, in this work, a novel RCFT is developed as a substitution of FT in the azimuth focusing. where ( ) is the phase term, and , , , , , and , respectively represent the initial phase, the centroid frequency, the chirp rate, and the change rate of chirp rate,. From (4) again, it is clearer that the azimuth signal in a range bin is a multicomponent QFM signal. In this case, the conventional FT is not suitable to process QFM signal, and the second and high-order terms in (4) will deteriorate the image if they are not properly handled. Therefore, in this work, a novel RCFT is developed as a substitution of FT in the azimuth focusing.

Proposed Algorithm Description
For radar imaging of non-uniformly rotating target, the radar echo signals in a range cell can be characterized as multicomponent QFM signals after migration compensation, which has been illustrated in Section 2. In this section, an efficient ISAR imaging of targets with complex motions based on a noise-resistance bilinear coherent RCFT is proposed.

Description of the Proposed RCFT
In the CPF approach proposed in [17,18], instantaneous autocorrelation function of (4) is defined as cross-terms auto-terms +R s,n-terms (t m ; τ m ) noise-terms (5) where τ m is the lag time variable. The R s,c-terms (t m ; τ m ) and R s,n-terms (t m ; τ m ) are the cross-terms and the noise terms, respectively, and detailed expressions can be found in [17]. Taking the Fourier transform (FT) to (5) along the lag time variable τ m , one can obtain the CPF auto-terms where f τ 2 m is the frequency variable corresponding to the lag variable τ m and δ(·) is the Dirac delta function. The CPF s,c-terms t m ; f τ 2 m and CPF s,n-terms t m ; f τ 2 m are the cross term and the noise term after the FT.
Note also that from (6), the CPF performs the discrete FT (DFT) in terms of the lag time variable τ 2 m whose sampling grid is non-uniform. As a result, the traditional uniform-sampling based efficient fast FT (FFT) cannot be directly applied. To calculate the CPF, the non-uniform discrete Fourier transform (NUDFT) is usually utilized, and the cost of the direct usage of NUDFT is as high as O N 2 τ m with the signal length N τ m . To reduce the computational cost of the NUDFT, the non-uniform FFT (NUFFT) [19] is preferred, and its computational cost of performing the FT along the τ m axis is O(2N τ m log 2 N τ m ) without the performance loss. For more information, the detailed implementation procedures on NUFFT can be found in [19]. Applying the NUFFT in (5) produces the CPF that is the same as in (6) as where NUFFT τ 2 m is NUFFT operator along the lag time variable τ 2 m . In (6), because of the nonlinear coupling between azimuth slow-time variable t m and the lag time variable τ 2 m , the energy of the auto-terms concentrates along the line of f τ 2 m = b k,2 + 3b k,3 t m in the slow time-Doppler frequency t m − f τ 2 m plane. Utilizing this line, the third-and second-order coefficients are extracted by the slope and the y-intercept [17]. However, in the case of multi-component QFM signal, the identifiability of the CPF is of problem because of cross-terms and spurious peaks [17,18]. To visually demonstrate this issue, the CPF of a two-component QFM signal with length N = 2024 with parameters Figure 2. From Figure 2a,b, a sharp spurious peak is of presence, and the energy of auto-terms is focused along the slope lines (i.e., f τ 2 m = b k,2 + 3b k,3 t m , k = 1, 2), whereas the energy of cross-terms is diffused in the t m , f τ 2 m domain since their positions vary with time t m . Generally speaking, for a K-component QFM signal, there are K 2 − K cross-terms and K 2 − K /2 spurious peaks [17], which pose a serious interference to the estimation and detection of the auto-terms, unless they are properly reduced. In (6), because of the nonlinear coupling between azimuth slow-time variable and the lag time variable , the energy of the auto-terms concentrates along the line of = , + 3 , in the slow time-Doppler frequency − plane. Utilizing this line, the third-and second-order coefficients are extracted by the slope and the y-intercept [17]. However, in the case of multicomponent QFM signal, the identifiability of the CPF is of problem because of cross-terms and spurious peaks [17,18]. To visually demonstrate this issue, the CPF of a two-component QFM signal with length N = 2024 with parameters σ = σ = 1, , = −0.1, , = 5 × 10 and , = −1 × 10 , = 0.1, , = −1 × 10 , and , = 2 × 10 is shown in Figure 2. From Figure 2a,b, a sharp spurious peak is of presence, and the energy of auto-terms is focused along the slope lines (i.e., = , + 3 , , = 1,2), whereas the energy of cross-terms is diffused in the , domain since their positions vary with time . Generally speaking, for a K-component QFM signal, there are − cross-terms and ( − )/2 spurious peaks [17], which pose a serious interference to the estimation and detection of the auto-terms, unless they are properly reduced. To overcome the identifiability problem of the CPF, the cross-terms and spurious peaks must be properly reduced. To this goal, two approaches of the Radon-CPF transform (RCT) in [18] and Hough generalized high-order ambiguous function (Hough-GHAF) method in [20] have been developed to suppress the cross-terms, spurious peaks for multicomponent QFM signal. To utilize the energy of auto-terms, for the RCT, the integration was performed along the slope line defined by the polar  To overcome the identifiability problem of the CPF, the cross-terms and spurious peaks must be properly reduced. To this goal, two approaches of the Radon-CPF transform (RCT) in [18] and Hough generalized high-order ambiguous function (Hough-GHAF) method in [20] have been developed to suppress the cross-terms, spurious peaks for multicomponent QFM signal. To utilize the energy of auto-terms, for the RCT, the integration was performed along the slope line defined by the polar distance ρ T (radius) from the origin and polar angle θ T formed by the perpendicular to the line in the radon domain. The RCT method S RCT (ρ T ; θ T ) is given by [16] where Ψ t m represents the RCT operator on CPF t m , f τ 2 m . From (8), the integral accumulates all the energy of auto-terms and suppresses the cross-terms and spurious peaks.
The similar idea based on multilinear function of fourth-order GHAF is also adopted by Hough-GHAF method in [20]. The Hough-GHAF is defined by where Θ t m represents the Hough-GHAF transform operator. The detailed definition of the GHAF(t m , f τ m ) can be found in [18]. It is worth mentioning that, the Hough-GHAF is based on a multilinear function of fourth-order, and thus its SNR threshold is higher than RCT. From both (8) and (9), although the energy of auto-terms is explored, the operations in RCT and Hough-GHAF are not coherent, and therefore, the suppression ability of the cross-terms and noise is still not adequate.
To perform coherent integration, the RCFT is developed that fully exploits the energy of the auto-terms. However, the difficulty of performing coherent integration comes from the fact that the cubic and quadratic power terms of t m are present in auto-terms of (6). They are must be eliminated first because the peaks of auto-terms would be unfocused if the direct integration along t m were utilized. In what follows, a novel RCFT algorithm is proposed to coherently integrate the energy of auto-terms along slope lines.
First, to eliminate the effect of the quadratic power terms of t m in (6), here, we intelligently exploit the idea of the sampling property of the Dirac delta function, which is where g(t m ) is a general function of the variable t m , and t c denotes a fixed time index. According to (10), in order to utilize sampling property of the Dirac delta function, an appropriate phase term function should be designed to construct a same expression with the delta function. With this thinking, a modified CPF (MCPF) by utilizing the sampling property of the delta function is designed as . It is now clear from (11) that the negative effects of the quadratic power of t m in auto-terms are removed by means of the sampling property of the Dirac delta function. This is the first important step in the proposed RCFT algorithm.
It is also found from (11) that the cubic power term of t m just corresponds to the slope of the auto-terms energy distribution in the time-frequency plane. Therefore, to eliminate the effect of the cubic power terms of t m in (5), inspired by the Radon-Fourier-transform (RFT) in [21,22], a novel RCFT algorithm is defined by where Γ t m denotes the RCFT operator. The (12) is a novel transform kernel function that is given by where f t m is the frequency variable with respect to t m . Note that in the case of zero cubic term, i.e., b 3 = 0, the proposed RCFT reduces to the coherent integrated CPF (CICPF) approach proposed in our previous work [9]. Moreover, when both the second-and three-order terms are zeros, i.e., b 3 = b 2 = 0, the proposed RCFT becomes the FT, which indicates that FT is a special case of the RCFT. Substituting (11) into (12) and after reassigning yields where σ 2 In (14), when the slope of the searching slope 3 , the cubic power of t m in auto-terms is eliminated, also demonstrated by RFT method [21]. Therefore, the proposed RCFT in (14) is capable of realizing the coherent integration for auto-terms while suppressing the cross-terms and spurious peaks. This is the second important step in the proposed RCFT algorithm. Moreover, when the searching slope line fully overlaps with the energy distribution slope line of the auto-terms, namely tan(θ T ) = 3b k,3 and ρ T = b k,2 cos θ T , the proposed RCFT maximizes the output energy of auto-terms and produces a distinct peak in which the maximal output energy E max is calculated by where G FT is the FT coherent integration gain. Meanwhile, when tan(θ T ) = 3b k,3 or ρ T = b k,2 cos θ T , the energy of the auto-terms integration S RCFT (ρ T ; θ T ) E max due to incoherent integration or the fact that only part of the auto-terms energy is accumulated. Figure 2c,d depicts the coherent where only the auto-terms are accumulated into peaks, while the cross-terms and spurious peaks are almost completely suppressed. Although Figure 2c shows that there also exist the cross-terms on the RCFT plane, they are much smaller compared to the auto-terms. Unlike the traditional quadratic time-frequency distributions that usually exploit the smoothing, optimal kernel design or nonlinear filtering techniques to moderately tradeoff between the cross-terms and resolution, the proposed RCFT, on its own, can greatly "suppress" the cross-terms without any resolution loss. By "suppress", we mean a relative suppression is achieved since the RCFT greatly strengthens the energy of the auto-terms instead of suppressing the cross-terms directly.
From (12), interestingly, the proposed RCFT has similar operations as RCT and Hough-GHAF, and they all use the energy of the auto-terms along time-frequency trajectory in t m − f τ 2 m domain. The main difference is the coherent accumulation developed in the proposed RCFT. Therefore, the RCFT will definitely outperform the existing methods in complex environments via coherent integration operation. The detailed procedure of the proposed RCFT is shown in Figure 3. RCFT will definitely outperform the existing methods in complex environments via coherent integration operation. The detailed procedure of the proposed RCFT is shown in Figure 3. The differences and advantages of RCFT compared with others approaches are briefly summarized as follows. Remark 1: The RCFT employs the merits of both RCT and FT, and it not only has the same integration time as RCT but also works well as a useful tool for nonstationary signals. Remark 2: The bilinear cubic phase function in (5) utilizes only one time correlation, which is viewed as a signal energy preservation because each additional one time correlation loses about 4 -5 dB in the SNR threshold [17]. In addition to that, the 2-D coherent integration realized in the proposed RCFT will further enhance the SNR. Therefore, the proposed RCFT algorithm provides a good performance, especially when the SNR is low, see simulation section. Remark 3: the NUFFT speeds up the Fourier transform along the non-uniformly spaced lag-time axis, which is helpful for algorithm real-time realization.

Numerical study of RCFT
In this section, a simulation with the same parameters in Figure 2 is provided to demonstrate the performance comparisons with the RCT, and Hough-GHAF methods under different SNR conditions. Figure 4a and Figure 4b depict the comparisons of RCT, Hough-GHAF, and the proposed RCFT with SNR being 5 dB and -5 dB, respectively. In Figure 4a, the cross-terms and spurious peaks are suppressed by all three approaches. However, both the RCT and Hough-GHAF are sensitive to noise. In Figure 4b, when SNR is low, say −5 dB, the spectrum of the Hough-GHAF is overwhelmed by the noise due to the non-coherent integration and multilinear function of fourth-order utilization. On the other hand, the RCT and the proposed RCFT are able to generate two distinct peaks at true locations. However, the sidelobes of RCFT are much lower than that of RCT due to the 2-D coherent integration, which indicates that the RCFT presents better cross-terms and noise suppression ability.  The differences and advantages of RCFT compared with others approaches are briefly summarized as follows.
Remark 1: The RCFT employs the merits of both RCT and FT, and it not only has the same integration time as RCT but also works well as a useful tool for nonstationary signals. Remark 2: The bilinear cubic phase function in (5) utilizes only one time correlation, which is viewed as a signal energy preservation because each additional one time correlation loses about 4-5 dB in the SNR threshold [17]. In addition to that, the 2-D coherent integration realized in the proposed RCFT will further enhance the SNR. Therefore, the proposed RCFT algorithm provides a good performance, especially when the SNR is low, see simulation section. Remark 3: the NUFFT speeds up the Fourier transform along the non-uniformly spaced lag-time axis, which is helpful for algorithm real-time realization.

Numerical Study of RCFT
In this section, a simulation with the same parameters in Figure 2 is provided to demonstrate the performance comparisons with the RCT, and Hough-GHAF methods under different SNR conditions. Figures 4a,b depict the comparisons of RCT, Hough-GHAF, and the proposed RCFT with SNR being 5 dB and -5 dB, respectively. In Figure 4a, the cross-terms and spurious peaks are suppressed by all three approaches. However, both the RCT and Hough-GHAF are sensitive to noise. In Figure 4b, when SNR is low, say −5 dB, the spectrum of the Hough-GHAF is overwhelmed by the noise due to the non-coherent integration and multilinear function of fourth-order utilization. On the other hand, the RCT and the proposed RCFT are able to generate two distinct peaks at true locations. However, the sidelobes of RCFT are much lower than that of RCT due to the 2-D coherent integration, which indicates that the RCFT presents better cross-terms and noise suppression ability.
are suppressed by all three approaches. However, both the RCT and Hough-GHAF are sensitive to noise. In Figure 4b, when SNR is low, say −5 dB, the spectrum of the Hough-GHAF is overwhelmed by the noise due to the non-coherent integration and multilinear function of fourth-order utilization. On the other hand, the RCT and the proposed RCFT are able to generate two distinct peaks at true locations. However, the sidelobes of RCFT are much lower than that of RCT due to the 2-D coherent integration, which indicates that the RCFT presents better cross-terms and noise suppression ability.

ISAR Imaging for Maneuvering Target Based on the Proposed RCFT Algorithm
In this section, an efficient ISAR imaging algorithm of targets with complex motions is developed, in which the FT in conventional RD algorithm is replaced by the proposed quasi-time-frequency analysis bilinear coherent RCFT. The main steps of the ISAR imaging algorithm are listed as follows.
Step 1: Perform the range compression and the translational motion compensation including envelope alignment and phase autofocus. RCFT cterms ( f t m ; ρ T ; θ T ) in the Doppler Centroid f t m -polar radius ρ T -polar angle θ T domain.
Step 4: Project the three-dimensional data matrix RCFT( f t m ; ρ T ; θ T ) onto the Doppler frequency axis along the polar radius ρ T and polar angle θ T , which is obtained by It is a common knowledge that the Doppler frequency of each scatterer is proportional to its cross-range position in the target. Hence, the cross-range ISAR image of the target can be obtained by its projection onto the Doppler frequency axis. In Figure 5, the projection results of the RCFT shown in the Figure 2 along the Doppler frequency dimension is presented. It is seen that two distinct peaks appear along the Doppler frequency axis, which corresponds to two scatterer cross-range positions.
It is a common knowledge that the Doppler frequency of each scatterer is proportional to its cross-range position in the target. Hence, the cross-range ISAR image of the target can be obtained by its projection onto the Doppler frequency axis. In Figure 5, the projection results of the RCFT shown in the Figure 2 along the Doppler frequency dimension is presented. It is seen that two distinct peaks appear along the Doppler frequency axis, which corresponds to two scatterer cross-range positions.

Figure 5. Maximization projection onto the Doppler frequency dimension.
Step 5: Set a proper extraction threshold or a filter to suppress the residual cross terms and noise in the Doppler centroid frequency dimension. In practice, the threshold is usually determined by subtracting −3～−4.5 dB from the maximal energy.
Step 6: Repeat the process of step1-step5 for all range cells, and the final high-resolution ISAR image is thus produced by regrouping all the range-Doppler frequency centroids. Since the proposed method does not require computations such as parameter estimation for each scatterer, it is computationally more efficient than the similar parameter estimation based algorithms. The flowchart of the proposed ISAR imaging algorithm is shown in Figure 6. Step 5: Set a proper extraction threshold or a filter to suppress the residual cross terms and noise in the Doppler centroid frequency dimension. In practice, the threshold is usually determined by subtracting −3~−4.5 dB from the maximal energy.
Step 6: Repeat the process of step 1-step 5 for all range cells, and the final high-resolution ISAR image is thus produced by regrouping all the range-Doppler frequency centroids. Since the proposed method does not require computations such as parameter estimation for each scatterer, it is computationally more efficient than the similar parameter estimation based algorithms. The flowchart of the proposed ISAR imaging algorithm is shown in Figure 6.

Components Computational Complexity Analysis
In this section, we analyze quantitatively the computational complexity of our proposed algorithm. For comparison purposes, the RCD method in [15] where the high-resolution ISAR image can be obtained based a new quasi-time-frequency transform named Lv's distribution using the second-order phase model under a low SNR environment. Moreover, the cross-terms suppression in the RCD method is better achieved with no time-frequency resolution loss. On the other hand, to compare with the parameter estimation-based ISAR imaging method, the CIGCPF-CICPF algorithm [9] where it is recently proposed for maneuvering target ISAR imaging and parameter estimation with third-order motion model in low SNR condition.
In this comparison, for the illustration conveniences, assume that the range compression and

Components Computational Complexity Analysis
In this section, we analyze quantitatively the computational complexity of our proposed algorithm. For comparison purposes, the RCD method in [15] where the high-resolution ISAR image can be obtained based a new quasi-time-frequency transform named Lv's distribution using the second-order phase model under a low SNR environment. Moreover, the cross-terms suppression in the RCD method is better achieved with no time-frequency resolution loss. On the other hand, to compare with the parameter estimation-based ISAR imaging method, the CIGCPF-CICPF algorithm [9] where it is recently proposed for maneuvering target ISAR imaging and parameter estimation with third-order motion model in low SNR condition.
In this comparison, for the illustration conveniences, assume that the range compression and translational motion compensation have been completed. The computational complexities of above-mentioned two methods and our proposed algorithm are quantitatively provided. In general, an N-point FFT or inverse FFT (IFFT) needs 5N log 2 (N) floating-point operations (FLOPs) and one-time complex multiplication needs 6N FLOPs. In what follows, N r and N a are respectively used to denote the number of range cells and the number of azimuth pulses, K l is the target scatterer number in the lth range cell, N t m is the signal length t m , and N τ m is used to represent the length of the lag variable τ m .
For the RCD method, its implementation steps mainly include performing the quasi-time-frequency distribution (Lv's distribution) to each range cell. Take a range cell processing procedure for example, the complex multiplication in constructing the symmetric instantaneous autocorrelation function matrix with computational complexity of O(6N a N τ m ), the FFT operation along the lag-time variable axis with computational cost of O(5N a N τ m log 2 N τ m ), a keystone transform adopted to remove the coupling terms with complexity O(2(2N ker − 1)N a N τ m ), where N ker is the length of the interpolation operation kernel, the FFT operation along the scaling slow time variable with computational cost of O(5N a N τ m log 2 N a ), and neglecting other relatively small computational complexity operation steps. Therefore, the total computational complexity of the RCD algorithm [15] is Compared with the parameter estimation-based method and our proposed method, which will be discussed later, the time-frequency analysis-based RCD method has a great advantage in terms of computational complexity. However, this method suffers from imaging performance degradation without considering the third-order phase effects.
For the parameter estimation-based CIGCPF-CICPF method where it needs to estimation each scatterer parameter, the computational load consists of the following steps. Take one scatterer estimation for example, to estimate the first-and third-order coefficient using the CIGCPF, the complex multiplication in constructing fourth-order multilinear GCPF function matrix with computational complexity of O(18N a N τ m ), the NUFFT operation along the lag-time variable axis with computational cost of O(40N a N τ m log 2 N τ m ), one time compensation function multiplication with complexity O(6N a N τ m ), the FFT operation along the slow time variable with computational cost of O(5N a N τ m log 2 N a ). Then one Dechirping operation is required, which needs one N a -dimensional complex multiplication. Second, to obtain second-order coefficient using the CICPF, the computational complexity requirement is similar to the CIGCPF operation. Finally, one N a -dimensional FFT is needed to estimate amplitude. Therefore, the total computational complexity of the CIGCPF-CICPF method [9] for maneuvering target imaging with third-order phase model is Similar to the RCD method, the proposed ISAR imaging algorithm is also based the quasi-time-frequency analysis named RCFT. According to the imaging steps and the flowchart of the proposed algorithm in Figure 6, the proposed algorithm implementation procedures mainly include applying the proposed RCFT to each range cell. Take a range cell processing procedure for example, in constructing bilinear CPF function matrix with computational complexity of O (6N a N τ m ) According to the above analysis, the computational complexity of the proposed ISAR imaging is higher than that of the RCD method, but still much lower that of parameter estimation-based CIGCPF-CICPF approach. In conclusion, the proposed method may well achieve a trade-off between the computational complexity and the imaging performance, see performance analysis section.

Simulation Results and Analysis
To confirm the validity of the proposed algorithm, simulation experiments are conducted now under two conditions of input SNR = 5 dB and SNR = −3 dB with the simulation parameters of the radar and the moving target listed in Table 1. The target scatterer model used in the simulation was a ship with 49 scatterers, shown in Figure 7. Effective rotational motion angular velocity 0.018 rad/s Effective rotational motion acceleration 0.008 rad/s 2 Effective rotational motion acceleration rate 0.002 rad/s 3 In Figure 8, the range compression and imaging result using RD algorithm without the migration compensation are presented. It is clear that from Figure 8a the energy for the target disperses several range bins. In addition, the image produced by the RD algorithm is unclear and blurry. In Figure 9, the motion compensation including the translational compensation (envelope alignment and phase adjustment) and the migration through resolution cells (MTRCs) correction [9,[23][24][25], is conducted. From Figure 9, it reveals that the range migration is fully rectified and the energy of the target is now concentrated into one range bin. It is therefore concluded that the migration compensation is effective must, and this step is always conducted first in the following experiments. In the following, the comparison results of different approaches are presented. In Figure 8, the range compression and imaging result using RD algorithm without the migration compensation are presented. It is clear that from Figure 8a the energy for the target disperses several range bins. In addition, the image produced by the RD algorithm is unclear and blurry. In Figure 9, the motion compensation including the translational compensation (envelope alignment and phase adjustment) and the migration through resolution cells (MTRCs) correction [9,[23][24][25], is conducted.
From Figure 9, it reveals that the range migration is fully rectified and the energy of the target is now concentrated into one range bin. It is therefore concluded that the migration compensation is effective must, and this step is always conducted first in the following experiments. In the following, the comparison results of different approaches are presented. In Figure 8, the range compression and imaging result using RD algorithm without the migration compensation are presented. It is clear that from Figure 8a the energy for the target disperses several range bins. In addition, the image produced by the RD algorithm is unclear and blurry. In Figure 9, the motion compensation including the translational compensation (envelope alignment and phase adjustment) and the migration through resolution cells (MTRCs) correction [9,[23][24][25], is conducted. From Figure 9, it reveals that the range migration is fully rectified and the energy of the target is now concentrated into one range bin. It is therefore concluded that the migration compensation is effective must, and this step is always conducted first in the following experiments. In the following, the comparison results of different approaches are presented.  Figure 9b. Since the Doppler frequency is time-varying, it is observed that the image is seriously blurred in the cross range. The images produced by STFT, WVD, and SPWVD are demonstrated in Figure 10a-c, respectively. Because of the cross-terms interference in WVD, the shape of ship is not visible. Suffering from low resolutions, the images from STFT and SPWVD are not focused well and each scatter spreads out in azimuth direction. In Figure 10d, the ISAR imaging result obtained by the RCD method in [15] is presented. As discussed earlier, since this method is only valid for LFM signal, the image obtained is seriously blurred in the cross range. To clearly demonstrate the advantages of the proposed approach over the existing methods, we also compare the proposed method with the parameter estimation-based algorithms. Figure 10e,f provide the 2-D SAR images obtained by the IHAF-ICPF method [7], CIGCPF-CICPF method [9] at SNR = 5 dB. It is observed that ship shape is generated by those two approaches, but they are not free of residual interference. On another note, parameter estimation-based algorithms are computationally extensive and suffer from model mismatch issue, which limit their usages. Figure 10g shows the high-quality ISAR image produced by the proposed quasi-time-frequency RCFT, where image is well focused in both range and azimuth directions, which agrees with our theoretical analysis. Figure 11 presents the  Figure 9b. Since the Doppler frequency is time-varying, it is observed that the image is seriously blurred in the cross range. The images produced by STFT, WVD, and SPWVD are demonstrated in Figure 10a-c, respectively. Because of the cross-terms interference in WVD, the shape of ship is not visible. Suffering from low resolutions, the images from STFT and SPWVD are not focused well and each scatter spreads out in azimuth direction. In Figure 10d, the ISAR imaging result obtained by the RCD method in [15] is presented. As discussed earlier, since this method is only valid for LFM signal, the image obtained is seriously blurred in the cross range. To clearly demonstrate the advantages of the proposed approach over the existing methods, we also compare the proposed method with the parameter estimation-based algorithms. Figure 10e,f provide the 2-D SAR images obtained by the IHAF-ICPF method [7], CIGCPF-CICPF method [9] at SNR = 5 dB. It is observed that ship shape is generated by those two approaches, but they are not free of residual interference. On another note, parameter estimation-based algorithms are computationally extensive and suffer from model mismatch issue, which limit their usages. Figure 10g shows the high-quality ISAR image produced by the proposed quasi-time-frequency RCFT, where image is well focused in both range and azimuth directions, which agrees with our theoretical analysis. Figure 11 presents the ISAR images obtained by the aforementioned methods under SNR = −3 dB. Compared with other methods, under the low SNR condition, our proposed ISAR imaging algorithm still produces high quality image. This means most scatterers are relocated correctly with less scatterers loss and less artifact appearances as presented in Figure 11g. This superior performance obtained is all because the 2-D coherent integrations and bilinear features utilized in the proposed methods to suppress the cross-terms and spurious peaks, and at the same time resolution is not sacrificed.
To further evaluate the performance, the entropy of 2-D ISAR images is utilized to measure the image quality. It is well known that a better quality image indicates a smaller entropy [26][27][28]. The entropy for an image g(m, n) is where S = ∑ M−1 m=0 ∑ N−1 n=0 |g(m, n)| 2 . Entropies obtained by these methods are listed in Table 2, where the smallest values are obtained by the proposed imaging algorithm. This is in the agreement with the conclusions obtained in Figures 10 and 11 since the proposed method produces the clearest ship model free of artifacts, and also demonstrate the effectiveness of the proposed ISAR imaging algorithm under low SNR condition. Table 2. Entropies of ISAR Images in Figure 10 and Figure 11.  (g) Figure 10. ISAR images of the simulated data under SNR = 5 dB. (a) STFT algorithm in [10]; (b) WVD algorithm in [12]; (c) SPWVD method in [13]; (d) RCD method in [15]; (e) IHAF-ICPF method in [7]; (f) CIGCPF-CICPF method in [9]; (g) our proposed method.  [10]; (b) WVD algorithm in [12]; (c) SPWVD method in [13]; (d) RCD method in [15]; (e) IHAF-ICPF method in [7]; (f) CIGCPF-CICPF method in [9]; (g) our proposed method.
In order to imitate the measured data environment, we have added a simulation analysis. In this simulation, the amplitudes of scatterers vary randomly from 0.1 to 1, and the corresponding results are provided in Figure 12 It clearly shows that the few weak scatterers are not prominent anymore, but most scatterers are still present to form the shape of the target. In fact, all the nonlinear timefrequency analysis approaches are actually challenged by the issue of weak scatterers.  [10]; (b) WVD algorithm in [12]; (c) SPWVD method in [13]; (d) RCD method in [15]; (e) IHAF-ICPF method in [7]; (f) CIGCPF-CICPF method in [9]; (g) our proposed method.
In order to imitate the measured data environment, we have added a simulation analysis. In this simulation, the amplitudes of scatterers vary randomly from 0.1 to 1, and the corresponding results are provided in Figure 12 It clearly shows that the few weak scatterers are not prominent anymore, but most scatterers are still present to form the shape of the target. In fact, all the nonlinear time-frequency analysis approaches are actually challenged by the issue of weak scatterers.
In order to imitate the measured data environment, we have added a simulation analysis. In this simulation, the amplitudes of scatterers vary randomly from 0.1 to 1, and the corresponding results are provided in Figure 12 It clearly shows that the few weak scatterers are not prominent anymore, but most scatterers are still present to form the shape of the target. In fact, all the nonlinear timefrequency analysis approaches are actually challenged by the issue of weak scatterers. (c) (d) Figure 12. ISAR images of the simulated data under SNR = 5 dB in the case of uneven strengths of scatterers. (a) WVD algorithm in [12]. (b) SPWVD method in [13]. (c) RCD method in [15]. (d) Our proposed method.
The implementation times of different methods are provided in Table 3, which are obtained using a computer with double core CPU 3.4 GHz and memory of 8 G. From the table, the running times of time-frequency analysis (STFT, WVD, and SPWVD methods)-based approaches are lowest and the parameter estimation-based ISAR imaging approaches of IHAF-ICPF method and CIGCPF-CICPF are highly time consuming. The time cost of proposed RCFT is higher than the time-frequency analysis based approaches, but much lower than the parameter estimation-based approaches. Taking both the performance and computation into consideration, to conclude, the proposed method well trades off them in producing clear image of the maneuvering target at low SNR environment.

Conclusions
In this paper, the RCFT method is first proposed for the analysis of multiple QFM signals. Because of 2-D coherent integration realization and the bilinear function feature, the better cross-term interference suppression is achieved at no loss of time-frequency resolution, and also higher signal processing gain is obtained. The RCFT is computationally efficient since no expensive threedimensional (3-D) parameter search is required. After that, the RCFT is applied to the ISAR imaging  [12]. (b) SPWVD method in [13]. (c) RCD method in [15]. (d) Our proposed method.
The implementation times of different methods are provided in Table 3, which are obtained using a computer with double core CPU 3.4 GHz and memory of 8 G. From the table, the running times of time-frequency analysis (STFT, WVD, and SPWVD methods)-based approaches are lowest and the parameter estimation-based ISAR imaging approaches of IHAF-ICPF method and CIGCPF-CICPF are highly time consuming. The time cost of proposed RCFT is higher than the time-frequency analysis based approaches, but much lower than the parameter estimation-based approaches. Taking both the performance and computation into consideration, to conclude, the proposed method well trades off them in producing clear image of the maneuvering target at low SNR environment.

Conclusions
In this paper, the RCFT method is first proposed for the analysis of multiple QFM signals. Because of 2-D coherent integration realization and the bilinear function feature, the better cross-term interference suppression is achieved at no loss of time-frequency resolution, and also higher signal processing gain is obtained. The RCFT is computationally efficient since no expensive three-dimensional (3-D) parameter search is required. After that, the RCFT is applied to the ISAR imaging problem for a maneuvering target in producing the clear image. Numerical results demonstrate that the proposed method outperforms existing ISAR imaging algorithms in terms of both visual inspections and objective performance measures.