Senior Member, IEEE Yurduseven
email: okan.yurduseven@qub.ac.uk

Masoud Amir

Philipp Del Molaei
email: a.molaei@qub.ac.uk

Vincent Hougne

Fellow, IEEE Okan Fusco
email: v.fusco@ecit.qub.ac.uk

Vincent Molaei

Fusco

Philipp Del
email: philipp.del-hougne@univ-rennes1.fr

Efficient Joint Estimation of DOA, Range and Reflectivity in Near-Field by Using Mixed-Order Statistics and a Symmetric MIMO Array

Keywords: DOA-range-reflectivity estimation, MIMO array, MOSA, near-field

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

I N recent decades, active multiple-input multiple-output (MIMO) array systems have been widely studied, especially in areas related to radar and sonar target detection and localization [START_REF] Xiong | Fda-mimo radar range-angle estimation: Crlb, mse, and resolution analysis[END_REF][START_REF] Pan | Robust time-reversal is combined with distributed multiple-input multiple-output sonar for detection of small targets in shallow water environments[END_REF][START_REF] Maresca | Photonics for coherent mimo radar: An experimental multi-target surveillance scenario[END_REF]. By using waveform diversity techniques, these systems can offer many advantages over traditional phased-array radar systems, including better practical parameter identifiability, higher accuracy of parameter estimation, and flexible transmit beampattern design [START_REF] Gong | Adaptive weight matrix design and parameter estimation via sparse modeling for mimo radar[END_REF]. There are generally two different categories for configuring antennas in a MIMO radar system. The first one is the distributed MIMO radar [START_REF] Liu | Persymmetric adaptive target detection with distributed mimo radar[END_REF], in which the distance between transmitting and/or receiving antennas is large. In this case, the different view-angles from the antennas to targets are used to improve the target detection performance by capturing the spatial diversity of radar crosssection with non-coherent processing [START_REF] Sadeghi | Target localization geometry gain in distributed mimo radar[END_REF]. Oppositely, the other one, adopted in this paper, is the co-located MIMO radars [START_REF] Zhao | Mimo transmit beampattern matching under waveform constraints[END_REF][START_REF] Haghnegahdar | Interference cancellation in co-located mimo radars using waveform optimisation in signal dependent clutter[END_REF], in which the spacing between antennas is comparable with the wavelength. In this case, the waveform diversity is exploited to improve the estimation performance with the large virtual aperture [START_REF] Korkmaz | High-resolution, contiguous sar imaging using co-located mimo arrays: Experimental proof of concept[END_REF].

Direction of arrival (DOA) estimation plays a key role in modern signal processing for various applications such as radar, sonar, microphone arrays, wireless communications, electronic surveillance, medical diagnosis and treatment, radio astronomy and seismology [START_REF] Abdelbari | Direction of arrival estimation of wideband rf sources[END_REF]. Numerous DOA estimation techniques for classical single-input multiple-output (SIMO) setup have been reported in the literature [START_REF] Zheng | Classification and localization of mixed near-field and far-field sources using mixed-order statistics[END_REF][START_REF] Molaei | Components separation algorithm for localization and classification of mixed nearfield and far-field sources in multipath propagation[END_REF][START_REF] Molaei | A pure cumulant-based method with low computational complexity for classification and localization of multiple near and far field sources using a symmetric array[END_REF][START_REF] Abbasi | Frequency-diverse multimode millimetre-wave constant-ϵr lens-loaded cavity[END_REF][START_REF] Qin | Improved two-dimensional doa estimation using parallel coprime arrays[END_REF]. In recent years, however, the development of MIMO systems has opened up new opportunities in DOA estimation. In [START_REF] Zahernia | Music algorithm for doa estimation using mimo arrays[END_REF], the performance of the multiple signal classification (MUSIC) algorithm in MIMO arrays compared to conventional arrays was investigated by considering the parameters of the number of snapshots, the number of array elements and signal-to-noise ratio (SNR) values. Hayashi et al. [START_REF] Hayashi | Doa estimation in mimo radar using temporal spatial virtual array with music algorithm[END_REF], to detect and classify indoor human activities, developed a DOA estimation method for Doppler radars with a MIMO system. To improve the accuracy of DOA estimation, they proposed a temporal-spatial virtual array based on the Doppler shift of a moving target. Lonkeng et al. [START_REF] Lonkeng | Two-dimensional doa estimation using arbitrary arrays for massive mimo systems[END_REF] considered an arbitrary array geometry for the MIMO system. They studied the two-dimensional (2D) Fourier domain line search MUSIC algorithm to mutually estimate azimuth and elevation angles. In [START_REF] Cao | Doa estimation for multiple targets in mimo radar with nonorthogonal signals[END_REF], the problem of DOA estimation in MIMO radar with non-orthogonal signals has been addressed. In [START_REF] Chen | A mimo radar-based doa estimation structure using compressive measurements[END_REF], to reduce the number of channels, a DOA estimation method based on compressive sensing is presented for MIMO radar with co-located antennas. In the study [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF], Liu et al. have addressed the problem of DOA estimation in a MIMO radar system with orthogonal waveforms for fast-moving targets. They formulated a least absolute shrinkage and selection operator (LASSO) with an atomic norm to denoise the received signals. However, the denoising problem cannot be solved efficiently and requires additional steps; because it needs to be formulated first as . a dual optimization problem and then transformed into a semidefinite problem. All studies [START_REF] Zahernia | Music algorithm for doa estimation using mimo arrays[END_REF][START_REF] Hayashi | Doa estimation in mimo radar using temporal spatial virtual array with music algorithm[END_REF][START_REF] Lonkeng | Two-dimensional doa estimation using arbitrary arrays for massive mimo systems[END_REF][START_REF] Cao | Doa estimation for multiple targets in mimo radar with nonorthogonal signals[END_REF][START_REF] Chen | A mimo radar-based doa estimation structure using compressive measurements[END_REF][START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] make a far-field (FF) assumption and focus only on estimating the DOA parameter. Under this assumption, one or more plane waves impinge on the array.

Although plane-wave assumption can simplify modeling and processing, it is not valid for near-field (NF) applications; because when the target is located in the Fresnel zone of the aperture array, the shape of the spherical wavefront changes nonlinearly with the position of the array and is determined by both DOA and range parameters [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF]. A joint estimation of range, DOA and frequency is presented in [START_REF] Chen | A new algorithm for joint range-doa-frequency estimation of near-field sources[END_REF] in which the impinging signals are assumed to be complex exponential sequences received from third-party transmitters in the environment. In [START_REF] Zhi | Near-field source localization via symmetric subarrays[END_REF], a NF source localization algorithm based on conventional second-order statistics (SOS) is presented. By dividing the ULA into two symmetric subarrays, the steering vectors of the subarrays yield the property of rotational invariance in the signal subspace. A full estimation of signal parameters via rotational invariant techniques (ESPRIT) approach is also used for localizing with low computational complexity in [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF]. Methods [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF][START_REF] Chen | A new algorithm for joint range-doa-frequency estimation of near-field sources[END_REF][START_REF] Zhi | Near-field source localization via symmetric subarrays[END_REF] are only compatible with a passive scenario. A 2D-MUSIC processing for the simultaneous estimation of target angle and range is presented in [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF]. However, it has a very high computational complexity due to the use of a 2D spectral search as well as a spatial smoothing operation. In [START_REF] Singh | Near field targets localization using bistatic mimo system with spherical wavefront based model[END_REF], by considering the spherical wavefront, a subspace-based NF localization method is presented for the bistatic system in which the transmitter and receiver are separated by a distance comparable to the target distance. Based on an improved blind calibration algorithm, a NF source localization is proposed in [START_REF] Cui | Joint estimation of near-field source parameters and array response[END_REF]. Although the iterative approach used in [START_REF] Cui | Joint estimation of near-field source parameters and array response[END_REF] does not require precise initial conjectures, the convergence process of this approach has a high computational cost. A joint estimate of the angle and range applicable in the NF is presented in [START_REF] Xie | Parameters estimation of mixed far-field and near-field sources via second order statistics[END_REF]. In the first step, two correlation matrices are constructed to obtain DOA and frequency estimates. Then two more correlation matrices are formed to obtain the range value. The mechanism used in [START_REF] Xie | Parameters estimation of mixed far-field and near-field sources via second order statistics[END_REF] requires additional operations to pair the parameters. In [START_REF] Cao | Complex resnet aided doa estimation for near-field mimo systems[END_REF], by considering the problem as a regression task, a DOA estimation framework based on complex-valued deep learning is presented for short-range MIMO communication systems. Solving this regression task containing a massive number of variables is challenging since datasets need to capture many complicated feature representations. On the other hand, the problem structure assumes a passive detection that is limited to receivers.

In summary, and given the above, it can be concluded that the development of parameter estimation techniques for MIMO systems is still in progress. The majority of the studies carried out in the literature are limited to estimating only the DOA parameter, conventionally under the FF assumption. However, as noted, the joint estimation of DOA and range is an important challenge for NF applications. Backscatter is the portion of the outgoing radar signal that the target redirects directly back towards the radar antenna. The normalized measure of the radar return from a distributed target is called the backscatter coefficient. Backscatter coefficient estimation, as one of the most important parameters of radar systems (especially in imaging radars [START_REF] Marks | Fourier accelerated multistatic imaging: A fast reconstruction algorithm for multiple-inputmultiple-output radar imaging[END_REF][START_REF] Molaei | An efficient waveform diversity based on variational mode decomposition of coded beatfrequency shifted signals algorithm for multiple-input multiple-output millimetre-wave imaging[END_REF][START_REF] Molaei | Fourier compatible near-field multiple-input multiple-output terahertz imaging with sparse non-uniform apertures[END_REF]) has been neglected in these works. On the other hand, in most of the existing methods, the estimation of the parameters is based on the construction of the covariance matrix of the received signals and the Gaussian noise affects the entries of this matrix. Even in the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] that does not require such a matrix, noise cancellation is still a major challenge. Table I summarizes the main features of the existing methods and their limitations compared to this work.

In this paper, for the first time, joint estimation of the DOA, range and reflectivity is presented for a MIMO array in NF, where the effect of the spherical wave is nonnegligible. By employing the symmetry property in the array, first, two special spatial expressions based on fourth-order and SOS from the received data are defined so that the output of the first expression lacks targets' range information. By collecting all the sensor lags based on the first expression, a cross-cumulant matrix of the sensors is constructed. In two different approaches, one using eigendecomposition and the other employing the propagator method (PM) principle, the vectors required to estimate the DOAs are extracted. Next, the DOAs of the targets are estimated by a 1D spectral search. The second expression, which includes both angle and range information, is used to construct several covariance submatrices and then to construct a merged covariance matrix with a size corresponding to the dimensions of the array. The ranges are then estimated by using the corresponding estimated DOAs, applying eigendecomposition, and performing several 1D spectral searches. Finally, the equations required to estimate the reflectivity coefficients are derived (along with mathematical proofs). This derivation does not require any prior information such as noise power and is based only on data obtained from previous steps. The use of fourthorder cumulant (FOC) in the proposed method, in addition to improving the estimation accuracy, causes insensitivity to Gaussian noises. The proposed method, by separating the parameter estimation procedure into different phases, eliminates the very heavy multidimensional search. This will be very effective in significantly reducing the computational effort for real-time systems. Moreover, the estimated parameters are paired automatically and no additional steps are required.

The main contributions and novelties of this paper are summarized below:

■ The novelty in the scenario: for the first time in a practical scenario (a NF active MIMO system), the joint DOA, range and reflectivity estimates are provided. Works in the literature are either limited to DOA-only estimates or ultimately DOA and range estimates (usually in passive or SIMO scenarios).

■ The novelty in problem-solving:

• Unlike most existing works that use SOS in estimating parameters, this work uses mixed second-and fourth-order statistics (FOS).

• Definition of new cumulants: it is necessary to explain that cumulants can be defined in numerous forms. By changing the index of each of the signals in the function cum 4 {.}, one can reach a different output. In fact, the proper design and selection .

TABLE I A SUMMARY OF THE MAIN FEATURES OF EXISTING WORKS

Reference Method Basis

The Main Limitations Compared to This Work [START_REF] Zheng | Classification and localization of mixed near-field and far-field sources using mixed-order statistics[END_REF] Spatial cumulant (the cumulant output includes both DOA and range parameters) Applicable for the passive scenario, no reflectivity estimates [START_REF] Molaei | Components separation algorithm for localization and classification of mixed nearfield and far-field sources in multipath propagation[END_REF] Spatial cumulants (with the aim of separating non-coherent and coherent components) Applicable for the passive scenario, no reflectivity estimates [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] LASSO LASSO cannot be solved efficiently, DOA estimates only [START_REF] Chen | A new algorithm for joint range-doa-frequency estimation of near-field sources[END_REF] Spatial cumulants; ESPRIT-like estimator Applicable for the passive scenario, no reflectivity estimates [START_REF] Zhi | Near-field source localization via symmetric subarrays[END_REF] SOS Applicable for the passive scenario, no reflectivity estimates [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF] Spatial/temporal full ESPRIT-like approach Applicable for the passive scenario, no reflectivity estimates [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] 2D-MUSIC processing Very high computational complexity, no reflectivity estimates [START_REF] Singh | Near field targets localization using bistatic mimo system with spherical wavefront based model[END_REF] Cross-covariance matrices between the output data blocks Applicable for the bistatic system, no reflectivity estimates [START_REF] Cui | Joint estimation of near-field source parameters and array response[END_REF] Covariance matrix Iterative approach with high computational cost, applicable for the passive scenario, no reflectivity estimates [ of the sensors number to calculate the correlation between them is the most important factor in defining a cumulant so that the desired output is obtained and the information is sufficient to the end of the algorithm. For this purpose, a new special cumulant is defined in this paper. Cumulant output specific to the scenario of this paper is calculated for the first time, which is detailed in Appendix A. Note that the cumulant definitions available in the literature are for passive scenarios in which there is only one array (receiver array). Therefore, the selection of indices is limited to the elements of one array; while this work considers an active scenario. Therefore, the selection of indices is based on the interaction of the elements of the two arrays of transmitter and receiver.

• Mathematical separation of the data required to estimate each parameter, so as to avoid costly multidimensional estimates.

• Mathematical derivation of covariance and kurtosis of reflectivity coefficients.

■ Improving the performance of existing methods: the proposed method has high performance compared to the existing methods in reducing computational complexity and increasing estimation accuracy, which analysis and the results of Sections III-E and IV confirm this.

The rest of this paper is organized as follows: the mathematical data model is stated in Section II; the proposed method is explained in detail in Section III; simulation results are presented in Section IV; and finally, conclusions and future work are presented in Section V. Also, the details of some equations and proofs are given in the Appendix section.

Notation: Throughout the paper, superscripts (.) 1. The transmitter and receiver symmetric uniform linear arrays (ULAs) consist of M = 2M + 1 and N = 2N + 1 sensor elements, which are numbered -M, -M + 1, ..., 0, ..., M and -N, -N + 1, ..., 0, ..., N , respectively. The distance of adjacent sensors in each ULA is d. We use symmetric arrays to be able to simplify equations to the form required in Section III. Therefore, for arrays to be uniform (equal inter-element spacings), the number of elements in both the transmitting and receiving arrays must be odd. In the case of an array equipped with even-numbered antennas, one of the end antennas can be easily switched off, or the corresponding information ignored. In this case, the phase reference will still be in the position of one of the intermediate elements. M orthogonal waveforms are transmitted in the MIMO radar system, and the waveform in the m-th transmitting antenna is denoted as s m (t, l) in the time domain, where l denotes the pulse index, the number of pulses is L, m = -M, ..., 0, ..., M and l = 1, 2, ..., L. Therefore, we have

t∈T l s m (t, l) s * m ′ (t, l) dt = 1, m = m ′ 0, m ̸ = m ′ , (1)
where T l is the pulse duration and m ′ = -M, ..., 0, ..., M . Assuming that there are K NF targets, the DOA and range for the k-th target (k = 1, 2, ..., K) are denoted as θ k and r k , respectively. The target backscatter coefficients (ρ k (l) ∈ C) are assumed to obey the Swerling II model, i.e., they remain constant during the duration of one radar pulse but change from pulse to pulse [START_REF] Cao | Doa estimation for multiple targets in mimo radar with nonorthogonal signals[END_REF][START_REF] Khabbazibasmenj | Efficient transmit beamspace design for search-free based doa estimation in mimo radar[END_REF][START_REF] Mao | New approach for doa estimation in mimo radar with nonorthogonal waveforms[END_REF]. Also, they follow the independent and identical distribution between pulses. Moreover, without the loss of generality, they are assumed zero-mean. Later in Section III-D, the bias compensation for the non-zero mean will be discussed. The received signals in the n-th receiving antenna can be expressed as

y n (t, l) = K k=1 M m=-M s m (t, l) e jT m, k ρ k (l) e jT n, k +ω n (t, l) , (2)
where ω n (t, l) represents the additive white Gaussian noise . in the n-th receiving antenna during the l-th pulse, and n = -N, ..., 0, ..., N . T n, k denotes the phase shift associated with the k-th target due to the propagation time delay between the reference sensor and the n-th sensor. By considering the array center as the phase reference, T n, k can be obtained by [START_REF] Molaei | A pure cumulant-based method with low computational complexity for classification and localization of multiple near and far field sources using a symmetric array[END_REF] T

n, k = 2π λ r k 1 + (nd/r k ) 2 -2nd sin θ k /r k -1 , (3)
where λ is the signal wavelength. According to the secondorder Taylor series expansion, in the Fresnel region, which is defined as 0.62 (D 3 /λ), 2D 2 /λ , where D symbolizes the array aperture size, a simplification of the model in ( 3) is used as a good approximation and given by [START_REF] Molaei | Passive localization and classification of mixed near-field and far-field sources based on highorder differencing algorithm[END_REF] T

n, k ≈ η k n + ϕ k n 2 , (4)
where the electric angles η k and ϕ k are obtained by

η k = -2π d λ sin θ k , (5)
ϕ k = π d 2 λr k cos 2 θ k . (6)
After the matched filter h m (t, l) = s * m (τ -t, l) for the waveform s m (t, l) and sampling at time τ , and considering (4), we can obtain the sampled signal as

y n, m (l) = t∈T l y n (t, l) s * m ′ (t, l) dt = K k=1 e jT n, k t M m=-M ρ k (l) s m (t, l) e jT m, k s * m ′ (t, l) dt + t ω n (t, l) s * m ′ (t, l) dt = K k=1 e j(η k m+ϕ k m 2 ) ρ k (l) e j(η k n+ϕ k n 2 ) + ω n, m (l) , (7)
where ω n, m (l) ∼ CN 0, σ 2 . By collecting all the signals in the n-th receiving antenna, the vector form of the received signal can be obtained as

y n (l) ∆ = [y n, -M (l) , ..., y n, 0 (l) , ..., y n, M (l)] T = K k=1 ρ k (l) b n (θ k , r k ) a (θ k , r k ) + ω n (l) ∈ C M ×1 , (8)
where the steering vector in the transmitter and the noise vector, respectively, are defined as

a (θ k , r k ) ∆ = [a -M (θ k , r k ) , ..., a 0 (θ k , r k ) , ..., a M (θ k , r k )] T ∈ C M ×1 , (9)
ω n (l) ∆ = [ω n, -M (l) , ..., ω n, 0 (l) , ..., ω n, M (l)] T ∈ C M ×1 , (10)
where a m (θ k , r k ) = e j(η k m+ϕ k m 2 ) . By collecting all the received signals into a matrix, we can obtain

Y (l) ∆ = y T -N (l) , ..., y T 0 (l) , ..., y T N (l) T = K k=1 ρ k (l) b (θ k , r k ) ⊗ a (θ k , r k ) + W (l) = Dρ (l) + W (l) ∈ C N M ×1 , (11)
where the steering vector in the receiver, the noise matrix, the steering matrix and the vector for target scattering coefficients, respectively are defined as

b (θ k , r k ) ∆ = [b -N (θ k , r k ) , ..., b 0 (θ k , r k ) , ..., b N (θ k , r k )] T ∈ C N ×1 , (12)
W (l) ∆ = ω T -N (l) , ..., ω T 0 (l) , ..., ω T N (l) T ∈ C N M ×1 , (13)
D ∆ = [b (θ 1 , r 1 ) ⊗ a (θ 1 , r 1 ) , b (θ 2 , r 2 ) ⊗ a (θ 2 , r 2 ) , ..., b (θ K , r K ) ⊗ a (θ K , r K )] ∈ C N M ×K , (14)
ρ (l) ∆ = [ρ 1 (l) , ρ 2 (l) , ..., ρ K (l)] T ∈ C K×1 . ( 15
)
where b n (θ k , r k ) = e j(η k n+ϕ k n 2 ) .

III. PROPOSED METHOD In this section, the details of the mixed-order statistics algorithm (MOSA) are presented in the form of two methods. First, in Section III-A, we define two spatial expressions based on FOS and SOS from the received data. These expressions are defined mathematically so as to provide the information needed to estimate the DOA and range parameters in Sections III-B and III-C, respectively. Also, in Section III-D, the information of both expressions is used to estimate the reflectivity coefficients. Section III-E presents the main steps in implementing the algorithms.

A. Definition of the Special Spatial Expressions

We define one special spatial cross-cumulant c 4y (u, v, p, q) and one special spatial cross-covariance r y (u, v, p, q) for the array output signals with different sensor lags in the following forms:

c 4y (u, v, p, q) ∆ = cum 4 y * u, p (l) , y v, q (l) , y * -v, -q (l) , y -u, -p (l) , (16)
r y (u, v, p, q) ∆ = cov y * v, q (l) , y u, p (l) , (17)
where

u, v ∈ [-N, N ] and p, q ∈ [-M, M ].
The cumulant and covariance of ( 16) and ( 17) are calculated as

. c 4y (u, v, p, q) = E y * u, p (l) y v, q (l) y * -v, -q (l) y -u, -p (l) -E y * u, p (l) y v, q (l) E y * -v, -q (l) y -u, -p (l) -E y * u, p (l) y * -v, -q (l) E {y v, q (l) y -u, -p (l)} -E y * u, p (l) y -u, -p (l) E y v, q (l) y * -v, -q (l) , (18)
r y (u, v, p, q) = E y * v, q (l) , y u, p (l) . ( 19
)
The results of the cumulant and covariance of ( 18) and ( 19) can be written as (for more details refer to Appendix A)

c 4y (u, v, p, q) = K k=1 c 4ρ k e j2η k (-u) e -j2η k (-v) e j2η k (-p) e -j2η k (-q) , (20)
r y (u, v, p, q) = K k=1 r ρ k e j(η k u+ϕ k u 2 ) e -j(η k v+ϕ k v 2 ) e j(η k p+ϕ k p 2 ) e -j(η k q+ϕ k q 2 ) + σ 2 δ (u -v) δ (p -q) , (21)
where c 4ρ k and r ρ k are the kurtosis and the covariance related to the k-th backscatter coefficient (reflectivity), respectively, and can be written as (for more details refer to Appendix A)

c 4ρ k = cum 4 {ρ * k (l) , ρ k (l) , ρ * k (l) , ρ k (l)} = E |ρ k (l)| 4 -2E 2 |ρ k (l)| 2 -E ρ 2 k (l) E ρ * k 2 (l) , (22)
r ρ k = cov {ρ * k (l) , ρ k (l)} = E |ρ k (l)| 2 . ( 23
)
Note that the cumulants defined in [START_REF] Molaei | Components separation algorithm for localization and classification of mixed nearfield and far-field sources in multipath propagation[END_REF] and [START_REF] Zheng | Classification and localization of mixed near-field and far-field sources using mixed-order statistics[END_REF] are based on a passive scenario, and the interactions between indexes are limited to elements of an array (receiver array); whereas in ( 18)-( 21), the interactions between the indices include the elements of the receiver and transmitter arrays. Also, the approach of selecting the indices in this paper is such that the output of the first expression lacks range information; while the cumulant output in [START_REF] Zheng | Classification and localization of mixed near-field and far-field sources using mixed-order statistics[END_REF] includes both DOA and range parameters, and in [START_REF] Molaei | Components separation algorithm for localization and classification of mixed nearfield and far-field sources in multipath propagation[END_REF], the definition of cumulants is done with the aim of separating non-coherent and coherent components [START_REF] Molaei | Efficient clustering of non-coherent and coherent components regardless of sources' powers for 2d doa estimation[END_REF]. So, they have different definitions and outputs. In addition, since this paper considers an active scenario, reflectivity coefficients have also appeared in the outputs; while in [START_REF] Molaei | Components separation algorithm for localization and classification of mixed nearfield and far-field sources in multipath propagation[END_REF] and [START_REF] Zheng | Classification and localization of mixed near-field and far-field sources using mixed-order statistics[END_REF], the output of spatial expressions includes the signals received from non-cooperative sources of illumination in the environment.

B. DOA Estimation

According to [START_REF] Chen | A mimo radar-based doa estimation structure using compressive measurements[END_REF], and by collecting all the sensor lags, we can construct the cross-cumulant matrix of the sensors in the following form:

C ∆ =      C -M, -M C -M, -M +1 ... C -M, M C -M +1, -M C -M +1, -M +1 ... C -M +1, M . . . . . . . . . . . . C M, -M C M, -M +1 ... C M, M      = A v1 C ρ A H v1 ∈ C M N × M N , (24)
where the (i, i ′ )-th entry of the matrix

C m, m ′ ∈ C N × N is equal to c 4y (i -N -1, i ′ -N -1, m, m ′ ), and i, i ′ ∈ [1, 2N + 1].
In this way, the virtual steering matrix A v1 ∈ C M N ×K and the matrix C ρ ∈ C K×K can be written as follows:

A v1 = [a v1 (θ 1 ) ⊗ b v1 (θ 1 ) , a v1 (θ 2 ) ⊗ b v1 (θ 2 ) , ..., a v1 (θ K ) ⊗ b v1 (θ K )] , (25)
C ρ = diag [c 4ρ1 , c 4ρ2 , ..., c 4ρ K ] , (26)
where the virtual steering vectors

a v1 (θ k ) ∈ C M ×1 and b v1 (θ k ) ∈ C N ×1 are obtained as a v1 (θ k ) ∆ = [ã -M (θ k ) , ..., ã0 (θ k ) , ..., ãM (θ k )] T , ( 27
) b v1 (θ k ) ∆ = b-N (θ k ) , ..., b0 (θ k ) , ..., bN (θ k ) T , (28)
where ãm (θ k ) = e -j2η k m and bn (θ k ) = e -j2η k n . As can be seen, C lacks range information. In fact, in the definition of cumulant c 4y , we considered the index of the sensor elements in such a way that the output of cumulant does not depend on the range. This trick will help us avoid complex multidimensional estimates with high computational costs.

To estimate the DOAs of the targets, we implement the eigenvalue decomposition (EVD) of the Hermitian matrix C as

C = QΛQ H = Q s Λ s Q H s + Q n Λ n Q H n = [q 1 , q 2 , ..., q K ] diag [χ 1 , χ 2 , ..., χ K ] [q 1 , q 2 , ..., q K ] H + [q K+1 , q K+2 , ..., q M N ] diag [χ K+1 , χ K+2 , ..., χ M N ] [q K+1 , q K+2 , ..., q M N ] H ( 29
)
where Λ is a diagonal matrix with eigenvalues arranged as

|χ 1 | ≥ ... ≥ |χ K | > |χ K+1 | ≥ ... ≥ |χ M N |, Q ∈ C M × N is the matrix of eigenvectors with column vectors q i ′′ (1 ≤ i ′′ ≤ M N ). Q s ∈ C M N ×K , spanning the signal subspace of C, consists of the eigenvectors related to the diagonal elements of Λ s ∈ R K×K . Similarly, Q n ∈ C M N ×( M N -K) consists of the eigenvectors related to the diagonal elements of Λ n ∈ R ( M N -K)×( M N -K) , which spans the noise subspace of C. Since Q n is orthogonal to a v1 (θ) ⊗ b v1 (θ)
, with respect to [START_REF] Zhi | Near-field source localization via symmetric subarrays[END_REF] and ( 27)-( 29), we can estimate the DOAs of targets by the 1D spectral search function as

θk = arg max θ (a v1 (θ) ⊗ b v1 (θ)) H Q n Q H n (a v1 (θ) ⊗ b v1 (θ))] -1 , k = 1, ..., K. (30)
Since the EVD operation has a relatively high computational complexity, to improve performance, another technique for estimating DOAs without eigendecomposition is presented here. We use the PM principle [START_REF] An | Propagator method for doa estimation using fourthorder cumulant[END_REF][START_REF] Ahmed | A higher-order propagator method for 2d-doa estimation in massive mimo systems[END_REF][START_REF] Que | Joint processing of doa estimation and signal separation for planar array using fast-parafac decomposition[END_REF][START_REF] Dong | Robust doa estimation for sources with known waveforms against doppler shifts via oblique projection[END_REF] as an alternative to the EVD approach to estimate the noise subspace with less computational complexity. C can be partitioned according to [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] as follows:

C = C ′ C ′′ = A ′ v1 A ′′ v1 C ρ A ′ v1 A ′′ v1 H , (31)
.

where

C ′ ∈ C α× M N , C ′′ ∈ C ( M N -α)× M N , A ′ v1 ∈ C α×K , A ′′ v2 ∈ C ( M N -α)×K
, and α is an arbitrary value of the interval 1, M N -1 . Thus, there will be a matrix Π ∈ [START_REF] Ahmed | Higher-order unitary propagator method for 2d-doa estimation of non-circular sources via uniform rectangular array[END_REF]. An estimate of Π can be obtained as follows (see [START_REF] Ahmed | A higher-order propagator method for 2d-doa estimation in massive mimo systems[END_REF][START_REF] Que | Joint processing of doa estimation and signal separation for planar array using fast-parafac decomposition[END_REF][START_REF] Dong | Robust doa estimation for sources with known waveforms against doppler shifts via oblique projection[END_REF] for details):

C K×( M N -α) such that A ′′ v1 = Π H A ′ v1
Π = C ′ C ′ H -1 C ′ C ′′ H . ( 32
)
If we define

F ∆ = Π T -I M N -α T , we have F H A v1 = Π H A ′ v1 -A ′′ v1 = 0.
In other words, the columns of matrix F form a basis of the null space of A v1 . This means that the subspace spanned by the columns of F is equivalent to the subspace spanned by the noise subspace obtained through an eigendecomposition of C. So, we can write

A H v1 F = 0. ( 33
)
To ensure more accurate DOA estimates and improve resolution, we orthonormalize [START_REF] Akkar | Improved reactance domain unitary propagator algorithms for electronically steerable parasitic array radiator antennas[END_REF] the components of the matrix F as follows:

Fo = F FH F -1 2 . ( 34
)
Finally, the DOAs of all signals can be estimated by the 1D spectral search function as:

θk = arg max θ (a v1 (θ) ⊗ b v1 (θ)) H Fo FH o (a v1 (θ) ⊗ b v1 (θ))] -1 , k = 1, ..., K. (35
) Remark 1: In the proposed method, to avoid phase ambiguity, λ/4 is considered as the upper bound of d (see Appendix B for more details).

C. Range Estimation

According to [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF], and by collecting all the sensor lags, we can construct the cross-covariance matrix of the sensors as

R ∆ =      R -M, -M R -M, -M +1 ... R -M, M R -M +1, -M R -M +1, -M +1 ... R -M +1, M . . . . . . . . . . . . R M, -M R M, -M +1 ... R M, M      = A v2 R ρ A H v2 ∈ C M N × M N , (36)
where the

(i, i ′ )-th element of the matrix R m, m ′ ∈ C N × N is equal to r y (i -N -1, i ′ -N -1, m, m ′ )
, and the virtual steering matrix A v2 ∈ C M N ×K and the matrix R ρ ∈ C K×K can be written as

A v2 = [a (θ 1 , r 1 ) ⊗ b (θ 1 , r 1 ) , a (θ 2 , r 2 ) ⊗ b (θ 2 , r 2 ) , ..., a (θ K , r K ) ⊗ b (θ K , r K )] , (37)
R ρ = diag [r ρ1 , r ρ2 , ..., r ρ K ] . (38)
To estimate the ranges of the targets, EVD is applied to R

R = UΨU H = U s Ψ s U H s + U n Ψ n U H n = [u 1 , ..., u K ] diag [ψ 1 , ..., ψ K ] [u 1 , ..., u K ] H + [u K+1 , ..., u M N ] diag [ψ K+1 , ..., ψ M N ] [u K+1 , ..., u M N ] H , (39)
Now, we define the function f (θ, r) as follows:

f (θ, r) ∆ = (a (θ, r) ⊗ b (θ, r)) H U n U H n (a (θ, r) ⊗ b (θ, r)) -1 . ( 40
)
With regard to ( 9) and ( 12), and using the following equation, we can estimate the range of targets:

rk = arg max r f θk , r , k = 1, ..., K. (41)
These 1D spectral searches are performed within the Fresnel region. In the above method, the corresponding DOAs and ranges are automatically paired.

D. Reflectivity Estimation

With the DOA and range pairs, θk , rk s, estimated in the previous sections, an estimate of the steering matrix can be reconstructed as follows:

D ∆ = b θ1 , r1 ⊗ a θ1 , r1 , ..., b θK , rK ⊗ a θK , rK . (42)
Therefore, by considering [START_REF] Zheng | Classification and localization of mixed near-field and far-field sources using mixed-order statistics[END_REF] and that the mean noise is zero, an estimate of the backscatter coefficients can be obtained in the following form:

ρ = D † Ȳ, (43)
where

Ȳ = L l=1
Y (l)/L . Although the mean obtained from ( 43) is sufficient to solve our problem, ( 44) and ( 45) also provide estimates of the covariance and kurtosis of the backscatter coefficients, respectively (see Appendix C for proof):

rρ k = a θk , rk ⊗ b θk , rk H U s Ψ s -σ2 I K U H s † a θk , rk ⊗ b θk , rk -1 , (44)
ĉ4ρ k = a v1 θk ⊗ b v1 θk H Q s Λ s Q H s † a v1 θk ⊗ b v1 θk -1 . ( 45
)
The above results are obtained by assuming that the distribution's mean of backscatter coefficients is zero. In the more general case, where this mean is not necessarily zero, the mean bias can be eliminated from the received data with a simple subtraction Ỹ = Y-Ȳ before processing. This action is valid because Dρ + W -E {Dρ + W} = D (ρ -E {ρ}) + W. In this case, (43) can be updated as follows:

ρ = D † Ỹ + Ȳ . (46)
Although [START_REF] Liu | Spatial differencing method for mixed far-field and near-field sources localization[END_REF] and ( 45) are valid for both complex and real backscatter coefficients, in the case of real coefficients, by simplifying ( 22) and ( 23), estimates of the second-order moment (SOM) and the fourth central moment (FCM) of the backscatter coefficients can be obtained from rρ k + ρ2 k and ĉ4ρ k + 3r 2 ρ k , respectively.

.

E. Steps of Implementing the Proposed Method

The main steps of implementing the proposed algorithms MOSA I and MOSA II can be found in Algorithm 1. Furthermore, their computational complexity is computed considering the major multiplications involved in statistical matrices construction, eigen-decomposition and spectral search [START_REF] Molaei | High-performance 2d doa estimation and 3d localization for mixed near/far-field sources using fourthorder spatiotemporal algorithm[END_REF]. Note that only steps 3 and 4 are different in the two methods and the other steps are the same in both algorithms.

Algorithm 1

The main steps of implementing the proposed algorithms MOSA I and MOSA II.

INPUT: Y OUTPUT: θk , rk , ρk , and SOM and FCM of ρ k , k = 1, 2, .., K

Step 1:

Ȳ → Ỹ and Y -Ȳ → Y.
Step 2: Estimate C based on L pulses by using ( 18), ( 24) and [START_REF] Monti | Mantle cloaking for co-site radiofrequency antennas[END_REF].

Step 3-I: Implement the EVD of C and form Q n from (29).

Step 3-II: Partition C as [START_REF] Molaei | An efficient waveform diversity based on variational mode decomposition of coded beatfrequency shifted signals algorithm for multiple-input multiple-output millimetre-wave imaging[END_REF]. Then form the matrix Fo by using [START_REF] Molaei | Fourier compatible near-field multiple-input multiple-output terahertz imaging with sparse non-uniform apertures[END_REF] and [START_REF] Mao | New approach for doa estimation in mimo radar with nonorthogonal waveforms[END_REF].

Step 4-I: Obtain θk s from [START_REF] Marks | Fourier accelerated multistatic imaging: A fast reconstruction algorithm for multiple-inputmultiple-output radar imaging[END_REF].

Step 4-II: Obtain θk s from [START_REF] Molaei | Passive localization and classification of mixed near-field and far-field sources based on highorder differencing algorithm[END_REF].

Step 5: Estimate R based on L pulses by using ( 19), ( 36) and [START_REF] Yang | Microstrip antennas integrated with electromagnetic band-gap (ebg) structures: A low mutual coupling design for array applications[END_REF].

Step 6: Implement the EVD of R and form U n from (39).

Step 7: Obtain rk s from ( 39) and [START_REF] Dong | Robust doa estimation for sources with known waveforms against doppler shifts via oblique projection[END_REF].

Step 8: With the pairs θk , rk , obtain D according to [START_REF] Akkar | Improved reactance domain unitary propagator algorithms for electronically steerable parasitic array radiator antennas[END_REF].

Step 9: Estimate ρ k , the covariance of ρ k and kurtosis of ρ k from ρ = D † Y + Ỹ , ( 44) and [START_REF] Liu | Two-stage matrix differencing algorithm for mixed far-field and near-field sources classification and localization[END_REF], respectively.

Note that the estimation of all three parameters in the proposed method is done under the same conditions in terms of SNR, number of samples, number of elements participating in the estimation and aperture size. Therefore, if the conditions are unfavorable (for example, low SNR), it affects all estimates, whether these estimates are performed in parallel (simultaneously) or sequentially (one by one). In some approaches in the literature [START_REF] Liu | Spatial differencing method for mixed far-field and near-field sources localization[END_REF][START_REF] Liu | Two-stage matrix differencing algorithm for mixed far-field and near-field sources classification and localization[END_REF] (and not in the proposed method), the estimation of the first parameter suffers from a sharp loss in the aperture, which can affect the estimation accuracy of the next parameter. The second point is that if the first parameter (e.g. DOA) is not estimated correctly, the correct or incorrect estimation of subsequent parameters will no longer matter; because the target location is detected with both DOA and range parameters and not alone. Another point is that in the proposed method, we have used FOS, which offers higher accuracy [START_REF] Liang | Passive localization of mixed near-field and farfield sources using two-stage music algorithm[END_REF], to estimate the first parameter (DOA) to provide higher reliability.

Both methods construct M 2 FOC matrices of size N × N , and also M 2 covariance matrices of size N × N . EVD operation is applied to R. Also, MOSA I requires an additional EVD operation before estimating DOAs. Furthermore, a 1D spectral search and K ones are required to estimate the DOAs and ranges, respectively. Therefore, the major computational complexities of MOSA I and MOSA II are given in ( 47) and [START_REF] Molaei | A one-step algorithm for mixed far-field and near-field sources localization[END_REF], respectively:

O 10 M N 2 L + 8 M N 3 /3 + 180/∆ θ + K 2D 2 /λ -0.62 D 3 /λ /∆ r M N 2 , (47)
O 10 M N 2 L + 4 M N 3 /3 + 180/∆ θ + K 2D 2 /λ -0.62 D 3 /λ /∆ r M N 2 , (48)
where ∆ θ and ∆ r are the size of the angle and range search steps, respectively.

IV. SIMULATION RESULTS In this section, the results of some numerical simulations to evaluate the performance of the proposed method are presented. All results are performed on the MATLAB R2020b of 64-bit Windows 10 operating system with 12GB of random-access memory and a Core-i7 central processing unit at 2.7GHz. We compare the performance of the proposed MOSA I and MOSA II algorithms with the methods [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF], [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] and [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF] as well as the corresponding Cramer-Rao bound (CRB) [START_REF] Khamidullina | Conditional and unconditional cramér-rao bounds for near-field localization in bistatic mimo radar systems[END_REF]. ∆ θ and ∆ r are assumed to be 0.1 • and 0.1λ, respectively. The parameters L and α are assumed to be 250 and K, respectively. The control parameter in the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] is considered equal to 2 M N log M N σ n . Where not mentioned, ρ k s are assumed to be drawn from a uniform distribution in [0.5, 1.5]. In methods [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF], [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF], [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF] and the proposed method, where not mentioned, d is considered equal to d = λ/2 , d = λ/2 , d = λ/4 and d = λ/4 , respectively. The number of time lags in method [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF] is assumed to be equal to 2. To evaluate the performance of estimation accuracy, we use the root mean square error (RMSE) criterion, which is defined as follows [START_REF] Molaei | A one-step algorithm for mixed far-field and near-field sources localization[END_REF]:

RM SE = 1 N T K K k=1 N T i=1 ζk, i -ζ k 2 , (49)
where N T denotes the number of independent Monte-Carlo executions, and ζk, i is an estimate of the parameter ζ k in the i-th trial.

Example 1: In the first experiment, we examine the estimated spatial spectrums for DOA and range. Here, the values M and N are equal to 5 and 11, respectively. The interval [2.45λ, 12.5λ] approximates the Fresnel zone. Targets are located at {θ 1 = -29 • , r 1 = 11λ}, {θ 2 = 1 • , r 2 = 3λ} and {θ 3 = 41 • , r 3 = 4λ}. The SNR is set to 20dB. Figs. 234show the estimated spatial spectrums for the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF], the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] and the proposed method. Also, the computational times of the different methods are given in Table II. Regarding DOA estimates (considering Figs. 2(b) and 3), it can be seen that the proposed method has provided the best performance in terms of estimation accuracy. Given the use of FOC in the proposed method and the fact that kurtosis of Gaussian noise is zero, so by reducing the noise effect, it was expected that the estimation accuracy of the proposed algorithms would be better than the other methods. Although in [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF], in the . final stage, a 2D-MUSIC processing is used to estimate the parameters, the input data of this estimator is based on the covariance of the smoothed received signals; while in the proposed method, FOSs are used for DOA estimation. An analysis of the superiority of estimation accuracy using cumulants is given in [START_REF] Liang | Passive localization of mixed near-field and farfield sources using two-stage music algorithm[END_REF]. Although in this example, the spectra of MOSA I and MOSA II are apparently similar, the slight differences can be found by looking closely at the vertical axes. This similarity should not be interpreted to mean the same functionality, as this is just one example, and changing parameters (e.g., the size of search steps or the number of elements) may lead to different results. The method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] provides better resolution than other methods. However, by examining Table II, it is found that it is by no means suitable in terms of computational time (especially for real-time applications). The reason for the very high computational time in the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] is the use of 2D search to simultaneously estimate the angle and range. In fact, it has to search a very large space. In the proposed method, we have reduced 2D search to 1D searches, which has significantly reduced computations. Among MOSA I and MOSA II, the latter has less computational time due to the elimination of an EVD operation. The long computational time in the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] is due to solving the minimization problem. Another important drawback of the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] is the creation of spurious peaks in the estimated spectrum. This may be even more acute at lower SNRs. Fig. 3 shows the normalized DOA spectrum estimated by method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] when d decreases from λ/2 to λ/4 . As can be seen, although the spurious peaks are removed in this case, the three peaks are more difficult to detect because the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] does not provide a suitable resolution. In addition, the accuracy of the estimates has decreased. Obviously, the larger the array aperture, the more accurate the results. On the other hand, the size of the array aperture is related to the number of array elements and the distance between adjacent elements in the array. Therefore, as the number of array elements and the distance between adjacent elements increase, the accuracy of the results increases. However, according to the analyzes in Appendix B, in the proposed method, when d is considered greater than λ/4 , although the accuracy increases, phase ambiguity occurs. Therefore, for the proposed method, the distance between two adjacent elements is considered λ/4 . Figs. 2(c) and4 show the spatial spectra related to the range estimates in the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] and the proposed method, respectively. As can be seen in the case of range estimates, the proposed method still provides higher accuracy. Since the estimated DOA values in this particular example were equal for MOSA I and MOSA II, the associated spatial spectra of the range are the same, as the range estimation mechanism is the same for both algorithms. Note that the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] provides only the DOA parameter estimation. However, since our proposed method estimates DOA separately from other parameters, the virtual steering vectors used in our methods contain pure angles (see [START_REF] Cui | Joint estimation of near-field source parameters and array response[END_REF] and ( 28)), similar to the steering vectors used in the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] (although the forms of these vectors are different in these methods).

Example 2: In the second scenario, we reduce the number of transmitting and receiving antennas to 3 and 9, re- -7 show the estimated spatial spectrums for the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF], the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] and the proposed method. Also in Table III, the computational times of different methods are given. It is clear that by reducing the number of transmitting and receiving elements, it is expected that the computational time in all methods will be reduced, which is confirmed by comparing Tables II andIII. By considering Figs. 567and Table III, it can be concluded that by reducing the number of elements, the proposed method still provides the best performance in terms of estimation accuracy and computational time. suffers from the problem of creating spurious peaks and low resolution. The reason for the reduction in computational time compared to the previous example is the reduction in the size of the processed data due to the reduction in the number of elements. Also, according to Table III, it can be seen that the computational time of the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] is more than 25 times the proposed method. Note that since the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF], as well as the proposed method, are based on subspace decomposition, they provide better resolution.

Example 3: In the third example, we calculate the RMSEs . As can be seen, at low SNRs, the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] provides the best accuracy. However, as the SNR increases, the superiority of the proposed method becomes more apparent, so that for SNRs higher than 12dB, the closest diagrams to the CRBs belong to the proposed method. At very low SNRs, the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] shows the worst performance. A major reason may be the negative impact of spurious peaks created in this method, which was mentioned in previous examples. By Comparing Figs. 8 and9, it can be seen that in general, CRB and MOSA II have shown the lowest and highest sensitivity to the change in the number of elements, respectively (mean point-to-point changes in algorithms [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF], [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF], MOSA I, MOSA II, and CRB It can also be seen that the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] shows a saturation behavior so that by increasing the SNR, no significant improvement in accuracy is achieved. This is because it requires a 2D spatial smoothing process on the spatial/time sampled data matrix. An insufficient number of time samples can lead to such saturation behavior. An analysis of this has been done in [START_REF] Liu | Localization of nearfield signals based on linear prediction and oblique projection operator[END_REF]. Note that although the RMSEs in Fig. 11 are calculated with a smaller number of elements than in Fig. 10 (and we expect larger RMSEs), the area of the range search, as well as the corresponding values, are smaller in the latter case (we expect smaller RMSEs). Therefore, given these two opposite cases, it is expected that the variance of the errors in the two cases will moderate each other to some extent, which the results confirm. Note that CRB represents an optimal and theory performance that is independent of the type of algorithm and depends only on the simulation parameters [START_REF] Molaei | Derivation of multidimensional stochastic cramer-rao bound for mixture of noncoherent and coherent near/far field signals[END_REF]. However, estimation algorithms can eventually exhibit a suboptimal behavior, in which the RMSE depends on the problem-solving method in addition to the value of the parameters. Example 4: In the fourth example, we examine the performance of the proposed algorithms in estimating the reflectivity coefficients (which have been neglected in similar works). Figs. 12(a the reflectivity, SOM and FCM of the backscatter coefficients in 100 independent experiments by MOSA I and MOSA II with the parameters of Example 1, respectively. As can be seen, although the mean of the target backscatter coefficients is non-zero, both methods have successfully estimated all three parameters. These estimates are also obtained with the parameters of Example 2, the results of which are given in Figs. 13(a) and 13(b). Figs. 12 and 13 confirm that the estimates in 100 independent experiments (marked with circular markers) are distributed around the correct values (lines). Table IV shows the RMSEs calculated in Figs. 12 and13. Obviously, with decreasing number of elements, accuracy is expected to decrease, which can be seen by comparing Figs. 12 and 13 as well as by examining Table IV. As can be seen, the estimated samples in Figs. 13(a) and 13(b) show more deviation than the corresponding samples in Figs. 12(a) and 12(b). Also, the results estimated by MOSA I are somewhat more accurate than MOSA II, which is consistent with the findings of previous examples. In order to show the compatibility of the proposed method with the complex backscatter coefficients, here we present the results of the estimates assuming the coefficients are complex and assuming the same distribution. targets. Now, similar to Example 2, we reduce the number of elements and repeat the above experiment. It can be seen that this time the angular resolution decreases slightly; so that all three peaks can be identified up to the angle of 9 • (purple diagram), but when the third target is at the angle of 8 • , the threshold of joining the two adjacent peaks can be identified (see green diagram). This relative decrease in performance is consistent with a decrease in the number of elements.

Example 6: In the sixth example, an analysis of computational complexity is presented. Fig. 17 shows a comparison between the computational complexity in the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] and the proposed method versus the number of elements for ∆ θ equal to 0.01 • , 0.1 • and 1 • . The computational complexities of the proposed method are given in [START_REF] Khamidullina | Conditional and unconditional cramér-rao bounds for near-field localization in bistatic mimo radar systems[END_REF] and [START_REF] Molaei | A one-step algorithm for mixed far-field and near-field sources localization[END_REF]. The major computational complexity in the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] is equal to

O l 1 l 2 M N -l 1 + 1 (L -l 2 + 1) L + 4(l 1 l 2 ) 3 /3 +180 2D 2 /λ -0.62 D 3 /λ (l 1 l 2 ) 2 /(∆ θ ∆ r ) , (50)
where l 1 × l 2 denotes the dimensions of the scanning window [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF]. It is clear that as the value of ∆ θ increases, the complexity decreases; because the search steps get bigger. This comes at the cost of reducing the estimation accuracy. As Fig. 17 shows, the computational complexities of the proposed method are much less than that of the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] (a) (b) (especially for a small or medium number of elements). This advantage holds for all three values of ∆ θ . This confirms the results of Tables II andIII. Also, as seen, MOSA II has less overall computational complexity. The dominant term in computational complexity is the spectral search, which is the same in both algorithms and varies with M N 2 . But the term related to eigendecomposition changes with M N 3 , so as the number of elements increases, the two graphs are expected to be further apart. A general conclusion and comparison of the performance of the two proposed algorithms are provided in Section V. proposed method have both higher estimation accuracy and much less scatter around the true location, which indicates a more reliable MIMO structure compared to a SIMO one. This was not unexpected because a MIMO system, by using spatial diversity, can provide better performance than a conventional SIMO structure. In the case of MOSA I and MOSA II, as expected, the results are almost identical, which again confirms the practicality of using the PM principle instead of eigendecomposition to estimate the noise subspace. The slightly better estimation results in MOSA I than in MOSA II are consistent with the results obtained in the previous examples. Also, as can be seen, the standard deviation values of the DOA estimates for all three targets in each method are close to each other. In such cases, according to the analyzes presented in [START_REF] Yuen | Performance analysis of higher order esprit for localization of near-field sources[END_REF], it is expected that the values of the standard deviation of the first target will be higher than those of the other targets because it has a greater range. The results presented in Fig. 18 curacy, causes insensitivity to Gaussian noises. The proposed method, by separating the parameter estimation procedure into different phases, eliminates the very heavy multidimensional search. This will be very effective in significantly reducing the computational effort. Moreover, the estimated parameters are paired automatically and no additional steps are required.

The performance of the proposed method was compared with similar methods as well as theoretical CRBs. The simulation results showed that the proposed method has much less computational time than similar methods, which is a major advantage for employing them in real-time applications.

It was also shown that in terms of estimation accuracy of all parameters, they provide high performance (especially for medium and high SNRs). Although the resolution of the proposed method is better than the method [START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF], the results showed that the method [START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] offers better performance in terms of angular resolution than other methods.

All in all, given all the examples, one can draw a general conclusion about MOSA I and MOSA II, that MOSA I as a whole provides a relatively better estimate of the parameters; however, it has relatively more computational and time complexities. Increasing the number of elements brings the accuracy of these two methods closer to each other, but the difference in calculations in them becomes more. Therefore, in low-element and high-element scenarios (massive MIMO scenarios), the use of MOSA I and MOSA II, respectively, is recommended. This is a logical trade-off that can be applied in different applications according to priority and need.

To put the proposed method into further applications, future works will be concentrated on developing the proposed approach for more complex scenarios and providing new solutions to the problem, as follows:

■ Dealing with a multipath environment [START_REF] Molaei | Closed-form expression of stochastic crb for mixed near-field and far-field sources in multipath propagation environments[END_REF]: In this scenario, the received signals are no longer necessarily independent [START_REF] Molaei | Two-dimensional doa estimation for multi-path environments by accurate separation of signals using k-medoids clustering[END_REF]. Therefore, not only will the simplifications in Appendix A no longer be valid (the equations will be much more complex), but the rank of the statistical matrices will decrease [START_REF] Molaei | Components separation algorithm for localization and classification of mixed nearfield and far-field sources in multipath propagation[END_REF]. The use of radio-based simultaneous localization and mapping process [START_REF] Wymeersch | Radio localization and mapping with reconfigurable intelligent surfaces: Challenges, opportunities, and research directions[END_REF][START_REF] Del Hougne | Deeply subwavelength localization with reverberation-coded-aperture[END_REF] and reconfigurable intelligent surfaces technology [START_REF] Wymeersch | Radio localization and mapping with reconfigurable intelligent surfaces: Challenges, opportunities, and research directions[END_REF][START_REF] Del Hougne | Robust position sensing with wave fingerprints in dynamic complex propagation environments[END_REF] may provide efficient solutions for more complex assumptions and more specific applications in this direction.

■ Considering the mutual coupling effect [START_REF] Zheng | Misc array: A new sparse array design achieving increased degrees of freedom and reduced mutual coupling effect[END_REF]: In general, there are two strategies to deal with the issue of mutual coupling:

• In the physical layer: The use of special antennas made or insulated with special materials, for example, using mantle cloaking [START_REF] Monti | Mantle cloaking for co-site radiofrequency antennas[END_REF]. The use of photonic bandgap structure between elements [START_REF] Yang | Microstrip antennas integrated with electromagnetic band-gap (ebg) structures: A low mutual coupling design for array applications[END_REF] and neutralization line [START_REF] Mansouri | Design and evaluation of a compact diversity antenna for lte application[END_REF] are other solutions. On the other hand, electrically small antennas [START_REF] Patel | A discussion on electrically small antenna property[END_REF] with a size of λ/10 or less can be used. In this case, the distance between the antennas can be less than λ/10 ; and to reduce mutual coupling in this case, metamaterials with cells less than λ/100 can be used.

• Another strategy could be to use a compensator matrix in the receiver's processing section. This method is a suitable technique to eliminate the effect of mutual coupling in dipole array antennas. By considering the effect of mutual coupling between the array elements, [START_REF] Zheng | Classification and localization of mixed near-field and far-field sources using mixed-order statistics[END_REF] changes to the following form:

Y (l) = MDρ (l) + W (l) , (51)
where M denotes the mutual coupling matrix. As discussed in [START_REF] Mao | Doa estimation of coherent signals based on direct data domain under unknown mutual coupling[END_REF], it is often sufficient to consider the ULA coupling model with finite non-zero coefficients, and a symmetric Toeplitz matrix can be used to model mutual coupling. Therefore, the output observed at the n-th sensor can be written as follows:

. y n (l) = 0 1×(n-P ) , m P -1 , ..., m 1 , 1, m 1 , ..., m P -1 ,

0 1×( N -n-P +1) Dρ (l) + ω n (l) , (52)
where when the distance between the two sensors is greater than P d, the mutual coupling coefficient is assumed to be zero [START_REF] Wang | Doa estimation under unknown mutual coupling and multipath with improved effective array aperture[END_REF]. Mutual coupling effect will exist even if the distance between the elements is λ/2 . Thus, this issue is usually addressed as an independent problem in the literature. On the other hand, the inter-element spacing of λ/4 is known as conventional spacing in some localization problems, especially in the NF [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF][START_REF] Wan | A fast group sparsity based phase retrieval algorithm for non-coherent doa estimation[END_REF][START_REF] Ma | Mixed far-field and near-field source localization using a linear tripole array[END_REF]. The development of the proposed method considering mutual coupling will be investigated in future works. ■ In this paper, a spectral search-based estimator was employed. Although the proposed method provides less computational complexity than multidimensional estimation techniques and optimization solvers, the use of ESPRIT-like estimators can further reduce computational complexity. The ESPRIT technique requires a structure with two shift-invariant subarrays from an antenna array to form two basic ESPRIT equations [START_REF] Jo | Direction of arrival estimation using nonsingular spherical esprit[END_REF]. Although the separation of terms corresponding to DOA and range in this paper helps to bring the steering matrix form closer to the basic ESPRIT equations [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF], since we are dealing with a MIMO (not SIMO) system, the structure of the virtual steering matrix is more complex. Specifically in our case, A v1 lacks the form required to apply the ESPRIT technique directly. The k-th column of A v1 is given below:

a v1 (θ k ) ⊗ b v1 (θ k ) = [e -j2η k (-M -N ) , e -j2η k (-M -N +1) → ×e -j2η k , e -j2η k (-M -N +2) , . . . , e -j2η k (-M +N ) , e -j2η k (-M +1-N ) → ×e -j2η k (1-2N ) , e -j2η k (-M +1-N +1) , . . . , e -j2η k (M +N ) ] T (53)
As can be seen, the translational invariance requirement for using the ESPRIT technique does not exist in the virtual steering vector. In other words, it is not possible to directly create the two shift-invariant vectors needed to form the two basic ESPRIT equations. Note that in the above equation only a small part of the k-th column of A v1 is given. In fact, shifts e -j2η k (1-2N ) are repeated alternately along the vector. The limitation mentioned above does not apply to the MUSIC estimator; because it is compatible with any type of array geometry (or steering vector). To use the ESPRIT technique, it is necessary to redefine the spatial expressions in ( 16) and [START_REF] Hayashi | Doa estimation in mimo radar using temporal spatial virtual array with music algorithm[END_REF] and/or the spatial matrices in [START_REF] Zhi | Near-field source localization via symmetric subarrays[END_REF] and [START_REF] Molaei | Efficient clustering of non-coherent and coherent components regardless of sources' powers for 2d doa estimation[END_REF]. We intend to consider this as one of the future works. Although there is a technical justification for not using ESPRIT in this work, it should be noted that under equal conditions, the aperture loss in ESPRIT-based techniques is slightly greater than in spectral search-based estimators [START_REF] Molaei | A pure cumulant-based method with low computational complexity for classification and localization of multiple near and far field sources using a symmetric array[END_REF]. Also, despite the less complexity of ESPRIT, MUSIC provides higher accuracy [START_REF] Liu | Off-grid doa estimation for colocated mimo radar via reduced-complexity sparse bayesian learning[END_REF].

APPENDIX A Proof of (20) and (21) and their estimation: According to (7) and starting from [START_REF] Lonkeng | Two-dimensional doa estimation using arbitrary arrays for massive mimo systems[END_REF], we can obtain [START_REF] Wymeersch | Radio localization and mapping with reconfigurable intelligent surfaces: Challenges, opportunities, and research directions[END_REF]. Due to the independence of reflectivity coefficients from each other as well as noise, (54) can be simplified as [START_REF] Del Hougne | Deeply subwavelength localization with reverberation-coded-aperture[END_REF]. Due to the noise being Gaussian, the terms related to the noise neutralize each other. Finally, by simplifying exponential terms and factoring the common expression e j2η k (-u+v-p+q) in the terms related to the signal, the following is obtained:

c 4y (u, v, p, q) = K i=1 [E {ρ * k (l) ρ k (l) ρ * k (l) ρ k (l)} -E {ρ * k (l) ρ k (l)} E {ρ * k (l) ρ k (l)} -E {ρ * k (l) ρ * k (l)} E {ρ k (l) ρ k (l)} -E {ρ * k (l) ρ k (l)} E {ρ k (l) ρ * k (l)}] e j2η k (-u+v-p+q) . (56)
Therefore, according to the definition of kurtosis in ( 22), ( 20) is obtained. In a similar way, and of course simpler, ( 21) can be obtained as follows:

r y = E K k=1 ρ * k (l) e -j(η k v+ϕ k v 2 ) e -j(η k q+ϕ k q 2 ) + ω * v, q (l) K k ′ =1 ρ k ′ (l) e j(η k ′ u+ϕ k ′ u 2 ) e j(η k ′ p+ϕ k ′ p 2 ) + ω u, p (l) = K i=1 E {ρ * k (l) ρ k (l)}e -j(η k v+ϕ k v 2 ) e -j(η k q+ϕ k q 2 ) e j(η k u+ϕ k u 2 ) e j(η k p+ϕ k p 2 ) + E ω * v, q (l) ω u, p (l) . (57)
Given that the number of samples is finite, c 4y (u, v, p, q) and r y (u, v, p, q) can be estimated as follows: ĉ4y (u, v, p, q) = 1 L L l=1 y * u, p (l) y v, q (l) y * -v, -q (l) y -u, -p (l)

- 1 L 2 L l
y * u, p (l) y v, q (l) L l=1 y * -v, -q (l) y -u, -p (l)

- 1 L 2 L l=1
y * u, p (l) y * -v, -q (l) L l=1 y v, q (l) y -u, -p (l)

- 1 L 2 L l=1
y * u, p (l) y -u, -p (l) L l=1 y v, q (l) y * -v, -q (l),

ry (u, v, p, q) = 1 L L l=1 y * v, q (l) y u, p (l).

APPENDIX B Upper bound of inter-element spacing: Consider ( 27), ( 28), [START_REF] Marks | Fourier accelerated multistatic imaging: A fast reconstruction algorithm for multiple-inputmultiple-output radar imaging[END_REF] and [START_REF] Molaei | Passive localization and classification of mixed near-field and far-field sources based on highorder differencing algorithm[END_REF]. From (5) we define ξ k ≜ ã (θ k ) = b (θ k ) = e j4πd sin θ k /λ . ξ k s are the constituent entries of virtual steering vectors for estimating DOAs. Without placing any conditions on d, we will have a series of ambiguous DOA estimates as follows:

θk (z) = arcsin (ν k + zλ/2d ) , k = 1, 2, ..., K, (

. c 4y (u, v, p, q) = E K k=1 ρ * k (l) e -j(η k u+ϕ k u 2 ) e -j(η k p+ϕ k p 2 ) + ω * u, p (l)

K k ′ =1
ρ k ′ (l) e j(η k ′ v+ϕ k ′ v 2 ) e j(η k ′ q+ϕ k ′ q 2 ) + ωv, q (l)

K k ′′ =1
ρ * k ′′ (l) e -j(-η k ′′ v+ϕ k ′′ v 2 ) e -j(-η k ′′ q+ϕ k ′′ q 2 ) + ω * v, q (l)

K k ′′′ =1
ρ k ′′′ (l) e j(-η k ′′′ u+ϕ k ′′′ u 2 ) e j(-η k ′′′ p+ϕ k ′′′ p 2 ) + ωu, p (l)

-E K k=1
ρ * k (l) e -j(η k u+ϕ k u 2 ) e -j(η k p+ϕ k p 2 ) + ω * u, p (l)

K k ′ =1
ρ k ′ (l) e j(η k ′ v+ϕ k ′ v 2 ) e j(η k ′ q+ϕ k ′ q 2 ) + ωv, q (l)

E K k ′′ =1
ρ * k ′′ (l) e -j(-η k ′′ v+ϕ k ′′ v 2 ) e -j(-η k ′′ q+ϕ k ′′ q 2 ) + ω * v, q (l)

K k ′′′ =1
ρ k ′′′ (l) e j(-η k ′′′ u+ϕ k ′′′ u 2 ) e j(-η k ′′′ p+ϕ k ′′′ p 2 ) + ωu, p (l)

-E K k=1
ρ * k (l) e -j(η k u+ϕ k u 2 ) e -j(η k p+ϕ k p 2 ) + ω * u, p (l)

K k ′′ =1
ρ * k ′′ (l) e -j(-η k ′′ v+ϕ k ′′ v 2 ) e -j(-η k ′′ q+ϕ k ′′ q 2 ) + ω * v, q (l)

E K k ′ =1
ρ k ′ (l) e j(η k ′ v+ϕ k ′ v 2 ) e j(η k ′ q+ϕ k ′ q 2 ) + ωv, q (l)

K k ′′′ =1
ρ k ′′′ (l) e j(-η k ′′′ u+ϕ k ′′′ u 2 ) e j(-η k ′′′ p+ϕ k ′′′ p 2 ) + ωu, p (l)

-E K k=1
ρ * k (l) e -j(η k u+ϕ k u 2 ) e -j(η k p+ϕ k p 2 ) + ω * u, p (l)

K k ′′′ =1
ρ k ′′′ (l) e j(-η k ′′′ u+ϕ k ′′′ u 2 ) e j(-η k ′′′ p+ϕ k ′′′ p 2 ) + ωu, p (l)

E K k ′ =1
ρ k ′ (l) e j(η k ′ v+ϕ k ′ v 2 ) e j(η k ′ q+ϕ k ′ q 2 ) + ωv, q (l)

K k ′′ =1
ρ * k ′′ (l) e -j(-η k ′′ v+ϕ k ′′ v 2 ) e -j(-η k ′′ q+ϕ k ′′ q 2 ) + ω * v, q (l) .

(54)

c 4y = K i=1
E {ρ * k (l) ρ k (l) ρ * k (l) ρ k (l)}e -j(η k u+ϕ k u 2 ) e -j(η k p+ϕ k p 2 ) e j(η k v+ϕ k v 2 ) e j(η k q+ϕ k q 2 ) e -j(-η k v+ϕ k v 2 ) e -j(-η k q+ϕ k q 2 ) e j(-η k u+ϕ k u 2 ) e j(-η k p+ϕ k p 2 ) -K i=1 E {ρ * k (l) ρ k (l)}e -j(η k u+ϕ k u 2 ) e -j(η k p+ϕ k p 2 ) e j(η k v+ϕ k v 2 ) e j(η k q+ϕ k q 2 ) K i=1 E {ρ * k (l) ρ k (l)} e -j(-η k v+ϕ k v 2 ) e -j(-η k q+ϕ k q 2 ) e j(-η k u+ϕ k u 2 ) e j(-η k p+ϕ k p 2 ) -K i=1 E {ρ * k (l) ρ * k (l)}e -j(η k u+ϕ k u 2 ) e -j(η k p+ϕ k p 2 ) e -j(-η k v+ϕ k v 2 ) e -j(-η k q+ϕ k q 2 ) K i=1

E {ρ k (l) ρ k (l)}e j(η k v+ϕ k v 2 ) e j(η k q+ϕ k q 2 ) e j(-η k u+ϕ k u 2 ) e j(-η k p+ϕ k p 2 ) -

K i=1 E {ρ * k (l) ρ k (l)}e -j(η k u+ϕ k u 2 )
e -j(η k p+ϕ k p 2 ) e j(-η k u+ϕ k u 2 ) e j(-η k p+ϕ k p 2 ) K i=1 E {ρ k (l) ρ * k (l)}e j(η k v+ϕ k v 2 ) e j(η k q+ϕ k q 2 ) e -j(-η k v+ϕ k v 2 ) e -j(-η k q+ϕ k q 2 ) + E ω * u, p (l) ωv, q (l) ω * v, q (l) ωu, p (l) -E ω * u, p (l) ωv, q (l) E ω * v, q (l) ωu, p (l) -E ω * u, p (l) ω * v, q (l) E {ωv, q (l) ωu, p (l)} -E ω * u, p (l) ωu, p (l) E ωv, q (l) ω * v, q (l) .

where ν k = λ∠ξ k /(4πd) and

⌈2d (-1 -ν k )/λ ⌉ ≤ z ≤ ⌊2d (1 -ν k )/λ ⌋ . (61)
For each θk and d, we have -π ≤ ∠ξ k ≤ π. To avoid any cyclic ambiguity in ∠ξ k , it must be limited to 0 ≤ ∠ξ k ≤ π and -π ≤ ∠ξ k ≤ 0 for 0 ≤ θk ≤ π/2 and -π/2 ≤ θk ≤ 0, respectively. Therefore, according to [START_REF] Patel | A discussion on electrically small antenna property[END_REF], to avoid phase ambiguity, λ/4 should be considered as the upper bound of d.

APPENDIX C

The derivation of (44) and (45): By considering A v1 , the kurtosis of the k-th backscatter coefficient can be computed as

c 4ρ k = ε H k diag [1/c 4ρ1 , 1/c 4ρ2 , ..., 1/c 4ρ K ] ε k -1 , ( 62
)
where ε k is the k-th column of the K×K-dimensional identity matrix. Due to the fact that A † v1 a v1 θk ⊗ b v1 θk = ε k , we have

c 4ρ k = (a v1 (θ k ) ⊗ b v1 (θ k )) H A † v1 H C ρ A † v1 (a v1 (θ k ) ⊗ b v1 (θ k ))) -1 . (63)
where a v1 θk ⊗ b v1 θk is the k-th column of A v1 .

Based on the property of A † v1 H = A H v1 † and due to the elimination of the noise effect in cumulant C, (45) can be reached. Equation ( 44) can be derived in a similar way. The difference is that unlike C, the sensor noise power must also be taken into account in the case of R. The amount of noise power can be estimated as follows:

σ2 = 1 M N -K M N k=K+1 ψ k . ( 64
)
APPENDIX D Given [START_REF] Chen | A new algorithm for joint range-doa-frequency estimation of near-field sources[END_REF], it is clear that the value of the covariance of the backscatter coefficient is real. Also, the first two terms in [START_REF] Molaei | High-performance localization of mixed fourth-order stationary sources based on a spatial/temporal full esprit-like method[END_REF] are real. Here we show that the third term is also real. Assume that ρ k (l) = a + jb, where a and b denote the real and imaginary parts of ρ k (l), respectively. So, we can write

. E ρ 2 k (l) E ρ * k 2 (l) = E (a + jb) 2 E (a -jb) 2 = E a 2 -E b 2 + j2E {ab} E a 2 -E b 2 -j2E {ab} = E 2 a 2 + E 2 b 2 + 4E {ab} , (65)
which is a real value. Therefore, the kurtosis of the backscatter coefficient will always be real.

Fig. 1 .

 1 Fig. 1. The system model of the co-located MIMO array for parameters estimation.

Fig. 2 .

 2 Fig. 2. The spatial spectrums for DOA and range estimations by using the method [25] ( M = 5 and N = 11); (a) the 2D spatial spectrum, (b) the spatial spectrum of DOA, (c) the spatial spectrum of range.

Fig. 3 .

 3 Fig. 3. The normalized spatial spectrum for DOA estimation by using the method [21] with d = λ/2 and d = λ/4 , and the proposed method (MOSA I and MOSA II). M = 5 and N = 11.

Fig. 4 .

 4 Fig. 4. The spatial spectrums for range estimations by using the proposed method. M = 5 and N = 11.spectively. The interval[1.75λ, 8λ] approximates the Fresnel zone. Targets are located at the same angles as Example 1 but with different ranges as {θ 1 = -29 • , r 1 = 6.1λ}, {θ 2 = 1 • , r 2 = 2λ} and {θ 3 = 41 • , r 3 = 3λ}. The SNR is still 20dB. Figs.5-7show the estimated spatial spectrums for the method[START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF], the method[START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] and the proposed method. Also in TableIII, the computational times of different methods are given. It is clear that by reducing the number of transmitting and receiving elements, it is expected that the computational time in all methods will be reduced, which is confirmed by comparing TablesII and III. By considering Figs.5-7and TableIII, it can be concluded that by reducing the number of elements, the proposed method still provides the best performance in terms of estimation accuracy and computational time. In the previous example, the difference in estimation accuracy between the proposed MOSA I and MOSA II algorithms was not very clear, but now we can see the superiority of MOSA I more clearly by examining Figs.6 and 7. Note that the ranges corresponding to the DOAs in Fig.6are marked with the same markers in Fig.7. This pairing of parameters is done completely automatically. Although the same DOAs are found in Figs.2(b) and 5(b), a closer look reveals that Fig.2(b) has a better resolution. Despite the reduction in the number of elements in this example, the very high computational time is still the major drawback of the method[START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF]. As Fig.6shows, the method[START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] still

Fig. 5 .

 5 Fig. 5. The spatial spectrums for DOA and range estimations by using the method [25] ( M = 3 and N = 9); (a) the 2D spatial spectrum, (b) the spatial spectrum of DOA, (c) the spatial spectrum of range.

Fig. 6 .

 6 Fig. 6. The normalized spatial spectrum for DOA estimation by using the method [21] with d = λ/2 and d = λ/4 , and the proposed method (MOSA I and MOSA II). M = 3 and N = 9.

Fig. 7 .

 7 Fig. 7. The spatial spectrums for range estimations by using the proposed method (MOSA I and MOSA II). M = 3 and N = 9.associated with the estimates in the previous examples and compare them with each other as well as with the corresponding theoretical CRBs[START_REF] Khamidullina | Conditional and unconditional cramér-rao bounds for near-field localization in bistatic mimo radar systems[END_REF]. N T is considered equal to 400. Figs.8 and 9show the RMSEs of DOA versus SNR with the assumptions of Examples 1 and 2, respectively. It is natural and obvious that with increasing SNR, the errors have decreased. As can be seen, at low SNRs, the method[START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] provides the best accuracy. However, as the SNR increases, the superiority of the proposed method becomes more apparent, so that for SNRs higher than 12dB, the closest diagrams to the CRBs belong to the proposed method. At very low SNRs, the method[START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF] shows the worst performance. A major reason may be the negative impact of spurious peaks created in this method, which was mentioned in previous examples. By Comparing Figs.8 and 9, it can be seen that in general, CRB and MOSA II have shown the lowest and highest sensitivity to the change in the number of elements, respectively (mean point-to-point changes in algorithms[START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF],[START_REF] Liu | A super-resolution doa estimation method for fast-moving targets in mimo radar[END_REF], MOSA I, MOSA II, and CRB

Fig. 8 .Fig. 9 .

 89 Fig. 8. The RMSEs of DOA versus SNR with the assumptions of Example 1.

  ) and 12(b) show the results of the estimates of .

Fig. 10 .

 10 Fig. 10. The RMSEs of range versus SNR with the assumptions of Example 1.

Fig. 11 .

 11 Fig. 11. The RMSEs of range versus SNR with the assumptions of Example 2.

Fig. 12 . 5 :Fig. 13 .

 12513 Fig. 12. Results of the estimates of the reflectivity, SOM and FCM of the backscatter coefficients in 100 independent experiments with the parameters of Example 1; (a) by MOSA I, (b) by MOSA II.TABLE IV RMSES CALCULATED IN IN FIGS. 12 AND 13

Fig. 14 .

 14 Fig. 14. Results of the estimates of the reflectivity, covariance and kurtosis of the backscatter coefficients in 100 independent experiments with the parameters of Example 1; (a) by MOSA I, (b) by MOSA II.

Example 7 :Fig. 15 .

 715 Fig. 15. Results of the estimates of the reflectivity, covariance and kurtosis of the backscatter coefficients in 100 independent experiments with the parameters of Example 2; (a) by MOSA I, (b) by MOSA II.

Fig. 16 .

 16 Fig. 16. Angular resolution results obtained by reducing the spatial spacing of the second and third targets in the proposed method (MOSA I); (a) M = 5 and N = 11, (b) M = 3 and N = 9.

Fig. 17 .

 17 Fig.[START_REF] Hayashi | Doa estimation in mimo radar using temporal spatial virtual array with music algorithm[END_REF]. Comparison of the computational complexity of the method[START_REF] Belfiori | 2d-music technique applied to a coherent fmcw mimo radar[END_REF] and the proposed method versus the number of elements for ∆ θ equal to 0.01 • , 0.1 • and 1 • .

.

Fig. 18 .

 18 Fig. 18. Estimated locations in 100 independent experiments by (a) method [22], (b) MOSA I, (c) MOSA II.

TABLE II COMPARISON

 II OF COMPUTATIONAL TIME OF DIFFERENT METHODS IN EXAMPLE 1

	Method	Computational Time
	[25]	1.08 hours
	[21]	9.18 seconds
	MOSA I	1.11 seconds
	MOSA II	0.82 seconds
	.

IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. XX, NO. XX, XXX 2021

This work was funded by the Leverhulme Trust under Research Leadership

Award RL-2019-019.